Calculating the gravity-free shape of sheet metal parts

This paper presents an iterative finite element (FE)–based method to calculate the gravity-free shape of nonrigid parts from an optical measurement performed on a non-over-constrained fixture. Measuring these kinds of parts in a stress-free state is almost impossible because deflections caused by their weight occur. To solve this problem, a simulation model of the measurement is created using available methods of reverse engineering. Then, an iterative algorithm calculates the gravity-free shape. The approach does not require a CAD model of the measured part, implying the whole part can be fully scanned. The application of this method mainly addresses thin, unstable sheet metal parts, like those commonly used in the automotive or aerospace industry. To show the performance of the proposed method, validations with simulation and experimental data are presented. The shown results meet the predefined quality goal to predict shapes within a tolerance of ± 0.05 mm measured in surface normal direction.

considering their real shape (compare W. W. Cai et al. [1]). With these simulation models, often referred to as digital twin, optimizations can be performed and more precise predictions can be made. Especially, for the automotive and aerospace industry, this principle is used for optimizing the assembly of thin parts, for example, sheet metal body parts as presented by A. Rezaei Aderiani et al. [2]. Before simulations with digital twins were performed in the quality assurance processes, only physical measurements were compared. The measurements had to fulfill special requirements. For example, the position and the clamping of the individual part during measurement need to be similar to the assembled situation. This is necessary to be able to compare the result with a measurement of the assembled situation. To achieve such a condition, the part needs to be over-constrained to prevent deflections caused by its weight. By that, the part becomes tensioned and its geometry is influenced. However, to be able to run correct simulations with the measurements, the gravityfree shape is required. To obtain this shape, a non-overconstrained measurement is necessary. This problem arises, as measurement fixtures are expensive, and performing two measurements on two different fixtures is inefficient and expensive. Nevertheless, over-constraint measurement is necessary for quality assurance.
This paper solves this problem by using an iterative FE based approach to calculate the gravity-free shape based on a non-over-constrained measurement. Therefore, an FE simulation model is set up from the measured geometry. The simulation result enables the opportunity to perform quality measures directly on the gravity-free shape without using expensive fixtures and eliminating the error caused by the over-constraining at the same time. As the gravity-free shape is represented as a simulation mesh, the part can be virtually put in any other situation by performing additional simulations. For example, virtual clamping can be applied to obtain the shape of the part, clamped into any measurement fixture. By applying this concept, only one measurement on one fixture has to be performed and the part can be virtually put in any desired situation to perform quality measures on it.
The paper is structured as follows: first, definitions established in Section 2, followed by a literature review about the topics affecting this paper, are given in Section 3. Next, a detailed description of the problem is given in Section 4. After that, the method of calculating the gravityfree shape is explained and an overview of the steps that need to be performed will be given in Section 5. Then, the method is applied to virtual and physical experiments in Section 6. In Section 6.1, the method is applied to simple test geometry and also an industrial use case and in Section 6.2, to a physical test setup to validate the virtual experiments. Last, the results are discussed in Section 6.4, and the paper is concluded in Section 7. In particular, the contributions of the following work are: -Algorithm to calculate the gravity-free shape of a simulation model that is suitable for almost all commercial FE solvers. -Validations with simulation data for a real-world use case. -Experimental validation with a simple, easy to understand sheet metal part.

Definitions
This section discusses the theoretical foundations affecting the work of this paper. In particular, phrases that are defined differently in the literature are clarified.

Free shape
The free shape is the geometry of the part situated in the free state. Unfortunately, there are different definitions of the phrase "free state." In the literature sometimes, the term is used to describe the situation of a weightless, load-and constraint-free state. This paper sticks to the definition from DIN EN ISO 10579 [3]: Free state (translated) State of a part, where nothing but gravity load caused by its own weight is applied.

Gravity-free state
3D scans with a non-over-constrained fixture obtain one possible shape of the part in the free state. As the method of this paper is aimed to obtain the shape without the influence of any load including gravity, the term "gravityfree" (hereinafter "GF") is used to describe this load-free state. In Fig. 1, a simple 3D example shows the GF shape and the free shape with a non-over-constrained fixture.

Nonrigid parts
The general sense of the term "nonrigid part" describes an object that deforms under its weight in a relevant magnitude. As every object deforms due to gravity, a classification is needed to be able to differentiate between rigid and nonrigid parts. G. N. Abenhaim et al. [4] introduced such classification by using the term "compliance behavior." This basically describes a ratio between the applied force and the occurring deflection of the part. Based on this measure, three classes of parts are introduced. DIN EN ISO 10579 [3] uses another definition that is more application dependent: Nonrigid part (translated) A part that deforms in the free state with an amount that exceeds the dimensional and/or geometrical tolerances defined in the technical drawing For this paper, both descriptions are applicable but, for consistency, the paper sticks to the definition from DIN EN ISO 10579.

Background
In the literature, two approaches for finding the GF shape of a measured part are present. The first approach is to use a simulation mesh of the measured state and set up an inverse FE formulation to iteratively calculate the gravity-free state. The second approach uses a FE simulation of a CAD part to calculate the deformations, caused by its weight when it is held by the measurement fixture. Then, the displacements of the simulation result are mapped to the measurement in the opposite direction. For both methods, the present work will be highlighted and discussed in the following. Besides that, the literature for fixtureless inspection methods is highlighted as the presented method also has common goals.

Inverse finite element formulation/inverse form finding
To obtain the stress-free shape from an FE model representing the loaded shape, an inverse FE simulation can be used. Besides the FE model of the loaded shape, the load case that led to the deformed shape must be known. With this information, an inverse FE problem can be formulated to calculate the unloaded shape. Inverse formulations can be used for a wide variety of mechanical problems. S. Germain and P. Steinmann [5] made a comparison between an inverse mechanical formulation and Limited-Broyden-Flechter-Goldfarb-Shanno method for shape optimization, showing that an inverse formulation is much faster in terms of computing time. Later, the same authors S. Germain et al. [6] proposed a recursive algorithm for solving inverse form finding problems in isotropic elastoplasticity. P. Landkammer et al. [7] again used the same algorithm for solving an inverse form-finding problem for orthotropic plasticity. A. Ask et al. [8] set up and validated an inverse formulation for the electro-elastic case. All presented publications are showing promising results; however, for using an inverse formulation in real applications, an implementation for a commercial simulation software is necessary. Besides that, all publications on inverse form finding have in common that no experimental validations are showing the usability for real use cases.

Mapping of gravity deformation
To bypass the problem of generating costly simulation models from measurement data, an intuitive way is to take advantage of the similarity between the measured part and its CAD model. Therefore, a simulation model from the CAD geometry is loaded with the constraints of the measurement fixture and gravity load. The displacements of the simulation result are then mapped to the measured point cloud to update the vertex positions of the measurement. C. Lartigue et al. [9] used a voxel filter for the point registration on the point cloud to map the nodes of the simulation mesh to the measurement. A commonly used and old method for this kind of registration problem is multidimensional scaling which is used in a variety of medical and mechanical applications; see [10][11][12]. A detailed description is given by I. Borg and P. Groenen [13] and a survey on the method is proposed by N. Saeed et al. [14]. H. Radvar-Esfahlan et al. [15] improved this registration method by calculating and comparing geodesic distances on the part surface to refine the mapping. F. Thiébaut et al. [16] presented a validation with experimental data showing the use of mapping simulation data for evaluating shape deviations of nonrigid parts. Although the problem of mapping simulation results in a measurement representing a similar geometry seems to be solved satisfactorily, this method is always based on the assumption that the geometries, CAD, and measured part behave in the same way.

Fixtureless inspection methods
To save the cost for expensive measurement fixtures, different approaches are dealing with the post-processing of optical measurements that do not need complex fixtures or no fixtures at all. The most common fixtureless approach is called virtual clamping. The basic idea of virtual clamping is to generate an FE model of the scanned part and apply the boundary conditions of any desired measurement fixture. The simulation results deliver a result that can be compared to the CAD model. V. Tuominen [17], for example, is using the principle of virtual clamping for the inspection of subframe components from an automotive application. More recently, V. Sabri et al. [18] present an algorithm where the CAD mesh is iteratively deformed toward the measurement by a continuous displacement field. By checking neighborhood conditions of the displacement field, areas, where profile defects like dents are occurring, can be identified. This work was extended by S. Karganroudi et al. [19]. A survey of fixtureless inspection methods is given by G. N. Abenhaim et al. [4]. The challenge of post-processing point clouds that were obtained from a measurement without fixture is to distinguish between deflection caused by its weight and geometrical differences compared to CAD. Although the presented methods can differentiate between profile defects and geometrical deviations caused by loads or variabilities, the methods are not capable of comparing the part with the CAD in a free state. Similar to the principle of "virtual clamping" is to perform a virtual assembly, based on measured geometries. In the work of I. Gentilini and K. Shimada [20], a full pipeline from generating the FE model of the measurement toward predicting the postassembly shape of the measured part is proposed. Here, also different materials and objects, for example, plastic lid of a storage box, or a metal sheet housing of a toaster, were investigated. Nevertheless, the authors are stating that to be able to use this method properly with weak parts, the influences of gravity and clamping must be removed in a previous step.

Research questions
The presented related work is showing some solutions or partial solutions for finding the GF shape of a measurement. All have in common that they either need a special implementation for the used solver or are based on similarity and need an implementation for solving the point registration problem. Also, the literature does not provide any validations with experimental data for the case where the CAD geometry cannot reflect the real behavior of the measured gravity-free state. To apply any kind of gravity compensation method to industrial use cases, both problems, the implementation to commercial software and the validation with experimental data for the general case instead of special cases, must be addressed.

Problem description
Before describing the method, the problem of measuring nonrigid parts is explained more in detail. Usually, measurements represent a compromise between deflections due to weight and ambiguity of the measurement. Figure 1 shows a simple 2D example for a measurement. The part is constrained at the right end and displaced due to gravity with maximal displacement influence toward S 0 and the displacement (vector) field u 0 illustrated via a few representative vectors. The free shape is the geometry that is measured and is the input of our proposed method. The goal is to determine the unknown displacement field u 0 and obtain the unknown GF shape from the measured part; see Eq. (1). The goal for the approximation is to achieve a deviation measured in surface normal direction of less than 0.05 mm. The value of 0.05 mm is the measurement uncertainty of typically used optical measurement systems for the quality assurance of sheet metal parts; see W. Boesemann et al. [21].

Deflections of the part on a measurement fixture
When a part is measured, it is usually placed in a measurement fixture. This fixture ensures that the part is always positioned in the same way to be able to reproduce measurements. Various optimization methods exist to determine an optimal amount and positions of locators for measuring or machining a part; see [22][23][24]. Especially, for measuring nonrigid parts, the goal of an optimization is to minimize Fig. 1 Problem definition. Example of a simple 2D beam constrained at the right end; S 0 represents the free shape resulting from the influence of gravity (u 0 ) deflections (caused by weight) by using as few locators as possible. Figure 2 shows an example demonstrating the problem. Situation (a) depicts the GF of the desired part geometry (CAD geometry) positioned on two rigid support elements. As the part is considered as a nonrigid, it deflects when applying gravity to it. This case is depicted in (b). To prevent large deflections during the measurement process, an additional support element is added in the center in (c). By adding the additional support element, the measurement obtained from (c) is comparable with the CAD geometry model.

Unambiguity of measurement
Adding additional support elements to the fixture to reduce deflections during the measurement process has one significant disadvantage. The measurement yields an ambiguous geometry. Figure 2(a) shows an actual part geometry, not undergoing the influence of gravity, that differs from the CAD geometry. After applying gravity (b) and adding the additional support element (c), the shape can no longer be differentiated from the measurement of the CAD geometry. The missing information is the distribution of the forces that act on the support elements. In Fig. 2, (a) and (b), the distribution is unambiguous since the forces on the support elements can be calculated merely by the center of mass relative to the contact points of the support elements. For the case depicted in (c), at least one force must be given to calculate the forces on the other two support elements. Unfortunately, the forces acting on the support elements are usually not measured during a 3D scanning process. Even if the forces were known, the measurement still would have to be post-processed to make the difference visible for an inspection engineer.

Fig. 2
Unambiguous boundaries used for comparison of ideal and actual geometry for different constraints/loads. a GF shape for the two geometries. b Ideal and actual geometry loaded with gravity influence. c Additional support used (blue) to compensate the part's deflection; the geometries now agree qualitatively As over-constrained measurements are often used in industrial applications, the relevance of this problem is demonstrated via a real-world example, shown in Fig. 3.
Here, a simulation model of an engine hood was used to show a worst-case scenario. The hood has dimensions 1290 mm × 1490 mm × 415 mm (X, Y, Z) and is made of seven individual steel sheet metal components. All components were meshed with first-order shell elements (number of nodes 365,700), joined together by modeling all spot welds and adhesives. The input for this demonstration is a fictitious badly shaped part, depicted on the left-hand side in the gravity-free state. This part has deviations up to 25 mm when compared to the desired CAD geometry. By applying the boundary conditions of a potential fixture using four support elements in Z-direction (perpendicular to the image plane) and load the structure with gravity, the shape shown in the middle is obtained. This shape is viewed as the measured geometry. The shown fixture setup is very similar to real fixture layouts used in automotive applications. An intuitive approach for gravity compensation would be this one: for example, one maps an inverted displacement field, obtained from an FE simulation calculating the deformation of the CAD geometry, on the same fixture. This step is shown in the right picture in the figure. As can be seen, the resulting shape is not even close to the known GFshape. The reason for this large difference of up to 24.6 mm was motivated in Fig. 2. The distribution of the weight to the support elements is usually not obtained during the measurement process, and it is assumed to be equally distributed in the simulation. In this particular case, the weight is distributed as unevenly as possible due to the bad input shape. The image shown in the middle of Fig. 3 depicts the forces on the support elements, indicated with two plus signs and two zeros. The two support elements with the plus signs carry together 100% of the weight. The other two support elements are merely touched and are used only to keep the part in balance. Thus, the real deformation due to the weight of the part is not recognized and is hidden by the over-constrained fixture. Most methods described in the literature point out this problem and only consider measurements that are not over-constrained.

Similarity between CAD and real geometries
FE simulations based on ideal CAD geometries are commonly used to predict deformations caused by the weight of the part when it is placed in a measurement fixture. Simulation results can be mapped to 3D scans to correct a measured shape as discussed by F. Thiébaut et al. [16]. These methods typically assume that a measured part geometry behaves similarly to the simulation model that predicts the deformation caused by the weight of the part. This assumption must be treated with caution. Already minor differences in shape can cause recognizable differences in simulation results. Figure 4 illustrates this problem. Here, a steel sheet metal part with dimensions 300 mm × 200 mm × 0.5 mm is simulated using a linear elastic material model. Two different geometries were used. The top row shows an ideal CAD geometry represented by a plane, and the bottom row shows a part with small shape deviations below 0.5 mm, compared to the CAD geometry. Both geometries were meshed with first-order shell elements (number of mesh nodes 15,251) and loaded equally with gravity-induced forces applied perpendicularly to the image plane (Z-direction). Boundaries are shown in the figure. To demonstrate the influence of the orientation of the part relative to the fixture, two simulations were performed-one with gravity in positive and one with gravity in negative Z-direction. One can also interpret this as two different parts with inverse shape deviations. For the ideal plane, both simulations lead, as expected, to the same displacement field. For the part with shape deviations, two different displacement fields resulted from the FE simulation, and neither one of them matches the result of the simulation obtained for the CAD geometry. The two displacement fields obtained for the part with shape deviations differ up to 0.2 mm, which is evident by direct comparison (right-hand side).
These results show that minor shape deviations can cause recognizable errors when assuming that real part geometries behave like ideal geometries. Also, different fixture layouts and shape deviations were tested, resulting in a sensitive behavior of the outputs. Thus, it depends greatly on shape deviation and fixture layout whether it is possible to use an FE simulation result, produced for an ideal CAD geometry, as a basis for performing shape corrections on a measurement by mapping simulated displacements. In real-world applications, fixture layout is known but shape deviations are unknowns. Therefore, it is difficult to verify the quality of a mapping approach as it can be different for individual parts.
To overcome the described shortcomings, the authors present an iterative FE simulation-based approach that calculates a valid GF shape of a measured sheet metal part by using measured geometry as simulation input. It is assumed that the input measurement is not over-constrained to prevent effects like the ones depicted in Fig. 3.

Method
This section is split into two parts. First, the method is explained by discussing mathematical aspects of calculating the GF shape. Second, the used algorithm resulting from the mathematical description is presented. To calculate the GF shape, the method needs the following inputs: (1) a mesh representation of the measured object, (2) the locator positions and directions, (3) material information, and (4) the gravity direction. The output of the computation will be a mesh representation of the GF-state of the measured object.

Generation of the simulation model
The first step after measuring a part is the generation of a simulation model from the measured geometry. There are many different methods available for the generation of a mesh from a measured geometry. I. Gentilini and K. Shimada [20] for example are applying a bubble packing method to the point cloud to reduce its density. After triangulating the coarsened point cloud, they define the elements as FE shell elements to run simulations. Another method that automatically creates volume FE meshes from STL files is proposed by Y. Liu et al. [25]. More recently, also a simulation method called finite cell method became more popular. This method enables running simulations completely without the use of meshing the measured object. It can handle three-dimensional point clouds from different kinds of scans. A detailed review and a description of the method are given by D. Schillinger and M. Ruess [26]. Another possible way to obtain an FE mesh from a measurement is to use commercially available software tools for reverse engineering to get a CAD representation of the measured object and then again meshing the new CAD model. Thus, the generation of the simulation mesh is treated as a problem on its own that can already be solved in a satisfying way for different kinds of applications.
To finalize the simulation model material properties, boundaries and loads must be added. The boundaries and loads can be derived from the measurement fixture. As mentioned, the fixture must be unambiguous to calculate the GF shape only with geometrical information as input (compare Section 4.2). For further explanations, the load case shown in Fig. 1 is used.

Iterative approach (Loop 1)
Starting from the initial situation, the most intuitive way of finding the GF shape would be to invert the direction of the load. This might work for really small magnitudes of u 0 and in some special cases but is wrong in general. As such, it is not applicable reasoned by the difference of geometrical stiffnesses. The pre-studies in Section 4.2 already discussed this problem. Motivated by the problem of different magnitudes by flipping the gravity direction, an iterative approach is used to approximate the displacement field u 0 . The basic idea of this approach is to subdivide the total gravity load into smaller increments Δg and applying the load incremental wise to the structure. Per iteration, two simulation runs are performed, one with the load of +Δg, and in a separate simulation, the same incremental load is applied in the inverse direction. Figure 5 shows the resulting displacement fields u for every node. The resulting shape S 1 is the updated mesh configuration after the first iteration, being the initial configuration for the next iteration. The shape S 1 can be expressed mathematically as: The displacements u ±Δg 0 can be written with: (3) in (2) results in: To approximate the target displacement field u 0 , all incremental displacements need to be summed up: In every iteration, only node positions of the mesh are updated. Tensions and forces are not considered for the next iteration. When updating the position of a node, the geometry of the simulation model is changed. Thus, the displacement field calculated in the next iteration changes as well. This iterative approach solves the problem of different geometrical stiffnesses in different directions of gravity. To obtain the GF shape approximation, iterations are performed until the sum of all Δg values equals the value of g, captured mathematically in Eq. 5. Note that one could think, performing only the simulation run with +Δg and flipping the displacement vectors might be sufficient, as the direction of displacement is very similar to the simulation run with −Δg, for small Δgs. Unfortunately, it turned out that this slight difference in direction is summed up while performing the iterative steps leading to recognizable errors. The user controls increments for this "Loop 1." By defining how many iterations are to be performed, the value of Δg is given as Δg = g/#iterations. Similarly to other iterative approaches, it is expected that the quality of the approximation improves with an increasing number of iterations, at the expense of computational cost. The ideal number of iterations is a trade-off between precision and computational cost and must be determined for every application individually via a sensitivity study. Once the ideal number of iterations is determined, this number can be used for the calculation of the GF shape for further measurements of the same part.
Iterative approaches can be used to approximate a solution. By reducing the increment size resulting in a higher number of iterations the solution should become more precise. It turned out that the presented approach has a weak spot that causes a convergence toward a shape that slightly differs from the desired GF shape.

Stress stiffening (Loop 2)
To find the reason for the error in the approach of Section 5.2, first, the solution of one iteration is verified. To do so, two simulations are compared. The first simulation is applying g + Δg to the known GF shape to obtain the configuration S +Δg 0 from Fig. 6. The second simulation is applying the load +Δg on the free shape S 0 . Note that the shape S 0 was calculated before by applying the load +g to the GF shape and exporting the mesh. The initial meshes of both simulations are considered stress-/strain-free. It turned out that these two simulations do not result in the same geometry. As a reason, an effect called stress stiffening is identified, which is not considered in the second simulation which causes the differences in the simulation results.
Stress stiffening, sometimes also named geometric stiffening, incremental stiffening, initial stress stiffening, or differential stiffening, is the effect that a loaded structure behaves differently from a stress-free structure. This effect usually occurs in simulation with thin structures like shells or beams. Depending on the geometry, the structure can be stiffened or weakened by this effect. The mathematical aspects are not discussed in detail; they can be looked up for example in the Ansys Theory Guide [27]. For further improvements, it is only important to know that the stress stiffening effect influences the behavior of a loaded structure.
As the stresses are lost after the mesh update in each iteration of Loop 1, the stress stiffening effect cannot be considered directly. To overcome this problem, a second iterative approach is used, based on the calculations made before as described in Section 5.2. As mentioned before, a shape that does not exactly match the desired GF shape can be calculated, but this shape is already very close to the GF shape. The shape obtained from Section 5.2 is further used as initial geometry for the following calculations. The basic idea for the next iterative simulations is to use the similarity of the GF shape and the already calculated shape. Both shapes result in similar displacement fields when applying the full gravity load. As the stress stiffening effect is considered in the forward simulation, the resulting displacement field will differ from the approximation calculated before in Section 5.2. This scenario is depicted in Fig. 7a, where S Loop1 is the shape initially used from the previous iterative calculations and S 0 the free shape. Now, the displacement field u g that has been obtained from loading S Loop1 with the full gravity can be used to displace the free shape S 0 in the opposite direction. The resulting shape S Loop1 is closer to the desired GF shape (see Fig. 7b). This behavior can be expressed as: As both geometries, S 0 and S Loop1 , are represented by the same mesh, the displacements can be applied per node. So a registration method for matching the geometries is not necessary. By running the forward simulation from Fig. 7 and updating the shape S Loop1 multiple times, the resulting geometry comes every step closer to the GF shape. So the following sequence can be established: The limit of S i is the desired GF shape: As the change of the shape from S i to S i+1 can be calculated, a termination criterion can be defined as: The maximal change of the shape converges to zero for an infinite number of iterations. The threshold can be viewed as a measure for the precision of an intermediate result.
This approach allows a user to control output precision for a specific application by specifying the value of . Note that the whole calculation of the GF shape might work by only applying "Loop 2" multiple times directly to the measured shape S 0 , but this would result in a large number of iterations that can be drastically reduced by using the geometry of Loop 1 as input for the second loop.

Algorithm
This section will describe an algorithmic implementation of the method described before. This algorithm highly relates to the method proposed by S. Germain et al. [6]. The difference in the method which is proposed in this paper Fig. 7 Two-step optimization of "Loop 2." a First, full gravity load g is applied to the shape S Loop1 ; b second, a displacement field u g is applied in inverse direction to the shape S 0 , resulting in S Loop1 is that it does not use any inverse FE formulation. This might result in way higher run times (see Section 6.3) but makes the method easier to implement with almost any simulation software and thus easier and capable to use in real applications.
To explain the algorithm in detail, Fig. 8 is used. As a prerequisite, a simulation model of the measured geometry, containing boundaries, material information, and direction of gravity, is necessary. The algorithm is split up into two loops. These loops refer to the iterations explained before in Sections 5.2 and 5.3. Before starting the computation, the user has to define the number of iterations for the first loop and the abort criteria for the second loop. At the beginning of Loop 1, a Δg is defined by dividing the full gravity load by the defined number of iterations. Then, two simulations can be started in parallel, one with the load vector +Δg and one with −Δg. From these two results, the magnitude and direction of the displacement for the node update can be obtained. This step refers to the principle explained in Section 5.2. After the node positions of the mesh nodes are updated, the next iteration of "Loop 1" can be performed with the updated mesh. When all iterations are done, the mesh geometry already is close to the GF shape. The updated mesh geometry is then used as input for Loop 2. This loop refers to the iterations explained in Section 5.3. Here, the following steps are performed. First, the full load of +g is applied to the resulting mesh from Loop 1. Then, the direction of the displacement vectors becomes flipped and applied to the original input mesh that represents the measurement of the free state. The resulting mesh is closer to the GF shape then the result from Loop 1 and is used for the next iteration of Loop 2. This loop is running until the maximum change of the displacement field from one iteration to the next one is below the user-defined value . The result of the algorithm is a simulation model of the part in its GF shape. From here, the user can run any forward simulation to virtually clamp or assemble the part in any desired situation to perform quality measures in other situations if needed.
As stated in the "Abstract," this algorithm does not need a CAD model of the part and also does not use any inverse formulations. The user merely needs access to the Fig. 8 Flowchart of algorithm to calculate the GF shape. The algorithm is divided into two loops: columns represent steps that are performed during one iteration. Loop 1, shown in blue, iterates a fixed number defined by the user. Loop 2, shown in red, terminates based on a threshold condition node position of the mesh and the displacement vectors of the results to implement the proposed algorithm. Commonly, these are fundamental functions provided by scripting interfaces of commercial simulation software. An alternative visualization of the method is given in Algorithm 1 expressed as pseudocode.

Validation
This section provides validations to show the performance of the proposed algorithm. There are two kinds of validations that are presented in this section. The first kind are virtual validations showing the performance of the algorithm at simple and complex parts. These validations have a ground truth (CAD) to measure the precision of the method. The second kind of validations are validations with experimental data to show the whole workflow and precision based on a physical experiment.

Validation with simulation data
To show whether the proposed method works in theory, virtual experiments were performed. In these simulations, a CAD model was deformed by a gravity load with unambiguous boundary conditions. The resulting shape then was exported. With this mesh representing the free state, a simulation model preserving the same material model and boundaries was set up. The mesh is considered as tensionand stress-free. This simulation model is the input for the algorithm (see Section 5.4). The result of the algorithm can then be compared to the CAD model which represents the GF shape.

Simulation setup
To show the performance of the proposed method, the algorithm from Fig. 8 is implemented as a plugin in the commercial simulation software Abaqus from Dassault Systèmes. The kernel of this implementation only takes about 175 lines of code. With this prototype plugin, the GF shapes of the two simulation models shown in Fig. 9 are calculated. On the upper row, a simple piece of sheet metal with the dimensions 200 mm × 300 mm × 0.5 mm is shown. As material, aluminum with a simple linear elastic material model is chosen. On the lower row, a more application relevant geometry is depicted. This geometry represents the roof of a car with the dimensions of about 2000 mm × 1000 mm, with 0.7 mm thickness. The material is chosen as steel, also with a linear elastic material model. Both geometries have meshed with the default Abaqus shell element, a first-order shell element. To avoid uncertainties caused by a coarse mesh resolution, both geometries were meshed densely (using 15,251 nodes for the simple sheet metal and 169,378 nodes for the car roof). Besides the geometry, the supports of a fictitious measurement fixture are shown in the first column "Initial setup." The direction of the gravity load is chosen perpendicular to the image plane (Z-direction). The next column "Applied gravity" shows the influence of gravity compared to the CAD geometry from the first column. This result is obtained by running a simulation that applies the gravity load to the meshed CAD geometry with the boundaries already shown. The resulting mesh is then considered stress-free and is used as an input simulation model representing the free state and retaining the same material model and boundaries.
The user-defined parameters for the algorithm were 4 for the number of iterations (n) of Loop 1 and = 0.004 for the abort criteria of Loop 2 (compare Section 5.4). These parameter values were obtained by performing simulations using the CAD geometry. There exists an optimal number of iterations for Loop 1 for every part. The value of the parameter n was varied between 1 and 10. In this particular case, we found that n = 4 was optimal for the number of iterations of Loop 1 for both geometries. When the user does not set the number of iterations properly, the second loop will take more iterations but will nonetheless find the correct GF shape. The second parameter = 0.004 is related to the precision of the output and can be chosen as desired.

Results
The last two columns of Fig. 9 show the simulation results of the proposed method. The results of both loops are shown separately in columns 3 and 4 of the figure compared with the CAD geometry. The differences between the geometries resulting from Loop 1 and their corresponding CAD geometry are below 0.5 mm which is a good first approximation but for quality assurance purposes way too coarse. After the second loop of the algorithm, the deviations between the CAD geometry, which represents the GF shape, and the corresponding simulation results were reduced drastically. Both results differ below a value of 0.019 mm which is less than half than the defined quality goal from Section 4.

Validation with experimental data
To show that the proposed method not only works in theory but in real applications as well, a physical experiment is demonstrated in this section. Therefore, the physical setup is explained first in Section 6.2.1. Next, the simulation setup is explained in Section 6.2.2. Lastly, the results are presented in Section 6.2.3. The goal of this validation with experimental data is to prove that the GF shape of a sheet metal part can be calculated from an optical measurement. This cannot be validated directly without sending the part together with its fixture and a measuring machine to space. So, this validation works indirectly by performing the steps shown in Fig. 10. First, the specimen is scanned in situation 1. From this scan, a simulation model is created by fitting a surface to the measured point cloud. Then, the fitted surface becomes meshed and boundaries, as well as loads from situation 1, are applied to the model. The next step is to calculate the GF shape by using the proposed algorithm. Then, another simulation step is performed, where the boundary conditions of situation 2 are applied to the model of the GF shape. The last step is to compare the simulation result with the measurement of situation 2. Only if the GF shape is calculated properly, the comparison succeeds.

Experimental setup
The setup of the physical experiment is depicted in Fig. 11. The specimen is a steel sheet metal part (structural steel St37 -1.0037) with dimensions of 200 mm × 300 mm and a thickness of 0.7 mm. As support, four balls that are in-plane are used for restricting the part in Z-direction (direction of gravity). For positioning the specimen in the X and Y directions, precision rods are used as boundaries.
The specimen and the rods are sprayed with a mat metal primer to reduce reflections that cause faulty measurement results when using optical measurement systems. As an optical measurement system, an HP Structured Light Scanner Pro S3 was used. This device provides scans with a precision of 0.05 mm. Note that the authors have used different experimental configurations (i.e., lighting conditions, calibrations, angle of view) for acquiring the measurement data. The highest quality point cloud is used for the validation presented in this section.
The specimen was scanned in the two situations shown in Fig. 11. Situation 1 has unambiguous boundaries. This is the measurement from which the initial simulation model is created. Situation 2 has an additional support in Z-direction which results in ambiguous boundaries. This is the situation that is simulated by calculating the GF shape of situation 1 and applying the boundaries of situation 2 in a second simulation step. Note that it is only possible to calculate the shape of situation 2 based on the shape of situation 1. Contrariwise, it is not possible because the measurement of situation 2 is not unambiguous. To ensure that no plastic deformation occurs when placing the specimen on the fixture, an FE simulation with an ideal geometry (plane) was performed to estimate the mechanical stresses for situation 1. Maximal stress is obtained at about 5% of the yield strength for the steel used.

Simulation setup
The first step of the workflow right after measuring the part is to create a simulation model of the measured geometry.
To do so, a surface fitting approach is used, where a parametrized surface rectangle with the correct dimensions is fitted to the point cloud. Therefore, the commercial software tool Alias from Autodesk is used. The fitted surface  With the simulation model of the measured geometry created, the simulations can be performed to obtain the GF shape.

Results
After performing the steps described in Fig. 10, the results are shown in Fig. 12. The first row shows the measurements of the same object in two situations (see Fig. 11) compared to a plane surface with the dimensions 200 mm × 300 mm which is considered as CAD for the following comparisons. All measurements are aligned by using the contact points of the part with its fixture as alignment boundaries for a 3-2-1 alignment (see Fig. 11). The measurement of situation 1 shows that the part behaves similarly to the part from the virtual validations shown in Fig. 9. The left picture of the second row shows the result of the surface fitting compared to the measurement of situation 1. It can be seen, that the measurement is noisy on a scale of ±0.05 mm but due to the surface fitting, the resulting surface is smooth. Based on this fitted surface, the simulation mesh is created. The right picture of the second row shows the result of the proposed algorithm which calculated the GF shape compared with the CAD model. Note that in this comparison the color transition is much smoother because the simulation mesh and the CAD can both be assumed as smooth. The last row shows the result of the simulation with the boundary conditions of situation 2 applied to the calculated GF shape, compared with the actual measurement of situation 2. Note that the measurement of situation 2 was also smoothed by fitting a surface to it as well. Ideally, the shown comparison would not show any deviations, but as actual measurement data is used for this comparison there will always be differences. The picture shows that the differences between measurement and simulation result are below ±0.05 mm measured in the surface normal direction. The results and run times of both the virtual and physical validations are summed up in Table 1.

Run times
To compare the proposed algorithm in terms of run time with other methods, the simulation time of the iterations were determined and summed up. The shown numbers (Table 1) are only covering simulation run time, as the pre-and post-processing times are highly dependent on the implementation. All simulations were performed on the same system: -CPU: AMD Ryzen 5 1600 @ 3.20 GHz (8 threads were used) -16 GB DDR4 RAM -Windows 10 64 bit LTSB -Abaqus Vers. 2018 Fig. 12 Result stages obtained from validation with measured data. The largest node error value for the final comparison, considering the entire mesh used to generate this image, is less than 0.05 mm. This value is less than the tolerance value published by the manufacturer of the 3D scanner used to digitize the object

Discussion
This section will take a deeper look into the results of the validations with simulation data presented in Section 6.1. By using a CAD model, a precise ground truth model is given to compare with. The only differences that can be measured are either a consequence of the numerical uncertainty of the FE solver or result from the proposed algorithm. The message file of the simulations contains displacement corrections of the order of 4e−05, which is far below the scale used to evaluate the quality of the GF shape. Therefore, the results, shown in Fig. 9, are trustworthy. The results show that, as long as the simulation model does reflect the measured geometry (unambiguous measurement) and behavior of a real part, the proposed method can be used for calculating the GF shape. The deviation of less than ±0.019 mm is sufficient for the considered application and meets the design goals for our method as defined in Section 4. In both results, it can be seen that the areas that are far apart from the fixture supports have the highest error which is plausible. Nevertheless, one can improve the results further by performing additional iterations of Loop 2; see Fig. 8.
The measured run times shown in Table 1 only consider simulation times. The time required for solving the problem of applying full or partial gravity load to a structure is the bottleneck for reducing the overall time needed for a solution. The pre-and post-processing steps can be performed fully parallel, in principle. To reduce the overall run time of the proposed method would require one to reduce either the number of iterations or the complexity of the simulation. In terms of run time, the proposed method was not designed to compete with inverse FE formulations as described by S. Germain et al. [6], for example, but it still provides a viable alternative approach for industrial applications where real-time support is not necessary. The advantage of the proposed method is the simplicity of the implementation of the algorithm compared to implementations using inverse FE calculations. Our method's overall run time cost is also relatively high when compared to approaches that only map a precalculated gravity deformation field to a measurement. This cost is caused mainly by the expensive generation/adjustment of a simulation model for each scan. By performing all calculations based on measured geometry, the assumption of similar behavior between the ideal and real geometry no longer applies. Furthermore, there is no need to register the measured point cloud to the CAD geometry. The GF shape, calculated with our method, is valid for all part distortions. This increased reliability and precision justifies the higher run time cost, especially in applications where high precision is crucial.
Our validation results show that our method is suitable for real applications. The chosen specimen geometry shows a displacement-per-length ratio that is above the car roof from the virtual validations; see Section 6.1. The small test setup reflects deformations similar to those obtained when measuring large sheet metal parts. The results of the validations are shown in Fig. 12. The deflections of the part caused by gravity are plausible when comparing the results with the used boundary conditions. The second row on the left side in the figure shows a comparison between the surface fitting and measurement (the fitting target). The deviations produce a stripe pattern when compared to the fitted surface. Structured 3D light scans and the resolution of the measurement system explain this effect. The last row of the figure shows the quality of the validation with experimental data. The measured data had been smoothed by fitting a surface to the original point-cloud. Deviations are below ±0.05 mm, meeting the design goal of < ±0.05 mm. This is sufficient for assessing the geometric quality of sheet metal parts; see Section 4. This statement holds for the whole surface being compared and not merely for points of interest. The experimental results presented show lower maximal error values than experiment results presented in the literature, e.g., F. Thiébaut et al. [16]. Two major factors explain smaller differences between simulation and measurement: First, every measurement system has uncertainties that have to be taken into account; second, every model used for computer simulation is an approximation of complex real-world physical behavior. The authors are keenly interested in the simulation model, and inaccuracies of the simulation could be reduced by making a model physically more realistic. In this particular case, the model has not considered residual stresses and varying material thickness of the part that might be present in a manufactured sheet metal part.

Conclusions
This paper has described a solution for computationally determining the GF shape of a measured sheet metal part for quality assurance purposes. Our solution employs an iterative approach to FE simulation. As only input and output node positions of a simulation are needed to be accessible, so one can implement our method with commercially available simulation software. The paper provides a detailed description of our algorithm and test implementation, using commercial simulation software. Furthermore, the implementation was validated by discussing a comparison of simulation and experimental data. Also, different geometries and materials were used for the validation to show the usability and performance of the method. The quantitative numerical results of our validations show deviations below the uncertainty of typical measurement devices used for assessing sheet metal part geometries, which is ±0.05 mm. Thus, the results produced by the presented method and its implementation meet the quality goals the authors had established for the method. In addition, the geometrical precision and run times of our method were analyzed and discussed for different problem sizes.
To improve run time behavior, it would be worthwhile to analyze the convergence behavior of the algorithm and consider, for example, an extrapolation scheme for faster convergence. Another improvement would be the calculation of the displacement direction used in Loop 1. Also, it would be very beneficial to check whether the method is capable of calculating the GF shape of an overconstrained measurement when the force distribution on the supports is known.
-Algorithm to calculate the gravity-free shape of a simulation model that is suitable for almost all commercial FE solvers. -Validations with simulation data for a real-world use case.
-Experimental validation with a simple, easy to understand sheet metal part.
Funding This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -252408385 -IRTG 2057. Open Access funding enabled and organized by Projekt DEAL.

Availability of data and material On request
Code availability The method proposed in this paper is implemented as a custom code plugin into the commercial software "Abaqus" from Dassault Systèmes.

Declarations
Competing Interests The authors declare that they have no known competing nancial interests or personal relationships that could have appeared to inuence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.