Abstract
In this paper, computer vision enables recommending a reduced order model for fast stress prediction according to various possible loading environments. This approach is applied on a macroscopic part by using a digital image of a mechanical test. We propose a hybrid approach that simultaneously exploits a data-driven model and a physics-based model, in mechanics of materials. During a machine learning stage, a classification of possible reduced order models is obtained through a clustering of loading environments by using simulation data. The recognition of the suitable reduced order model is performed via a convolutional neural network (CNN) applied to a digital image of the mechanical test. The CNN recommend a convenient mechanical model available in a dictionary of reduced order models. The output of the convolutional neural network being a model, an error estimator, is proposed to assess the accuracy of this output. This article details simple algorithmic choices that allowed a realistic mechanical modeling via computer vision.
1. Introduction
In biomechanics, computer vision and mechanical testing have been coupled to obtain patient-specific simulation approaches, as proposed in [1]. At the same time, with the growth of industry 4.0, imaging techniques are more and more widespread in factories. When combined with artificial neural networks, digital images enable the classification of products to obtain the best possible process, as proposed in [2] for olive batches classification in oil extracting process or as shown in [3] for composite materials manufacturing. Nowadays, we have the possibility of extending these methods to the classification of mechanical parts produced in industry, in order to develop part-specific decision approaches. For mechanical parts, the quality of manufacturing processes has a direct influence on the ultimate mechanical properties of the manufactured parts. For example, the way the fracture is initiated in a specimen often reveals defects in the material whose origin can be tracked back to the manufacturing process [4]. In general, the numerical computation of mechanical stresses in a given manufactured part allows the predictive evaluation of the link between the ultimate mechanical properties of this part and the manufacturing process. The reader can find an example of how to optimize a process for curing composite parts in [5] according to this paradigm. The mechanical modeling of manufactured parts has for purpose to verify if defects induced by a manufacturing process are tolerable, if an observed part must be rejected, or if the manufacturing process must be improved.
In this paper, we restrict our attention to the stress prediction in a part under an observed loading environment by a digital image, while including all its geometrical defects. We propose a hybrid approach that simultaneously exploits a data-driven model and a physics-based model, in mechanics of materials. The reader can find a review on hybrid modeling in [6] for remaining useful life predictions of engineering systems. We show that the strength of the proposed hybrid modeling is its ability to incorporate an error estimator related to the modeling chain with computer vision and convolutional neural networks (CNN) [7, 8].
As explained in [9], computer vision with deep convolutional neural networks has achieved state-of-the-art performance on standard recognition datasets and tasks. In this paper, we explore the capabilities of a CNN as a recommender system for the mechanical modeling of structures submitted to various loading. The proposed hybrid modeling couples a noncentered principal component analysis (PCA) and a CNN, in order to preserve an accurate description of spatial information.
Image processing for computer vision is usually very fast. It does not make sense to couple computer vision with numerical simulations of mechanical stresses that take hours of computation. Hence, we couple computer vision with reduced order modeling of structures, in order to get fast mechanical predictions of stresses.
A reduced order model is a surrogate model obtained by the projection of high dimensional equations on a reduced space, and it also involves a reduced approximation space for the variables of the high dimensional problem. When they are the same, the surrogate model is a Galerkin reduced order model [10]. In hyperreduced order models these two reduced spaces are different [11]. Then, because we consider projection of physics-based equations, hyperreduced order models preserve the physical parameters involved in the high dimensional equations. The reduced spaces involved in the hyperreduced modeling are spanned by empirical modes extracted from simulation data, by using the proper orthogonal decomposition [10] of known finite element predictions. This procedure is similar to a noncentred principal component analysis (PCA). Hence, the proposed modeling via computer vision exploits both simulation data and observational data, which are, respectively, finite element predictions and digital images of mechanical tests.
In general, the empirical modes obtained by noncentred PCA are very sensitive to the loading environment imposed when computing the simulation data. If the variety of loading conditions considered to calculate simulation data is too wide, the number of empirical modes becomes too large. They can no longer reduce the numerical complexity of the mechanical balance equations. Clustering methods have been applied for model-order reduction in [12, 13], in order to preserve small reduced-bases of empirical modes. Moreover, cluster-based reduced order modeling (CROM) has been proposed in [14] to define a small subset of critical data to learn an efficient (CROM) with a sparse approximation space [15]. In this paper, a dictionary of hyperreduced order models is generated by considering clusters of possible loading environments in the observed mechanical tests. Then, the identification of the hyperreduced order model is done via recognition of a class of mechanical loads by a convolutional neural network. In practice, each item of the dictionary is not directly a hyperreduced order model. In order to face a possible variability on the geometry of the observed structures, it is more robust to define an item of the dictionary as a set of finite element solutions for various ideal geometries and for a given class of mechanical loads. Hence, the proposed workflow is robust enough to face geometrical defects in the observed mechanical parts. We assume that the mesh for the finite element modeling of the parts is obtained by using image-meshing techniques, as proposed in [16], of segmented 3D digital images obtained by X-ray computed tomography. We refer the reader to [17, 18] for more details on finite element modeling of 3D images obtained by X-ray computed tomography. An example of X-ray computed tomography applied to manufactured parts can be found in [19]. The ideal geometries involved in the proposed workflow were obtained by using computer-aided-design (CAD). Such CAD models are usually used to find an optimal design of the manufactured parts, with parametric finite element simulations [20].
2. Materials and Methods
Prior to the stress prediction, a reduced order model is setup for the projection of the mechanical balance equations. Here the reduced order model concerns the displacement in the observed mechanical part. Usually, for the computation of reduced approximations in nonlinear problems, we formally consider all possible situations in a given parameter space [21]. The parameters aim to describe all the possible mechanical problems, in advance, before the observation of a realization of one situation. This approach defines a tensor for the description of all possible displacements. The order of this tensor is the number of scalar parameters involved in the parametric equations, plus one. For instance, if a single parameter is introduced then we need two indices i and j to have access to the value of the degree of freedom of a finite element model for the value of the parameter. When D parameters are introduced, we need D+1 indices: i, , …, to have access to a scalar value in the tensor containing all the possible displacements. This tensor formalism aims to introduce a sampling procedure of the parameter space in order to get an estimation of the reduced approximation, or a reduced basis, for the displacements. For instance, this sampling procedure can be achieved by the proper generalized decomposition (PGD) [21] or the Tensor Train decomposition [22, 23]. In this paper, the loading environment is depicted by an image of 3968 × 2976 pixels. Then the parameter space dimension is around 12 millions (the number of pixels) for the description of all possible loading environments. Hence the tensor formalism for model-order reduction would require the decomposition of a tensor of order 12 millions. To our knowledge, no tensor decomposition method has been applied to such a huge tensor order. A purely tensor approach seems to be unaffordable. In this paper, we do not pretend to model all the possible solutions of mechanical equations related to a huge parameter space. We do not follow the usual paradigm of low rank approximations. The proposed image-based modeling aims to exploit available data for fast approximate predictions with fast error estimation.
The workflow of the proposed modeling via computer vision is shown in Figure 1. Four kinds of inputs are required:(i)a 2D digital image of the part in the test machine, this image is denoted by ;(ii)a database where are saved all simulation data, ordered with respect to a cluster index and the index of the available ideal mesh in the list , these are the meshes used when generating the simulation data by the finite element method;(iii)a 3D voxel image of the part alone, we assumed that this 3D image is obtained by X-ray computed tomography;(iv)a measurement of the load magnitude at the end of the mechanical test and a measurement of the displacement at one fixed point of on the part, at the end of the test.
The CNN network aims to recognize the index β of the class of loading environment. It gives access to the simulation data required to create on the fly a hyperreduced order model. The measurement of the displacement magnitude helps to get a precise location of the load by solving an inverse problem with hyperreduced equations as proposed in [24].
The stress predicted by the hyperreduction is obtained via constitutive equations in the framework of elasticity. It depends on the spatial position in the part, denoted by , and nodal values of the displacement, denoted by the vector . This vector has components. is the order of the finite element model of the part. The superscript is introduced for the variables related to the observed mechanical part in the experimental setup. The mesh of the mechanical part is denoted by and the stress is denoted by . The reduced approximation of the displacement reads , where is the matrix form of the reduced basis (it has less columns than rows) and is the solution of the hyperreduced balance equations. Since the finite element method is a numerical scheme for partial differential equations, Dirichlet boundary conditions are applied on the boundary of the domain occupied by the mesh. Here these conditions are null. Hence the displacements belong to a vector space, which is a subspace of a Hilbert space. The columns of span a subspace of this vector space. They fulfill the Dirichlet boundary conditions applied on . The hyperreduced balance equations are set on a reduced mesh , which is the restriction of to the finite elements connected to a given list of degrees of freedom, denoted by F. The residual of the finite element equations is denoted by . Hence the hyperreduced balance equation reads [20] The larger the set F the higher the computational complexity of the projection of the equations when considering the observed geometry. In the sequel, F is the set of degrees of freedom (dof) indices near the loading areas on supplemented by the list of dof in a region of interest. When is the identity matrix and F contains all the dof of the mesh, then (1) returns to the original finite element equations, .
Property 1. The following property holds: if the finite element solution is unique, if the solution of the hyperreduced equation is unique, and if the reduced basis is exactly reproducing the finite element solution, the hyperreduced solution is exact. Hence the following expression holds: .
In hybrid hyperreduced order models proposed in [24, 25], the reduced basis is extended with few finite element shape functions available inside . By construction, these shape functions are not connected to the remaining elements of . The set of indices of these shape functions is denoted by P. It is a subset of F, by construction. For the sake of simplicity, we order the degrees of freedom such that , where card(P) is the number of elements in the set P. Hence the reduced basis of the hybrid hyperreduced order model, denoted by , is the following block matrix: where is the identity matrix in dimension card(P). By substituting for in the hyperreduced balance equations, one obtains the hybrid hyperreduced equations, of which solution is denoted by . If the projection of on is exact, then it has also an exact projection on the larger subspace spanned by . According to Property 1, if is exact and if the hybrid hyperreduced equations have a unique solution, then the first components of should be zero and the last one should be equal to . It turns out that if is exact and if the hybrid hyperreduced equations have a unique solution, then the stress predictions and are equal. Hence, the following error estimator is proposed, in order to assess the accuracy of the hyperreduced order model that has been recognized by the computer vision workflow:where is the spatial domain occupied by the reduced mesh and c is a constant. The larger , the more complex the error estimation. When fulfills the finite element equilibrium equations, this error indicator is similar to the error indicator proposed in [26] for standard materials. If , if is the domain occupied by the full mesh , and if P contains all the dof indices, then is the true error. The constant c can be evaluated by following the procedure proposed in [25]. Here, we assume that .
As shown in Figure 1, is obtained by a noncentred PCA applied on simulation data . These data are displacement fields remapped on mesh from the mesh and restrained to the loads in the class β of the load clustering. This clustering of loading environments in the mechanical tests is presented below. A simple interpolation of the data in is done for the remapping of the displacement fields on the mesh for nodes in the domain occupied by . For nodes in that are not in the domain occupied by the mesh , the remapping is done via Laplace’s equation with an enforced continuity at the boundary of and outside . In practice, a robust model reduction is achieved if the meshes are restricted to elements that are not connected to nodes submitted to a concentrated load or a Dirichlet boundary condition. In the remapping procedure, the Dirichlet boundary conditions are enforced as Dirichlet boundary conditions of the Laplace’s equation. Then fulfills the Dirichlet boundary condition on the mesh .
The clustering of loading environments in the mechanical tests is performed by using the simulation data mentioned in the workflow in Figure 1. Unfortunately, these data do not have the same dimension since they are supported by different meshes . Then, they are remapped on the mesh of a bounding box that surrounds all the meshes , so comparisons are easier. The extrapolation of the data outside the meshes follows the Laplace’s equation again. The remapped simulation data on the bounding box are saved in the tensor of order 4. Four indices are introduced to have access to scalar values saved in . This value is denoted by , where i is the load case index, p is the dof index in the bounding box, j is the index of the mesh in the list , and n is the index related to additional parameter variations. For instance, local variations of the mechanical properties have the capability to enrich the simulation data for model-order reduction as proposed in [11]. When considering simulation data related to finite element models, a high resolution in spatial fields is achieved by high dimensional finite element space. Hence , the second dimension of , is often larger than 105 and can reasonably be up to 107 in industrial applications. In very high-dimension spaces, all the data are “far away” from the centre. Hence, a feature extraction is required in the framework of mechanical modeling, prior to clustering the simulation data. Several tensor decomposition methods are available in the literature for feature extraction. For instance, the k-PCA has been coupled to the proper generalized decomposition method in [27, 28] for extracting hidden model parameters. We refer the reader to [29] for a review on feature extraction. Here, we adopt a hierarchical Tucker format [30]:where is obtained by the following truncated singular value decomposition:where is the reshape of as a second order tensor: We restrict the clustering of data to features in dimension 2, for the ease of results visualization, and on the two following average data:This two-dimensional feature space enables visualizing clusters on loading environments. Here, we arbitrary select K clusters by using the k-means method. In the future we would experiment other clustering methods such as graph-based clustering, for instance. Then, we obtain a partition of the simulation data into K sets , …. Each set is a list of loading indices in . Then the original simulation data related to the cluster and the mesh are , for i in , p in and n in .
The CNN architecture chosen for this work is based on the layer composition initially described in [31]. The input to the CNN is a fixed-size RGB image of 246 × 246 pixels, denoted . This input image is obtained by downscaling the original 3968 × 2976 images denoted by in Figure 1.
The image goes through a set of convolutional layers followed by max-pooling layers. The convolution and pooling stack is repeated 3 times. The rectified linear unit (ReLu) activation function [8] is used for the convolutional layers. This function is denoted by . Each layer generates a feature map, denoted by after the layer. Each feature map is a tensor of order 3, which dimensions are denoted by . The first feature map is . The convolution layer applies linear filters to . Each linear filter is determined by the weights and the bias , such thatwhere , is the number of filters in the layer, , , and . The first convolution layer uses 32 filters; the next two use 64 filters. Here, the convolution kernel size is set to 3 × 3 and the stride to 1 pixel, no padding is used. The max-pooling is performed with windows of 2 × 2 pixels and a stride of 2. It readswhere , . The output of the last pooling layer is flattened and fed through 3 fully connected layers: the first one has 512 nodes and the second one 64 nodes; those first two layers use the ReLu activation function. The third layer is a soft-max layer with 4 nodes thus performing the 4-way classification task, such thatIn order to reduce overfitting, dropout regularization is implemented for the first fully connected layer with a dropout rate of 50%.
3. Results and Discussion
3.1. Experimental Setup
We consider a very simple mechanical test on a part in order to check its manufacturing process via the stress distribution in the part and the related response of the part submitted to various loading environments on the top of it. An image of the experimental setup is shown on Figure 2. The modeling via computer vision aims to recognize the loading environment applied on the part, in order to predict stresses by using a hyperreduced order model. The magnitude of the load and the vertical displacement at a fixed point are measured precisely during the mechanical test. But the location of the load has to be determined by the computer vision approach. We are considering L=18 possible loading cases regularly spaced on the top of the part.
The region of interest shown in red in Figure 3 is hidden in Figure 2 by the experimental setup. If this region of interest would be visible on each image of the experimental design, a digital image correlation [32] approach could have provided an estimate of the mechanical stresses in the part. But it is not appropriate here, because the part is partially masked in some digital images.
3.2. Simulation Data, Feature Extraction, and Clusters of Mechanical Loading
Prior to starting the manufacturing process, the experimental setup has been designed by using finite element simulations on four ideal geometries (M=4). The material of the part is elastic, for each mesh . The mechanical constitutive equations are the following:where is the symmetric part of the finite element shape function related to the degree of freedom of the mesh, E is Young’s modulus, and is Poisson’s coefficient. The mechanical properties of the material are E = 1,600 MPa, . Here, the residual of the finite element equilibrium equation reads Here, is the dof index where the load F is applied downward, is the spatial domain occupied by the mesh , and is the Kronecker delta. The observed magnitude of the mechanical load is 159.9 N.
For each geometry, three local variations of the Young modulus of -20% have been simulated by the finite element model, in order to enhance the space spanned by using the simulation data. The regions affected by these variations are shown in Figure 3. It turns out that we have done = 4 predictions of the displacement fields, for each geometry and each loading case. The region of interest is shown in red in Figure 3. It covers the subdomain where stress concentrations are expected to be highest during the mechanical test.
An example of the remapping of a vertical displacement field onto the bounding box is shown on Figure 4. The mesh of the bounding box has = 8. 105 dofs. The first feature mode (the first column of ) is shown in the right of Figure 4. This feature mode is consistent with the bending simulation of the mechanical parts.
We have arbitrary chosen K = 4 centroids to cluster the mechanical loadings. The L points are shown in Figure 5 with the 4 clusters presented by red circles.
3.3. Observational Data
About 250 high resolution digital images for each class of loading have been generated before starting any mechanical experiment. Examples of such images are shown on Figure 6. It takes approximately 10 minutes to get 1000 images. The size of each image is 3968 × 2976 pixels.
3.4. Training and Testing the Convolutional Neural Network
The CNN has been implemented with Keras library in TensorFlow. The train/test set was built following a 90/10 ratio upon the 1000 digital images of loading environments. The volume of training data was artificially increased by using a synthetic data augmentation strategy. Each “original” image was transformed using a combination of random sheer, zoom, and horizontal flip values. The sheer was limited to a maximum of 10% and the zoom to a maximum of 30%.
The training was performed by optimizing a multinomial cross entropy loss using a minibatch gradient descend approach with the RMSprop adaptive learning rate method; batch size is set to 32 and the model is trained on 120K steps (60 epochs).
The performance of the CNN is assessed on the test set; a top-1 error value of 1.9% was achieved. This great performance is explained by the easiness of the classification task due to no large variations between input images being observed, since they are all related to the same experimental setup.
3.5. Hyperreduction of Finite Element Equations
The modeling via computer vision has been applied to a mechanical test. The mesh of the part is in red in Figure 7. The smallest Hausdorff distance between the ideal meshes and is achieved by G4, so α = 4. This is the grey mesh superposed to in Figure 7. There is =8. 105 dof in .
The digital image of the mechanical part in the experimental setup is shown in Figure 8. The convolutional neural network recognizes the class number 3 of loading environment. Then β=3, and we can extract the simulation data from the database in order to compute the reduced basis , by using the noncentred PCA. We restrict this reduced basis to five empirical modes. Hence, the projection error of the selected simulation data on this reduced basis is less than 0.1%.
The reduced mesh is generated on the fly. It is shown in Figure 9. Here, F is a set of dof in the region of possible loading for class #3, plus the dof close to the displacement measurement, plus the dof in the region of interest. The elements of are the elements connected to F.
3.6. Stress Prediction and Error Estimation
A finite element simulation takes 45 min. The stress prediction by the hyperreduced order model shown in Figure 10 is obtained in less than 10 min. 99% of this computational time is spent for the solution of (1). By choosing a smaller set of degrees of freedom F, one can obtain prediction in less than 2 min. The larger the set F, the longer the simulations. As we can see in Figure 10, the elements of reduced mesh are not necessarily in a continuous domain. As explained in [11], the boundary conditions at the interface between and the remaining elements of are similar to Dirichlet boundary conditions. They are enforced by the empirical modes of the reduced basis. Then, the mechanical coupling between the discontinuous parts of the mesh is enforced by the empirical modes.
The exact error on the average stress in the region of interest is 0.1%. The map of the error estimator is shown in Figure 11. The dark grey elements in the reduced mesh (see Figure 9) are connected to the dof in set P. According to (1) by substituting the hybrid reduced basis for , the first finite element residuals are null, for i in . Hence the error estimator accounts partially for errors in finite element equilibrium equations. As shown in [26], these errors explain the discrepancy between the hyperreduced prediction and the finite element prediction of the stress.
4. Conclusions
The proposed reduced order modeling is related to very huge parameter space of dimension twelve millions, mainly due to the input image of the loading environment. The accuracy of the stress prediction is satisfactory to assess the quality of the process with mechanical considerations, even if the region of interest is hidden by the experimental setup. It is also reasonably fast in order to be inserted into a manufacturing process, aiming for part-specific decisions.
The output of the proposed workflow has a high spatial resolution. This is achieved by coupling a PCA, a clustering and a convolutional neural network. A local error estimator aims to indicate the discrepancy between the output and the stress that a finite element simulation would give corresponding to the loading environment recognized by the convolutional neural network. But this error indicator does not evaluate recognition errors, neither error on the mechanical behaviour of the observed material. So, the hyperreduced order model may not be the best that the available data could give.
In this paper, the inputs of the reduced order model are nonparametric loading conditions. They are defined solely by images of the loading environment. Since digital colour images are third order tensors, we expect a possible generalization of the workflow to more complex thermomechanical loading environments and more complex variations for the geometry of the observed parts.
Here, no Big-Data is required to train the proposed reduced order modeling via computer vision. The training starts with a nonsupervised machine learning by using a noncentred PCA and a clustering procedure, both on simulation data. Then, the CNN is trained on digital images by supervised machine learning upon the classes defined by the clustering procedure. Obviously, this approach can be implemented with larger sets of data.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
One of the authors is a Safran employee.
Acknowledgments
The experimental setup and the computer used for this publication have been funded by Mines ParisTech PSL Research University.