Case study on the influence of kinematic hardening within a parameter-free and non-invasive form finding approach

Inverse form finding – as a type of shape optimization – aims in determining the optimal preform design of a workpiece for a specific forming process, whereby the desired target geometry is known. Recently, a novel parameter-free and heuristic approach was developed to tackle this nonlinear optimization problem. Benchmark tests already delivered promising results. As a particular note-worthy feature of the approach, a coupling to an arbitrary finite element software is feasible in a non-invasive fashion. The aim of this contribution is to investigate the effect of kinematic hardening and cyclic loading on the convergence behavior of the algorithm.


Introduction and related works
Within the metal forming industry, trial-and-error methods of earlier days have since long been replaced by simulative predictions and simulation based optimization strategies. This leads not only to a cost reduction and an acceleration in the development process of new products, but is nowadays also a cornerstone for the manufacturing of highest quality and near-net-shaped functional components. Optimization in order to improve product properties is thereby often the main purpose of numerical forming simulations and consequently a key point in the field of future research, see Andrietti et al. [1]. Thereby, the main issue of optimization, besides the inverse problem of parameter identification, is the determination of an optimal tool and workpiece design, see Chenot et al. [2].
The latter constitutes the investigated type of optimization in this contribution and is also frequently denoted as inverse form finding. It aims in determining the optimal preform design of a workpiece for a specific forming process, whereby the desired target shape is known a priori. Especially in the context of sheet and bulk metal forming, or forging respectively, several concepts have been developed since decades to solve this non-linear optimization problem. Since path-dependency is occurring in plasticity, a similar iterative procedure of all shape optimization approaches can be recognized. It consists of updating the geometry of the material space, before re-starting a new simulation and comparing again the forming results of the spatial space to the desired target shape.
Form finding strategies are commonly distinguished between parameter-free and parameter-based approaches. In the former, a computer aided design (CAD) geometry is given and its control points serve as design variables, as initially proposed by Braibant and Fleury [3]. For example, Fourment and Chenot [4] already proposed a parameter-based optimization approach for determining an optimal workpiece and die design in non-steady-state forming processes. Within parameter-free approaches, as considered in this contribution, the geometry is given in a discretized setting as a finite element (FE) model and its nodal positions serve as design variables. This leads to a larger number of design parameters. Furthermore, nodalbased mesh moving strategies often suffer from distorted elements when updating the material configuration so that regularization and filter techniques have to be deployed, see Bletzinger et al. [5] and Scherer et al. [6].
As a further distinction, different strategies can be found for sheet and for bulk metal forming. Within sheet metal forming, the sought optimal material configuration renders commonly a planar sheet pre-cut and the aim is to determine an optimal cut-out contour, e.g. described by spline functions as in Kim et al. [7]. Before re-starting the forming simulation in the iterative optimization procedure, the surface inside the cut-out contour is automatically re-meshed in most cases, which is a feasible task for two-dimensional structural elements like shells. A different strategy which repositions the nodal coordinates of a discretized planar sheet pre-form without remeshing is presented by the inverse approach of Gou and Batoz [8]. In contrast, within bulk metal forming simulations three dimensional geometries can only be remeshed automatically with great difficulties, so parameter-free mesh moving strategies are predominated in this context.
It should be emphasized that several types of sensitivity analyses (see Haftka and Adelmann [9] for an overview) exist within both, parameter-based and parameter-free approaches, and have already been successfully adopted to the context of inverse form finding. Michaleris and Tortorelli [10] for example provided a general framework for the calculation of elasto-plastic tangent operators or Archajee and Zabaras [11] already presented a sensitivity approach for three-dimensional contact problems. However, in most cases it is cumbersome to compute the required sensitivity derivations of the system responses.
A further interesting and promising approach is proposed by Germain et al. [12] for elasto-plastic problems. Here an inverse mechanical formulation, parametrized with spatial coordinates, is solved with respect to the material coordinates, as introduced by Govindjee et al. [13] for elasticity. For the extension to plasticity, the internal plastic variables are recursively transferred to circumvent the path-dependency.
Recently, a novel heuristic and parameter-free approach was developed for inverse form finding problems, see Landkammer and Steinmann [14]. The updating operation is thereby performed by a reverse tangent mapping after analyzing the nodal distances to the target shape in the final deformation stage. As a benefit of the algorithm, a coupling by subroutines to any desired external FE code is feasible. Benchmark tests already demonstrated an excellent numerical efficiency. Also the applicability to metal forming processes have been already successfully shown [15]. The main goal of this contribution is to evaluate, whether a sophisticated hardening model or cyclic loading with reverse strain path changes affect the convergence behavior. This is realized by a case study with computational examples.
The paper is structured as follows: Sec. 2 explains the functional principle of the non-invasive form finding algorithm. Afterwards, Sec. 3 reiterates some basics from plasticity regarding isotropic and kinematic hardening, before the benchmark model for the computational examples is presented. A subsequent parameter study regarding the constitutive behavior and the loading path is conducted to demonstrate the applicability to combined hardening and cyclic loading. Finally, Sec. 4 draws a conclusion to summarize the findings.

Functional principle of the iterative update procedure
In this section, the basic framework of the update operation is explained, while a more detailed description is presented in Landkammer and Steinmann [14]. The iterative optimization algorithm is purely based on the nodal positions of the material and the spatial configurations. Almost every FE software provides options to export these data by subroutines. After the determination of improved material nodal positions within the update procedure, they are imported again and the next re-computation of the forming simulation starts automatically until suited spatial nodal positions are computed. Due to the ability of coupling the optimization program as a black box to the user's habitual forming software without interfering with the FE code of the simulation itself, the approach is denoted as non-invasive.  Mathematically, the forming process of a body renders a mapping from the material to the spatial configuration. With respect to a FE simulation, this can be formulated by a map of the coordinates ୦ of a discretized material configuration ℬ ୦ to the coordinates ୦ of the spatial configuration ℬ ୲ ୦ : As an objective for the optimization, a discretized target configuration with spatial coordinates ୲ ୦ has to be defined a priori. The optimal material coordinates are then determined, if the nodal distances between both spatial configurations are close enough for all element nodes ‫ܫ(‬ = 1, ⋯ , ݊ ୬୭ୢୣୱ ). This is evaluated by the objective function ߜ ୰୫ୱ , which renders the root mean square over all spatial nodal distances, whereby ூ refers to the nodal difference vectors in the spatial setting.
For a suited repositioning of the material coordinates, the final deformation state after the forming simulation has to be analyzed. The idea is to exploit the linear tangent map ୦ , which arises from the derivative of the deformation map with respect to the design parameters ୦ , for a node-wise fixed-point iteration to minimize the objective function ߜ ୰୫ୱ ( ୲ ୦ , ୡ୰ ୦ ).
However, the tangent map cannot be computed directly at the nodal positions ூ , but the deformation gradient is only given at the Gauss points. The computation of the deformation gradient is feasible with the Jacobian matrices of the material and the spatial configurations, whereby ܰ ( ) refers to isoparametric shape functions evaluated at the ݆-th Gauss point: To achieve an approximation ෩ of the deformation gradient at the nodal positions, a least squares smoothing is performed, as described by Hinton and Campbell [16]. This smoothing technique is also commonly used in FE codes for mapping the entries of a stress tensor from the Gauss points to the element nodes. With the inverse of the smoothed nodal deformation gradient, a difference vector ூ can be mapped from the spatial configuration ℬ ௧ ୦ to the discretized material configuration ℬ ୦ .
The update of the material configuration is then performed by adding all mapped difference vectors ூ to their material nodal positions: In order to control the update procedure in the case of mesh distortions, a global relaxation parameter α is introduced in Eq. (7). Before importing the updated material nodal positions and restarting the FE simulation, the mesh quality is analyzed by an elementwise evaluation of discrete Jacobian matrices as proposed in Knupp et al. [17]. If needed, the relaxation parameter is then automatically reduced, which helps to render the algorithm more robust.

Example: Tension and compression of a cuboid with a cylindric hole
The investigated model is a cuboid with a cylindric hole in its center and has already been used for benchmark tests within inverse form finding problems in [14] and [18].
The aim in this contribution is to analyze the effect of a material law with kinematic hardening to the convergence behavior of the algorithm. The dual phase steel DP600, which is a typical material in automotive industries or for Sheet-Bulk metal formed components [20] respectively, is used for the simulative investigations.

The material model of DP600 steel with isotropic and kinematic hardening
Some basics from plasticity theory are firstly reiterated. Within an elasto-plastic analysis, a yield function has to be evaluated for the current stress state to prove if plastic flow occurs. Herein, the equivalent stress ߪ ത is compared to the current yield stress ߪ ୷ and the yield condition Φ = 0 has to be satisfied at any time in case of plasticity. The equivalent stress renders e. g. for the commonly used orthotropic Hill yield criterion with the Cauchy stress tensor , the back stress tensor and the fourth order Hill tensor ℍ (respectively ℍ = ॴ ୢୣ୴ ୱ୷୫ for the isotropic v.Mises yield criterion with the symmetric deviatoric fourth order unit tensor ॴ ୢୣ୴ ୱ୷୫ ). For the sake of an accurate modelling of the physical behavior of elastoplastic materials, two main types of hardening are usually considered in the context of metal forming. Firstly, isotropic hardening is used to model an expansion of the yield surface. Thereby, the hardening function ߪ ୷ (߳ ୮ ) of the current yield stress with respect to the current equivalent plastic strain ߳ ୮ is evaluated. A frequently used model for isotropic hardening is for example the Voce type saturation law: Secondly, kinematic hardening is used to model a translation of the yield surface in the stress space. Here, the back stress tensor has to be determined simultaneously by an additional evolution equation: Kinematic hardening has to be considered to capture the so-called Bauschinger effect. This refers to a differing yield stress in case of a reverse strain path change when loading is applied in the opposite direction than before, as illustrated in Figure 2.
The parameters for the DP600 dual phase steel, see Table 1, are determined by Bouvier et al. [19] by weighting the results of several mechanical tests. Oliveira et al. [21] further utilized the parameter set for analyzing the springback prediction of sheets, in which the modelling of the Bauschinger effect has a major impact. As a feature, the parameter sets are thereby separately identified in [19] on the one hand for a pure isotropic hardening ( = ) with the Voce law, see Eq. (10) and on the other hand for a combined hardening model. In the latter, the Voce law is again used for the isotropic hardening and the Frederick-Armstrong law, see Eq. (11), is additionally applied for modelling the kinematic hardening effect. The Bauschinger effect is additionally demonstrated for the material parameters of the isotropic and the combined material law by a cyclic loading (tension / compression / tension) applied to a single element within a displacement controlled computation. In the case of a  Figure 2.

Simulation
For evaluating the influence of an isotropic and a combined material modelling on the convergence behavior of the form finding algorithm, the benchmark model of a cuboid with a cylindric hole is analyzed. Figure 3 depicts the dimension of the workpiece and illustrates the symmetry planes, which corresponds to the boundary conditions. The geometry is discretized in the FE model with 150 solid elements (264 nodes). In the subsequent optimization the material start configuration equals the desired spatial target configuration. The material model of the DP600 steel is adopted according to Table 1. The computation is carried out with multiplicative large strain plasticity with the commercial FE software MSC.MarcMentat.
For a comparison of the convergence behavior, the simulations are performed four times combining both hardening models with two different loadings, as listed in Table 2. Firstly, a monotonic loading of a pure compression (displacement of -4 mm in x-direction, 0-1 s) is applied for the isotropic and the combined material model. Secondly, a cyclic loading by a changing compression and tension (displacements of -4/ 4/ -4 mm in x-direction, 0-5 s) is applied for the hardening models. The according load path of the displacement controlled forming simulation is depicted in Figure 4.  Isotropic Cyclic (0-5 s) 4 Combined Cyclic (0-5 s)

Optimization
The optimization algorithm is implemented in Matlab.
For the update of the material configuration of the workpiece, a constant relaxation parameter (ߙ = 1) is used in all iterations and for all four modelled variants. Figure 5 shows the determined optimal material configuration for the pure compression (monotonic loading, 4 mm in x-direction) with the underlying isotropic (A) and the combined hardening model (B). Inputting the determined optimal material configurations, subsequent forming simulations compute the desired spatial nodal positions of the given target configuration (C, D).
For comparison, Figure 6 shows the determined optimal material configuration for the cyclic loading with reverse strain path changes (-4/ 4/ -4 mm in x-direction) with the underlying isotropic (A) and the combined hardening model (B). Subsequent forming simulations after inputting the determined optimal material configurations compute again the desired spatial nodal positions of the given target configuration (C, D). The equivalent plastic strains ߳ ୮ , the equivalent stresses ߪ ത and the equivalent back stresses ‫ܤ‬ ത are uniformly scaled and depicted. The computed maximal equivalent plastic strains vary between 0.265 and 0.287, whereas the maximal equivalent stresses range from 772 MPa to 848 MPa. In all cases, the algorithm continuously decreases the objective functions ߜ ୰୫ୱ , which can also be interpreted as the nodal standard deviation between the computed and the target configuration.
Furthermore, nearly linear convergence rates can be observed, see Figure 7. The investigation demonstrates that applying a combined hardening model and a cyclic loading does not poses problems to the optimization approach. A minimal error value ߜ ୰୫ୱ close to 10 ିଵହ is achieved after 25 (case 1: isotropic hardening / monotonic loading) to 30 iterations (case 3: isotropic hardening / cyclic loading). In the case of employing the combined hardening model no tendency regarding the efficiency is apparent (slightly more iterations for monotonic loading, but less iterations for cyclic loading).
In contrast, a slight difference can be observed, when applying a cyclic loading (more iterations are needed) in comparison to the case of a monotonic loading. A clear relation between the amount of the equivalent stresses, or equivalent plastic strains respectively, and the convergence rates is not recognizable.  Table 3 lists additionally the error values of the objective function and shows again clearly the excellent efficiency of the non-invasive approach.

Conclusion
Kinematic hardening is of crucial importance when the forming of a workpiece includes strain-path changes, as typically occurring e. g. in bulk forming of a sheet-metal. Therefore, the convergence behavior of a non-invasive iterative approach for inverse form finding, whose operating principle is briefly summarized, is investigated by a case study with a pure isotropic and a combined (isotropic and kinematic) hardening law. It is demonstrated that different hardening models does not affect the performance of the algorithm. Furthermore, only slightly more iterations are needed when applying a cyclic instead of a monotonic loading.