A Buckling Instability Prediction Model for the Reliable Design of Sheet Metal Panels Based on an Artificial Intelligent Self-Learning Algorithm

Abstract

Sheets’ buckling instability, also known as oil canning, is an issue that characterizes the resistance to denting in thin metal panels. The oil canning phenomenon is characterized by a depression in the metal sheet, caused by a local buckling, which is a critical design issue for aesthetic parts, such as automotive outer panels. Predicting the buckling instability during the design stage is not straightforward since the shape of the component might change several times before the part is sent to production and can actually be tested. To overcome this issue, this research presents a robust prediction model based on the convolutional neural network (CNN) to estimate the buckling instability of automotive sheet metal panels, based on the major, minor, and Gaussian surface curvatures. The training dataset for the CNN model was generated by implementing finite element analysis (FEA) of the outer panels of various commercial vehicles, for a total of twenty panels, and by considering different indentation locations on each panel. From the implemented simulation models the load-stroke curves were exported and utilized to determine the presence, or absence, of buckling instability and to determine its magnitude. Moreover, from the computer aided design (CAD) files of the relevant panels, the three considered curvatures on the tested indentation points were acquired as well. All the positions considered in the FEA analyses were backed up by industrial experiments on the relevant panels in their assembled position, allowing to validate their reliability. The combined correlation of curvatures and load-displacement curves allowed correlating the geometrical features that create the conditions for buckling instability to arise and was utilized to train the CNN algorithm, defined considering 13 convolution layers and 5 pooling layers. The trained CNN model was applied to another automotive frame, not used in the training process, and the prediction results were compared with experimental indentation tests. The overall accuracy of the CNN model was calculated to be 90.1%, representing the reliability of the proposed algorithm of predicting the severity of the buckling instability for automotive sheet metal panels.

1. Introduction

Denting resistance is an important issue in the design and manufacturing of sheet metal panels and structures and it is strongly influenced by the panel thickness and the manufacturing process conditions [1]. In addition to that, thanks to recent developments in computer-aided design and advanced manufacturing techniques, thinner and lighter complex surface topologies are being designed, even for commercial and mass-production vehicles [2]. Panel thickness, residual stresses induced by the manufacturing process as well as the surface curvature of panels [3] are all influencing factors for the denting resistance and, due to the strong interest of the automotive industry in this issue, several scholars focused their effort on studying this issue from the experimental as well the numerical points of view, as hereafter summarized.

From the experimental point of view, denting resistance can be assessed through indentation tests carried on sheet metal specimens, made of the material that is aimed to be utilized, or directly on the final product. Johnson and Schaffnit [4] investigated the denting resistance of cold-rolled low-carbon steels, used for the manufacturing of automotive outer panels, employing laboratory experiments where the dent depth at specific impact energy was measured. From the results of their analysis, they discovered that the denting resistance is proportional to the material yield strength and the square of the thickness, thus thinner plates require materials with higher yield strength to prevent denting.

In more recent work, Lu et al. [5] investigated the static and dynamic dent resistance of automotive body panels by means of laboratory experiments linking it to several designs and process variables, such as the material properties, the residual stress distribution after the forming process and the shape of the panel. One of the most interesting results of their analysis is that denting resistance is not only influenced by the material properties of the metal sheet but also by its surface curvature. As concerns the influence of the material properties and hardening behavior on the denting resistance, Shih and Horvath [6] employed dome test experiments and shown how a single denting loading, rather than an incremental ones, better replicate the denting conditions in real applications. In addition to that, in an interesting work by Holmberg and Thilderkvist [7], the influence of the original metal sheet properties and the stamping process parameters on the denting resistance were investigated on double-curved panels. From the results of this contribution, it can be inferred that the denting resistance is, as it can be imagined, strongly influenced by the original material’s properties of the metal sheet but it is also influenced by the blank holding force applied during the stamping process due to its influences on the residual stresses in the final product. Finally, Ekstrand and Asnafi [8] and Asnafi [9] defined a testing method to assess the denting resistance of automotive body panels and, in their analysis considered both major and minor curvatures of the panel as the main influencing parameters and that the boundary conditions, namely, the panel fixing jig’s geometry, have a strong influence on the results of the test.

However, the denting resistance of the original material and that of the manufactured panel are not equal due to the influences of the panel geometry and the manufacturing process. In addition to that, the application of strong fixing devices during the test, such as bolting or bead holders, may influence the results of the rest which may differ from the denting resistance of the final panel. For these reasons, in the research presented in this paper, denting experiments were carried out on an automotive panel already installed on a vehicle, and by considering different locations characterized by different geometrical features.

As concerns the development of numerical models for the prediction of the denting resistance of sheet metal panels, Holmberg and Nejabat [10] developed and validated a numerical model based on an automotive side-door, showing good accordance between numerical and experimental results in the low strain region, whereas the deviation increases for high strains. In [11], as also carried out in the research presented in this paper, the numerical model was developed neglecting the strain history in the panel caused by the manufacturing process. Similarly, concerning the finite element model development for the assessment of the denting resistance, Shen et al. [12] coupled different hardening models with the Baushinger effect to consider their influence on the denting resistance, showing its great influence. However, both the experiments and the numerical model were defined on test plates and not applied to real automotive panels, thus its reliability in real applications cannot be assessed. Finally, Park et al. [13] developed a numerical model for the estimation of the surface deflection of automotive body panels based on the surface curvature that allowed for the plotting of regions of the model where buckling instability is more likely to arise.

In the research presented in this paper, 20 numerical models relevant for 5 different panels belonging to 4 different vehicles were defined in ABAQUS 6.14/Standard (Dessault Systemes). The parts are defined as fender, front and back doors, hood, and roof and, for the sake of confidentiality, their relevant model is identified as M#1, M#2, M#3 and M#4, respectively. Numerical indentation analyses were carried out on different locations on these twenty parts, with 3 different thicknesses, resulting in a total of 1733 load stroke curves that were analyzed to investigate whether, in the considered load range of 0–20 kgf, buckling instability occurred. Afterward, through image analysis, the major, minor and Gaussian curvatures on the indentation points were extracted and correlated to the buckling instability [3,14]. Finally, the data set correlating the buckling instability to the considered curvatures was fed to a convolutional neural network (CNN) algorithm composed of 13 convolution layers and 5 pooling layers for training. Of the 1733 cases, 1386 (80%) were utilized for the training, whereas the remaining 347 (20%) were used for the first testing phase. Self-learning models have shown the capability of improving the performances of regression models [15] since their response can be continuously trained and updated, such as in the cause of machine learning [15,16]. For the case of pattern recognition, the application of deep learning in neural networks [17] allows correlating actions and effects with high accuracy, thus it was applied in this research. In recent years, artificial intelligence models were applied to various engineering topics, such as for the optimization of the friction welding parameters [18], the estimation of the radial-axial ring rolling process [19] and detection of cracking defects on stainless steels [20].

To validate the finite element model implementation procedure, experiments were carried out on some indentation positions of the above-mentioned panels and an additional part (door), completely neglected during the training and the first validation process of the CNN algorithm. By comparing the experimental and FEA load-stroke curves relevant for the same indentation points the average deviation was estimated to be 9.82%, defined as the ratio between experimental and simulated load-stroke curve area integrals. In addition to that, the proposed CNN-based model was utilized for the prediction of the buckling instability tendency on the above-mentioned additional door part, showing that the proposed algorithm can predict the higher or lower risk of instability buckling with an accuracy equal to 90.1%.

Thanks to the defined correlation between the denting resistance to the panel geometry, considered in terms of the major, minor and Gaussian curvatures, the proposed procedure can already be applied in the early stages of the design process and can be utilized by process engineers for the optimizations of the denting resistance with small changes in the aesthetical design. Moreover, the database can be extended by either adding more numerical, experimental or both, results, widening its application capabilities as well as its reliability.

Finally, to demonstrate the improvement in the prediction capability of the buckling instability phenomenon granted by the utilization of the proposed methodology, in comparison to previously published algorithms, the models proposed by Jung [1] and by Kim et al. [21] were applied to the experimental results cases carried out by the authors. From the results of this last cross-validation, it was possible to conclude that the developed methodology is also capable of predicting the buckling instability when its effect is not pronounced or, in other words, it has a higher prediction sensitivity.

In order to provide a clear insight into the main capabilities of the CNN-based model presented in this paper, the main findings of the research are highlighted as follows:

• The buckling instability of sheet metal panels can be estimated in the early design stages based on the modeling curvatures, avoiding expensive and time-consuming redesigns of the forming die;
• Both experimental and finite element results can be included in the training and validation data sets, allowing extension of the range of validity and applicability of the developed algorithm;
• The implemented image-based CNN methodology proved that machine learning algorithms can also already be utilized for optimization during the early stages of the design process. Moreover, although the methodology proposed in this paper was applied to sheet metal panels for the prediction of the buckling instability, it can be extended to different processes by accounting for the desired target function by applying the same implementation procedure presented in this paper.

2. Numerical Model Implementation

To train the CNN model, presented in Section 5, indentation finite element simulation models were implemented in ABAQUS/Standard considering the fender, front and back doors, hood and roof geometries. For each part, four different vehicle models, named M#1, M#2, M#3 and M#4, were considered for a total of 20 geometries. In Figure 1, the panels’ geometries, used for the training of the CNN model, are reported along with the M#5-Door, utilized only for validation purposes. For the sake of confidentiality, the dimensions of the parts, as well as specific details of the geometries, cannot be disclosed. Even so, since the proposed methodology is based on the local curvatures and not on the overall panel’s geometry, this latter information is not strictly related to the research. On the other hand, the local curvatures of the experimentally tested panels, also used for the finite element model validation, are summarized in the Appendix B of the paper.

All the panels are modeled considering the three initial thicknesses of the steel sheet, equal to 0.5, 0.6 and 0.7 mm, as normally employed in vehicle manufacturing. Moreover, as previously mentioned, since the proposed approach is aimed to be used in the early stages of the design process, where no details of the manufacturing process are available, the thickness variation related to the stamping process is not included in the analysis.

In the numerical simulations, all models were meshed with S4R elements, a 4-node general-purpose shell, reduced integration with hourglass control, for the smooth curvature portions of the model, and with S3 elements, a 3-node triangular general-purpose shell, for the remaining portions of the model. The need for utilizing two different types of elements is given by the complex geometries of the aesthetical details, where the utilization of S4R mesh may result in bad element aspect ratios. The average element side lengths are equal to 3.95 mm for the fender, 5.32 mm for the front door, 5.11 mm for the rear door, 8.07 mm for the hood and 10.0 mm for roof panels. The details of the mesh sensitivity analysis, for one of the considered roof panels, are described in Appendix A.

As concerns the boundary conditions applied to the model, the vehicle full-closed configuration was considered, thus fixed-boundary conditions were applied to the nodes belonging to the contours of doors, hood and roof, whereas the constraints introduced by the assembly bolts were considered for the fender, as shown in Figure 1 (red lines). During the indentation, both in the FEA model as well as in the experiments presented in Section 6, the indenter has a cylindrical shape with a radius equal to 12.5 mm and a fillet radius equal to 0.5 mm. In all the finite element models, the indenter was made to contact the panel surface with an orthogonal vector to the panel surface in the considered spot. Friction was modeled considering a Coulomb friction model with a friction coefficient equal to 0.2 in all the simulation models.

Since all the considered panels are manufactured with AISI-1008 galvanized steel, the same elastic-plastic material properties were utilized as reported in Equation (1), Swift-type hardening model, and shown in Figure 2, where experimental and numerical hardening curves are reported. As concerns the elastic mechanical properties, Young’s modulus of 206 GPa and Poisson’s ratio of 0.33 were considered. Elastic and plastic mechanical properties data were acquired from the results of tensile tests carried out on ASTM-E8 plate-type specimens carved from the undeformed AISI-1008 plate. Since the considered plates are manufactured through a rolling process, a slight anisotropy arises during the manufacturing operation, causing the rolling direction to have a higher yield and hardening than the transversal one. However, as previously mentioned, since the manufacturing process is not considered in this research, the anisotropy cannot be taken into account and, for this reason, the intermediate curve, relevant for the 45° direction, was considered for the estimation of the Swift model constants, Equation (1), as reported in Figure 2.

$${\sigma }_f=486.7{\left(0.008+\overline{\epsilon }\right)}^{0.202}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(\mathrm{MPa}\right)$$

(1)

3. Curvature Calculation Procedure

The algorithm proposed in the paper is aimed to be utilized during the early stages of the design process where no detailed information relevant for the manufacturing is normally available. For this reason, the authors decided to link the buckling instability to the major, minor and Gaussian curvatures, and to the initial sheet thickness; those can all be estimated from CAD files or, as in the case of this research, from the FEA model mesh files. The curvatures and thicknesses extracted from the panels’ mesh files were utilized as training input for the developed CNN algorithm, as presented in Section 5. The procedure presented in this section of the paper details the approach adopted for the extrapolation of the nodal curvature for the whole model based on its original mesh and is made necessary since the curvatures cannot be directly determined from the mesh file.

3.1. Quadrilateral Mesh into Triangular Mesh Conversion

The numerical models of panels are composed of the quadrilateral mesh (S4R) and the triangular mesh (S3) elements. To calculate the curvatures in the triangular mesh, the quadrilateral mesh is converted to the triangular mesh. The connected four nodes (1, 2, 3, 4) in the element are extracted from the mesh file and the 4-node elements are subdivided into two triangular meshes, (1, 2, 3) and (1, 3, 4) shown in Figure 3a. If the element has multi point constraint (MPC) conditions, the mesh is divided according to the number of MPCs, as shown in Figure 3b.

3.2. Vertex Normal Vector Calculation

After converting the quadrilateral mesh into the triangular mesh, the connectivity between the three nodes and the corresponding triangular element is checked and the list with each element and the corresponding three vertices is created. Afterward, for each one of the mesh elements (triangular) the element-wise normal vector ($N_i$), Figure 4a, is computed and is used for the calculation of the normal vector of the central vertex among the elements ($N_P$) in its surroundings according to Equation (2), as schematized in Figure 4b [22].

$$N_P\hspace{0.17em}=\hspace{0.17em}\frac{cN}{\left|cN\right|}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}where\hspace{0.17em}\hspace{0.17em}cN\hspace{0.17em}=\hspace{0.17em}{\displaystyle {\displaystyle \sum }_{i=0}^{n-1}\frac{N_i\sin {\alpha }_i}{\left|V_i\right|\left|V_{i+1}\right|}}={\displaystyle \sum _{i=0}^{n-1}\frac{V_i\times V_{i+1}}{{\left|V_i\right|}^2{\left|V_{i+1}\right|}^2}}$$

(2)

3.3. Nodal Curvature Calculation

According to the estimation of the normal vector for each vertex (node) of the triangular mesh, the relevant vertex curvature is estimated according to Equation (3) for the unit-length vector $\left(\begin{array}{cc}s& t\end{array}\right)$ in the plane with an orthogonal coordinate of $\left(\begin{array}{cc}u& v\end{array}\right)$ [23]. The unit-length vector is used to normalize the considered vector length to a $-1\le s\le 1$ and $-1\le t\le 1$.

$$\begin{array}{ll}{\kappa }_n\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}(s\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t)\prod \left(\underset{t}{\overset{s}{ }}\right)=(s^{\prime }\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{t}^{\prime })\left[\underset{0 \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\kappa }_2}{\overset{{\kappa }_1 \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}{ }}\right]\left(\underset{t^{\prime }}{\overset{s^{\prime }}{ }}\right)& \hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}{\kappa }_1{\left(s^{\prime }\right)}^2+{\kappa }_2{\left(t^{\prime }\right)}^2\\ \prod \hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}(D_{u}n\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{D}_{v}n)\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}\left(\underset{\frac{\partial n}{\partial u}v\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{\partial n}{\partial v}v}{\overset{\frac{\partial n}{\partial u}u\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{\partial n}{\partial v}u}{ }}\right)\end{array}$$

(3)

From the vertex normal vector information, as calculated in Section 3.2, the $\left(\begin{array}{cc}u& v\end{array}\right)$ coordinates system is constructed and the second fundamental tensor $\prod$ is computed.

The principal major (${\kappa }_1$) and minor (${\kappa }_2$) curvatures are estimated by calculating the eigenvalues of the second fundamental tensor whereas the Gaussian curvature (${\kappa }_G$) is determined by the multiplication of the major and minor curvatures $\left({\kappa }_1\cdot {\kappa }_2\right)$.

4. Buckling Instability Definition and Calculation

The buckling instability phenomenon is characterized by a change from positive to negative derivate in the load-displacement curve, as highlighted in Figure 5. In order to quantify the magnitude of the buckling, the integral of the area beneath the buckling region is divided by the area integral of the load-displacement reference curve, provided by the panels’ manufacturer, as shown in Equation (4).

The reference curve, as shown in Figure 5, defines the maximum displacement that the panel should show when a 20 kgf load is applied on its surface. In order words, the reference curve defined the minimum stiffness requirement for the panel and can be adjusted according to the requirements. Based on the ratio reported in Equation (4), the three classes of responses are defined as no canning ($AR=0$), soft canning ($0<AR<1$) and hard canning ($AR\ge 1$), respectively.

$$AR\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}{A}_{CANNING\hspace{0.17em}}\hspace{0.17em}\hspace{0.17em}/\hspace{0.17em}\hspace{0.17em}{\left({\displaystyle \int Lds}\right)}_{REFERENCE}$$

(4)

5. CNN Deep-Learning Algorithm Development and Training

This section of the paper is organized following the same flow utilized during the research for the development of the proposed neural network (NN)-based algorithm, thus, after an introduction of the theoretical backgrounds of the model, the focus will be placed on the methodology for the calculation of the considered three curvatures, utilized as input for the CNN algorithm. Afterward, the definition of the buckling instability, calculated from the results of the load-displacement curve, is provided along with the explanation of the pre-processing operations utilized for the construction of the training and validation data set for the CNN algorithm. Finally, the developed CNN model structure and validation procedure are presented, allowing to have an overall picture of the model background, implementation training and validation activities.

5.1. Neural Network (NN) Model Development

Neural network models consist of: (i) input layers, (ii) hidden layers and (iii) output layers, as can be seen in Figure 6a. Input layers are connected to the hidden layers by the weight functions (${\theta }_{i,j}$), which are calculated during the training of the NN algorithm. For the case of a 1D problem, for each one of the nodes, the inputs coming from the previous layer are defined as $x_i$, are multiplied by the weight functions (${\theta }_{i,j}$) and summed out to the bias values ($b_i$). The outputs of the layer are derived through the activation function ($f$), as shown in Equation (5) and Figure 6b.

$${\widehat{y}}_i=f\left({\displaystyle \sum {\theta }_{i,j}}\cdot x_i+b_i\right)$$

(5)

The logic for the calculation of the weights functions is based on Equation (6), where $E$ is the loss function and $\eta$ the learning rate. The standard NN model computes the weights to minimize the loss function, $E$, by applying the gradient descent and the backpropagation method [24,25]. However, to achieve a quicker convergence in the determination of the components of the filter matrix, ${\theta }_{i,j}$, the ADAM optimizer, presented in Equation (7), was utilized [26] in this research. The ADAM optimizer utilized a two-moment vector, ${\beta }_1$ and ${\beta }_2$, for the update of the filter matrix elements and it grants a higher convergence in comparison to the gradient descent function, which is based on a single moment vector.

Finally, in the developed algorithm, a categorical cross-entropy loss function is defined as in Equation (8), where $y_i$ is the true value and ${\widehat{y}}_i$ the predicted one, defined in Equation (5). The loss function is utilized to evaluate the quality of the prediction carried out by the trained NN algorithm during the training procedure.

$${\theta }_{i+1}={\theta }_i-\eta \frac{d{E}_i}{d{\theta }_i}$$

(6)

$$\begin{array}{l}{\theta }_{i+1}={\theta }_i-\eta \frac{{\widehat{m}}_{1,i+1}}{\sqrt{{\widehat{m}}_{2,i+1}}+\epsilon }\hspace{0.17em}\\ where\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{m}_{n,i+1}={\beta }_n{m}_{n,i}+(1-{\beta }_n)\frac{dE}{d{\theta }_i}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\&\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\widehat{m}}_{n,i+1}=m_{n,i}/(1-{\beta }_n^{timestep})\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}n\hspace{0.17em}=\hspace{0.17em}1,\hspace{0.17em}2\end{array}$$

(7)

$$E\left(y_i,\hspace{0.17em}{\widehat{y}}_i\right)=\hspace{0.17em}-{\displaystyle \sum _{i=1}^n{y}_i\log \left({\widehat{y}}_i\right)}$$

(8)

To avoid ‘vanishing gradient’ problems, arising if $dE/d{\theta }_i\le 0$, a ReLU (rectified linear unit) function $f$ was employed, Equation (9), and allows always obtaining a positive (>0) output from each one of the nodes of the hidden layer. If $x_i\le 0$, then the ‘answer’ of the relevant node is not considered. The reasoning behind this choice is given by the fact that, in the input layers, the above-mentioned three curvatures are always inputted as positive values and a negative or null number does not have any physical meaning and should be therefore neglected.

$$f(x)\hspace{0.17em}=\hspace{0.17em}\{\begin{array}{c}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}\hspace{0.17em}x\le 0\\ x\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}\hspace{0.17em}x>0\end{array}$$

(9)

For the output layer, the softmax activation function, defined in Equation (10), was applied to normalize the output values from 0 to 1, thus considering the highest number as the correct predicted answer. This strategy allows having a single answer from the trained NN algorithm which is subdivided into three classes: non-buckling, soft buckling and hard buckling.

$$f(x_i)\hspace{0.17em}=\hspace{0.17em}{e}^{x_i}\hspace{0.17em}\raisebox{1ex}{$$}\!\left/ \!\raisebox{-1ex}{$$}\right.\hspace{0.17em}{\displaystyle \sum _{i=1}^3{e}^{x_i}}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}i=1,2,3$$

(10)

In this research, the input layers for the NN algorithm are defined as the 2D image arrays of local major, minor and Gaussian curvatures, whereas the three output layers are set to predict non-buckling, soft buckling and hard buckling. The detailed explanation of the physical interpretation for the non-buckling, soft buckling and hard buckling outputs were provided in Section 4 of the paper.

The CNN methodology was utilized to convert the curvature distribution images into the input layers for the NN algorithm. The CNN algorithm is utilized for image recognition [27] and allows transforming an image into an array, whose features are then extracted by applying a convolution filter. The employed convolution filter has a 3 × 3 size and, as shown in the example reported in Figure 7, the filter was defined to extract only the principal diagonal values of the input image array. The three diagonal values, defined as ${\theta }_{i,j}\cdot x_{i,j}$ when $i=j$, are fed into Equation (5), where the subscript $j$ is added to the input value $x$ to account for the 2D nature of the inputted image, not considered in the 1D representation of Figure 6. Accordingly, Equation (5) becomes Equation (11), where the applied convolutional filter is shown as well.

$${\widehat{y}}_{i+k,j+k}=f\left({\displaystyle \sum _{i=1,j=1}^3{\theta }_{i,j}\cdot x_{i,j}+b_i}\right)$$

(11)

According to the model defined in Equation (11), the weights, the cross-entropy function of Equations (7) and (8), and the activation function of Equations (9) and (10), a CNN model composed of an input layer made of 3 channels, 22 hidden layers, and 3 channels for the output layer were implemented in Python 3.7.4 and Keras 2.3.1 with TensorFlow backend framework.

The details relevant for the pre-processing of the curvature images extracted from the mesh file, as presented in Section 3, are reported in Section 5.2, whereas the structure, training and cross-validation of the developed CNN algorithm are reported in Section 5.3, respectively.

5.2. Pre-Processing

The curvature distribution images of the major, minor and Gaussian curvatures are regarded as grayscale images in the graphical GUI where white is considered as 255 whereas black as 0. This conversion allows inputting the curvature distribution as a number included in a predefined range. The contour intervals are considered separately along with each curvature distribution, where the major curvature has 45 discrete sections of 0.09 to −0.009, the minor curvature is sorted out in 45 sections of 0.009 to −0.09 and the Gaussian curvatures has 36 intervals between 0.009 and −0.009.

These ranges were chosen according to the results of the three considered curvatures, all of which fall inside the above-mentioned ranges. In all images, a black background is added outside the original picture contour to avoid false readings during the conversion from image to channel input values. By considering the boundary projected area (top view) of the considered panel, the curvature distributions are converted into grayscale BMP images and a pixel-to-mm conversion ratio is calculated. This conversion ratio is not constant but depends on the size of the panel in the x–y plane (Figure 8).

In order to correlate the above-mentioned curvatures with the buckling instability data, a series of indentation points were chosen on the panel surface and the surrounding 50 mm × 50 mm indentation area, identified according to the indenter dimensions. Since the number of pixels in the indentation area is not constant, it varies from panel to panel, the original number of pixels is upgraded to 256 pixels by 256 pixels image. By following this procedure, all the images have the same number of pixels along x- and y-directions and can be used for a uniform and reasonable training of the NN algorithm. On the selected indentation points, finite element indentation analyses, as summarized in Section 2, were carried out and, according to the procedure detailed in Section 3 and Section 4, allowed linking the local curvatures distributions with the three classes of responses, namely, no canning, soft canning and hard canning.

By following this procedure, the training set, made of 1733 local curvatures distributions and relevant class assignment (569 no canning, 412 soft canning and 752 hard canning), along with the relevant thickness, was constructed and utilized for the training of the proposed CNN algorithm, as detailed in Section 5.3 of the paper. Furthermore, for the case of the validation panel, the same procedure was utilized but the results were not utilized for the training of the CNN algorithm.

5.3. CNN Model Structure, Training, and Cross-Validation

To predict the buckling instability phenomenon arising in thin metal panels according to the surface curvatures, as summarized so far, a combination of a pre-trained weight matrix in the VGG-16 model [28] and a self-developed classification layer was utilized according to the architecture shown in Figure 9. The architecture consists of a features extraction part and a classification part.

In the feature extraction part, the model consists of 13 convolutional layers and 4 max-pooling layers. The self-developed classification part has a global average pooling layer and dense layers, where the last dense layer has 3 nodes that are designed to provide a non-buckling, soft buckling and hard buckling classification. The drop-out method in the classification section is utilized after a dense layer to overcome possible overfitting problems.

The layers highlighted with the symbol ‘*’ in Figure 9 possess pre-trained weight matrices calculated with the ImageNet dataset [27] and are not updated during the model training. These non-trainable weight matrices are based on a huge classification of datasets and help the trained model to recognize the images by pre-trained weights matrix in the early stages, improving the classification results even if the images are different from that of the ImageNet [27,29,30]. In the convolutional layers, a progressively increasing number of filters is applied to recognize ever-smaller features in the images. Hence, by following this procedure, the number following the word ‘Conv’ in Figure 9 represents the number of filters in the relevant layer, the higher the number of extracted features the lower the size of the relevant feature. This results in the progressive zoom-in of the images, as schematized in Figure 10, where the results are relevant for the output of the first convolution of each convolutional block, from Conv,64 to Conv,512, for the case of the filter matrix reported in the detail of Figure 10. In Figure 10, the extracted features from each of the images are highlighted in yellow and green.

After the convolution layers, the featured images are resized in half by the max-pooling layer to amplify the features and reduce the size of the computed data. After the feature extraction, the 2D image arrays are flattened to an array of 1 row by the global average pooling (GAP in Figure 9) layer allowing matching with the averaged values of a 2D image into a single value for one node.

After exiting from the Conv,512 layer, and after the GAP application, the features are inputted in the Dense, 512 layer and are reduced to 256 features in the dropout layer with a 0.5 dropout ratio [31]. Because the dropout layer selects the 256 nodes randomly at each step during the training process, the model does not depend on specific nodes, allowing overcoming of the overfitting problem. In this step, the selection of the 256 features to be inputted into the Dense, 256 layer is carried out and the randomly selected 128 nodes are linked to 3 nodes in the last layer, thus allowing obtaining of the prediction of no canning, soft canning or hard canning. The hyperparameters are tuned as 0.00005 for the learning rate, with 10⁻⁸ for $\epsilon$, 0.9 for ${\beta }_1$ and 0.999 for ${\beta }_2$ in Equation (7), 128 batch size and 300 learning steps according to computing power.

The learning procedure is carried out by using the k-fold cross-validation process. In order to evaluate the prediction model for the non-trained data, the total dataset is randomly subdivided into k sets. One set is utilized as the test set for estimation of the training results whereas the remaining four sets are used for the training of the CNN algorithm. In this study, the 5-fold cross-validation process was utilized and the best model among the 5 split models was selected to predict the buckling instability of the target panel geometry.

The final results for the test set are evaluated by the average value of classification accuracy and $F_1\hspace{0.17em}$ score of the five split model results. The $F_1\hspace{0.17em}$ score is calculated by following Equation (12) as a harmonic mean value of precision and sensitivity, as in Equations (13) and (14). In Equations (13) and (14), $\mathrm{TP}$, $\mathrm{FP}$ and $\mathrm{FN}$ are defined as true positive, false positive and false-negative prediction, respectively.

$$F_1\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}2\times \frac{Precision\hspace{0.17em}\times \hspace{0.17em}Sensitivity}{Precision\hspace{0.17em}+\hspace{0.17em}Sensitivity}$$

(12)

$$Precision\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}\frac{\mathrm{TP}}{\mathrm{TP}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}\mathrm{FP}}$$

(13)

$$Sensitivity\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}\frac{\mathrm{TP}}{\mathrm{TP}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}\mathrm{FN}}$$

(14)

6. Denting Resistance Validation Experiments

To validate the finite element model presented in Section 2, indentation experiments were carried out on 5 different panels belonging to the five categories presented in Figure 1. Moreover, 22 indentation tests were carried out on the M#5-Door, only used for the validation of the developed NN algorithm. The six tested parts, with the relevant considered indentation points, are reported in previous Figure 1.

For the measurement of the load-displacement curves during the indentation experiments, a 6-axis robot, Figure 11, having the indenter as the final end-point, was employed. The dimensions of the indenter cylinder are reported in Figure 11 and are the same utilized in the FEA simulations, described in Section 2 and Appendix A.

The robot is equipped with a load cell, allowing measuring of the load applied during the test, whereas the displacement of the end-point is automatically measured by the robot arm encoder. The indentation points and their relevant coordinates on the panel surface are summarized in

Table A1. Numbering and coordinates of the validation point experiments for 0.7 mm thickness (dimensions in mm).

Model	P#	x	y	z	Maj.Curv.	Min.Curv.	Gauss.Curv.	RES
M#3 Fender	#1	650	100	30	−1.80 × 10−4	1.23 × 10−3	−2.23 × 10−7	N
	#2	550	0	68	1.76 × 10−4	1.66 × 10−3	2.91 × 10−7	N
	#3	650	−50	70	1.37 × 10−4	1.20 × 10−3	1.65 × 10−7	N
M#1 Front Door	#1	370	−265	135	−1.19 × 10−3	−1.37 × 10−5	1.64 × 10−8	N
	#2	485	−97	105	−2.56 × 10−5	1.65 × 10−3	−4.21 × 10−8	S
	#3	795	−95	105	−3.59 × 10−5	7.26 × 10−4	−2.60 × 10−8	H
M#2 Roof	#1	500	0	112	−9.08 × 10−3	−2.90 × 10−4	2.63 × 10−6	N
	#2	465	320	108	−9.65 × 10−3	−5.03 × 10−5	4.85 × 10−7	S
	#3	600	320	119	−8.60 × 10−3	−5.38 × 10−5	4.62 × 10−7	H
M#4 Rear Door	#1	240	−200	105	−1.81 × 10−3	−1.62 × 10−5	2.92 × 10−8	N
	#2	860	−230	100	−9.16 × 10−4	2.56 × 10−5	−2.35 × 10−8	N
	#3	450	−450	103	−5.69 × 10−4	1.86 × 10−4	−1.06 × 10−7	N
M#3 Hood	#1	450	750	245	−2.49 × 10−4	−8.62 × 10−5	2.15 × 10−8	N
	#2	−200	700	251	−1.90 × 10−4	−1.14 × 10−4	2.17 × 10−8	H
	#3	0	450	210	−2.29 × 10−4	−1.96 × 10−4	4.50 × 10−8	H
M#5 Door	#1	130	−130	88	−1.78 × 10−3	−3.64 × 10−5	6.46 × 10−8	N
	#2	350	−80	74	−4.43 × 10−2	4.67 × 10−5	−2.07 × 10−6	N
	#3	550	−120	80	1.61 × 10−5	8.40 × 10−3	1.35 × 10−7	N
	#4	220	−190	100	−3.87 × 10−3	−2.48 × 10−5	9.60 × 10−8	H
	#5	350	−200	100	−3.05 × 10−3	4.51 × 10−5	−1.38 × 10−7	N
	#6	470	−210	100	−1.51 × 10−3	1.20 × 10−4	−1.80 × 10−7	N
	#7	580	−240	104	−1.61 × 10−3	1.56 × 10−4	−2.51 × 10−7	N
	#8	680	−200	99	−1.81 × 10−3	1.05 × 10−4	−1.90 × 10−7	N
	#9	130	−330	104	−9.75 × 10−4	2.23 × 10−5	−2.17 × 10−8	N
	#10	290	−330	103	−9.41 × 10−4	1.20 × 10−4	−1.13 × 10−7	N
	#11	450	−330	105	−7.78 × 10−4	2.56 × 10−4	−2.00 × 10−7	N
	#12	590	−330	110	−6.26 × 10−4	2.87 × 10−4	−1.80 × 10−7	N
	#13	210	−420	100	−6.80 × 10−4	1.44 × 10−4	−9.82 × 10−8	N
	#14	370	−420	100	−5.59 × 10−4	4.75 × 10−4	−2.65 × 10−7	N
	#15	530	−420	107	−4.29 × 10−4	5.25 × 10−4	−2.25 × 10−7	N
	#16	130	−510	88	1.30 × 10−4	1.65 × 10−3	2.13 × 10−7	N
	#17	290	−510	88	−3.96 × 10−4	6.17 × 10−4	−2.44 × 10−7	N
	#18	450	−510	98	−6.66 × 10−4	8.57 × 10−4	−5.71 × 10−7	N
	#19	460	−35	55	−3.16 × 10−3	5.37 × 10−5	−1.70 × 10−7	N
	#20	595	−10	42	−2.26 × 10−5	8.93 × 10−3	−2.02 × 10−7	N
	#21	750	−20	25	−3.57 × 10−4	9.15 × 10−3	−3.26 × 10−6	N
	#22	875	−20	54	−6.11 × 10−3	−2.66 × 10−4	1.63 × 10−6	N

(Appendix B of the paper) whereas the load-displacement results are reported in the result section where they are compared with the results relevant for the implemented FEA models. In Table A1, in the ‘RES’ column, the letters ‘N’, ‘S’ and ‘H’ stand for no canning, soft canning and hard canning, respectively.

7. Results

To prove the reliability of the implemented FEA model, the results of the indentation experiments, in terms of load-displacement curves, were compared with those of the relevant numerical simulation model. The results of the 18 load-displacement curves are reported in Figure 12 whereas the experimental and numerical area integrals, along with the percentage errors, are reported in Figure 13 and detailed in

Table A2. Experimental vs. numerical load-displacement curve’s area integrals and errors.

Model	P#	EXP. Area (mm2)	FEA Area (mm2)	Error (%)	AR	Model	P#	EXP. Area (mm2)	FEA Area (mm2)	Error (%)	AR
M#3 Fender	#1	44.45	47.38	6.59	0	M#5 Door	#5	158.65	172.12	8.49	0
	#2	21.26	22.83	7.40	0		#6	131.62	143.57	9.08	0
	#3	50.01	54.72	9.43	0		#7	84.09	93.45	11.13	0
M#1 Front Door	#1	49.95	54.65	9.40	0		#8	80.75	72.60	10.09	0
	#2	79.94	89.57	12.05	0.07		#9	96.79	89.22	7.82	0
	#3	72.49	81.27	12.11	2.59		#10	77.05	86.17	11.84	0
M#2 Roof	#1	313.01	354.07	13.12	0		#11	101.89	115.26	13.12	0
	#2	229.83	248.36	8.06	0.54		#12	121.12	107.23	11.47	0
	#3	204.45	219.63	7.43	2.52		#13	76.76	68.94	10.18	0
M#4 Rear Door	#1	94.43	103.66	9.77	0		#14	222.01	192.16	13.44	0
	#2	100.60	90.23	10.31	0		#15	326.91	355.21	8.66	0
	#3	97.51	86.13	11.67	0		#16	99.85	109.74	9.91	0
M#3 Hood	#1	132.75	145.46	9.58	0		#17	265.83	244.01	8.21	0
	#2	244.89	264.16	7.87	15.95		#18	376.46	344.23	8.56	0
	#3	340.96	379.95	11.44	2.89		#19	328.76	360.36	9.61	0
M#5 Door	#1	72.19	79.31	9.87	0		#20	370.89	402.18	8.44	0
	#2	177.16	159.32	10.07	0		#21	268.97	301.27	12.01	0
	#3	135.57	124.34	8.29	0		#22	312.28	285.94	8.43	0
	#4	144.81	132.39	8.58	3.09

(Appendix B of the paper).

As concerns the developed buckling instability prediction model based on the CNN architecture, its reliability was verified in two different steps. In the former one, the k-fold cross-validation process was employed and the relevant accuracy and $F_1\hspace{0.17em}$ score were calculated and the results are summarized in

Table 1. Results of test-set in 5-fold cross validation.

n-Fold	1	2	3	4	5	Average
Accuracy	0.908	0.886	0.899	0.914	0.897	0.901
$F_1\hspace{0.17em}$ score	0.909	0.898	0.900	0.921	0.907	0.907

. The calculations were carried out considering a total of 347 points relevant for the panels shown in Figure 1. As previously mentioned, the M#5-Door points, Figure 1, were not used for any of the model training sets since they are aimed to be used only for validation purposes.

For the calculation of both accuracy and $F_1\hspace{0.17em}$ score, 80% of the whole data set was used for the CNN algorithm training, whereas the remaining 20% for its validation. All five sets reported in Table 1 are different from each other and were randomly selected to avoid any bias in the validation process.

The latter validation was carried out considering the prediction for buckling instability on the 22 points of the M#5-Door (Figure 1), where, from the results of the experiments, only point #4 showed hard canning, whereas the remaining 21 showed no canning. Moreover, to verify the reliability of the buckling instability prediction in the considered thickness range of 0.5–0.7 mm, a 0.55 mm thickness M#1-Hood part was numerically tested as well. As concerns the 0.55 mm thickness M#1-Hood, since no experiments are available, the comparison was carried out considering the reference curve criteria for the definition of the relevant canning class, as shown in Section 4. The detailed comparison is reported in

Table A3. Coordinates and results of the validation points for 0.55 mm thickness M#1-Hood (dimensions in mm).

Model	P#	x	y	z	Maj. Curv.	Min. Curv.	Gauss. Curv.	Pred. Res.	FEA Res.
M#1 Hood	#1	750	450	177	−9.90 × 10−3	−2.16 × 10−5	2.14 × 10−7	S	S
	#2	700	−200	187	−2.17 × 10−4	−1.09 × 10−4	2.36 × 10−8	H	H
	#3	450	0	466	−2.60 × 10−4	−2.54 × 10−4	6.62 × 10−8	H	H
	#4	250	−600	43	−5.60 × 10−4	3.79 × 10−3	−2.12 × 10−6	N	N
	#5	550	650	121	−2.38 × 10−3	−3.16 × 10−4	7.53 × 10−7	H	H
	#6	700	−700	137	−1.00 × 10−2	−2.63 × 10−4	2.64 × 10−6	H	H
	#7	850	750	137	−3.76 × 10−4	4.88 × 10−3	−1.84 × 10−6	S	N
	#8	250	−450	68	−1.03 × 10−3	−3.32 × 10−4	3.43 × 10−7	N	N
	#9	450	550	112	−1.09 × 10−3	−3.53 × 10−4	3.84 × 10−7	H	H
	#10	450	−450	114	−7.45 × 10−4	−2.72 × 10−4	2.02 × 10−7	S	S
	#11	650	500	146	−6.96 × 10−4	−2.25 × 10−4	1.56 × 10−7	N	N
	#12	900	−550	175	−1.02 × 10−4	5.83 × 10−4	−5.96 × 10−8	H	H
	#13	250	350	73	−4.68 × 10−4	8.77 × 10−3	−4.10 × 10−6	N	S
	#14	600	−400	148	−1.43 × 10−2	−2.29 × 10−4	3.27 × 10−6	N	N
	#15	200	−200	82	−6.59 × 10−4	−3.66 × 10−4	2.41 × 10−7	H	H
	#16	500	200	147	−3.60 × 10−4	−1.95 × 10−4	7.02 × 10−8	N	N
	#17	250	0	102	−5.83 × 10−4	−3.80 × 10−4	2.22 × 10−7	H	H
	#18	850	0	466	−1.52 × 10−4	−1.35 × 10−4	2.05 × 10−8	N	N

(Appendix B of the paper).

Before the prediction of the buckling instability of the M#5-Door, the whole 60 panels’ data set’s (20 different panels and 3 thicknesses) were fed into the CNN algorithm and 3 output channels were identified as red for hard canning, green for soft canning and blue for no canning. The results of the buckling instability for the 22 points of the M#5-Door (Figure 1) and the 18 points for the 0.55 mm thickness M#1-Hood are reported in Figure 14a,b, respectively. The load-displacement curve for point #4 of the M#5-Door, red points in Figure 14a is that of Figure 12.

Finally, to demonstrate the improvement in the prediction of the buckling instability phenomenon in comparison to previously published models, the two contributions of Jung [1] and Kim et al. [21] were considered. The two models detailed in these two contributions were applied to the panels’ geometries presented in Figure 1 and, along with the predictions carried out by the developed model, were compared with the experimental results summarized in previous Figure 12. The comparison among these three prediction models, and the relevant experimental results, are reported in

Table 2. Buckling instability prediction results comparison among developed and literature models for the prediction of the buckling instability of thin metal sheets.

Panel Model	P#	Experiments	Proposed Model	Jung [1]	Kim et al. [30]	Panel Model	P#	Experiments	Proposed Model	Jung [1]	Kim et al. [30]
M#3 Fender	#1	N	N	N	N	M#5 Door	#5	N	N	N	N
	#2	N	N	N	N		#6	N	N	N	N
	#3	N	N	N	N		#7	N	N	N	N
M#1 Front Door	#1	N	N	N	N		#8	N	N	N	N
	#2	S	S	N	N		#9	N	N	N	N
	#3	H	H	N	N		#10	N	N	N	N
M#2 Roof	#1	N	N	N	N		#11	N	N	N	N
	#2	S	S	N	N		#12	N	N	N	N
	#3	H	H	N	N		#13	N	N	N	N
M#4 Rear Door	#1	N	N	N	N		#14	N	N	N	N
	#2	N	N	N	N		#15	N	N	N	N
	#3	N	N	N	N		#16	N	N	N	N
M#3 Hood	#1	N	N	N	N		#17	N	N	N	N
	#2	H	H	H	H		#18	N	N	N	N
	#3	H	H	N	H		#19	N	N	N	N
M#5 Door	#1	N	N	N	N		#20	N	N	N	N
	#2	N	N	N	N		#21	N	N	N	N
	#3	N	N	N	N		#22	N	N	N	N
	#4	H	H	H	H

where N, S, and H stand for no canning, soft canning and hard canning, respectively.

The overall accuracy of the two considered literature models is good but they fail to recognize the soft-canning case, where the canning area, defined in Section 4, Figure 5, is not pronounced. This is clear for the case of the M#1-Front Door point #2 and M#2-Roof point #2, where the two literature models predict no canning whereas both experiments and developed CNN model predict experienced and predicted soft canning. Moreover, even if the considered canning area leads to a hard canning evaluation if the result is close to the soft canning, prediction inaccuracies have also been observed, such as the case of the M#1-Front Door point #3, M#2-Roof point #3, M#3-Hood point#3.

Through this additional validation, it is possible to conclude that the proposed methodology, based on the CNN method, is not only capable of accurately predicting the buckling instability phenomenon but it is also more sensitive to the transition between soft and hard canning, further enhancing its reliability.

8. Discussion

Based on the comparison between the experimental and FEA results concerning the load-stroke curves reported in Figure 12, some differences are present between the two curves. These differences are mainly caused by the two simplifications introduced in the FEA model, as hereafter summarized. On one hand, for the case of the real panels, the forming process causes a non-uniform thickness distribution, which also slightly alters the local denting resistance. On the other hand, the forming process also causes the arising and relief of residual stresses, which influence the local mechanical properties, thus also the denting resistance. However, since the proposed approach is aimed to be utilized during the early stages of the design process, only the conceptual model might be available, thus the consideration of the forming process, normally defined in a later stage, might not be feasible. For this reason, in the proposed approach, the thickness was modeled as constant considering three different levels, a representative for the sheet thicknesses mostly used in the automotive industry. In addition to that, in order to avoid any possible bias, isotropic mechanical properties were considered as well.

Although these two simplifications introduce some errors into the calculation, the comparison between experimental and FEA load-stroke curves results in the ratio between the integral areas equal to 9.82% being reliable. For the sake of conciseness, not all the results relevant for the tested panels can be included in the paper but the above-mentioned average deviation between experimental and numerical load-stroke curves was calculated on the whole result set and not only on the results included in the paper.

As concerns the conversion strategy from mesh file to grayscale images, as shown in Figure 15, the higher profiles, protruding from the panel surface, are considered to have negative major and minor curvatures, thus a positive Gaussian curvature. For this reason, in Figure 15 for the case of point #2 (minor curvature), the light grey at the bottom of the picture represents a higher z-coordinate (orthogonal to the panel surface) in comparison to the top of the picture (black). In the region in the middle between the light grey and the black, the change in curvature is visible as a progressive grey shading.

Similarly for point #2, the consideration of the Gaussian curvature allows avoiding loss of information when the images are supplied to the CNN algorithm, overcoming the difficulty in identifying surface features in cases of high-darkness images. By employing this approach, the training and validation sets, as presented in Section 5.3, were constructed and the high accuracy and $F_1\hspace{0.17em}$ scores allow concluding that the considered dataset construction procedure is balanced and not biased. This important aspect should be considered if the proposed algorithm is aimed to be employed to different components or in case the implemented dataset is aimed to be expanded by considering additional FEA or experimental results.

9. Conclusions

In the research presented in this paper, a methodology for the estimation of the buckling instability in sheet metal panels based on the CNN theory was defined and validated on a wide data set and shown to be accurate and reliable even though some simplifications were introduced in the development of the FEA model utilized for the construction of the training and validation data sets. From a global perspective, the developed image conversion approach and CNN algorithm for the prediction of the buckling instability of panels was shown to be a reliable method which, thanks to its reasonable simplifications, can also be utilized in the early stages of the design process, allowing to correct the panels’ geometries if and where required. The overall accuracy in the estimation of the buckling instability, calculated to be 90.1%, was estimated by comparing the predictions carried out by the developed model with the observations relevant for 1733 cases, proving the strong correlation between major, minor and Gaussian curvatures with the panel’s denting resistance. In addition to that, the proposed model can also be extended to different panels’ geometries and utilized in a more advanced manner by considering a combination between sheet metal forming and structural simulations, allowing to account for the local thickness distribution and residual stresses.