Fig. 1(c) schematically illustrates the two-step identification process introduced in this work. Prior to the FE-VFM analysis, uniaxial tensile tests for the base sheet metals are conducted, and the mechanical properties of the base materials are obtained. Then, virtual work residuals are calculated using the base material's constitutive parameters in each sub-zone, in STEP1. If materials are affected during the FSW, then residuals will be larger compared to those in not affected area, and boundaries of the WAZ can be easily identified. In STEP2, distributions of the hardening law parameters are optimized by minimizing virtual work residuals, for each sub-domain independently. The distributions of constitutive parameters within each sub-zone are quadratically interpolated,

$$ \bar{k} = Q(Y)\bar{k}^N = \begin{bmatrix} 1 & Y & Y^2 \end{bmatrix} \begin{bmatrix} 0 & 1 & 0 \\ \frac{1}{3} & 0 & \frac{1}{3} \\ \frac{2}{9} & -\frac{4}{9} & \frac{2}{9} \end{bmatrix} \begin{bmatrix} \frac{-1}{k} \\ \frac{2}{k} \\ \frac{-3}{k} \end{bmatrix} \quad (4) $$

where $\bar{k}$ is constitutive parameters normalized by its based material values, and Y is coordinate value in the length direction of the specimen which has zero values at the center of each sub-zone, and the nodal position in parameter interpolation is located Y=1.5, 0, and -1.5 for node number (uppercase in Eq. (4)) 1,2, and 3.

<figure><img src="image_1678.png" /><figcaption>Fig. 1. (a) Specimen geometries for applying the FE-VFM. (b) Cubic quadrilateral high-order finite element (FE) meshes used in the FE-VFM. (c) Two-step identification procedure.</figcaption></figure>

The isotropic constitutive model was used in this research, the generalized Hooke's law and the von Mises yield criteria were used. The Swift hardening law was used,

$$ \bar{\sigma} = K(\bar{\varepsilon} + e_0)^n. \quad (5) $$

Material constants of the base material AA6061-T6 sheets are listed in Table 1. In inverse measurements using the FE-VFM, only Swift hardening law parameters were optimized in STEP2, considering that the elastic constants are not significantly changed owing to the thermal histories. Also, note that normalized-nodal constitutive parameters of the Swift law were identified in STEP2 (denoted $\bar{k}$ in Eq. (4)), and real parameter distributions can be described by multiplying the values of the base materials in Table 1 and Eq. (5).
<table><caption>Table 1. Elastic constants and the Swift law parameters of the AA6061-T6 sheets.</caption><thead><tr><th colspan="2">Elastic constants</th><th colspan="3">Swift law parameters</th></tr><tr><th>E [GPa]</th><th>ν</th><th>K [MPa]</th><th>e<sub>0</sub></th><th>n</th></tr></thead><tbody><tr><td>70</td><td>0.33</td><td>508.1</td><td>0.03602</td><td>0.1577</td></tr></tbody></table>

Meanwhile, quadratic interpolation of the constitutive parameters makes it hard to uniquely identify when conventional types of virtual fields are employed. It is because of that the three nodal-normalized constitutive material parameters are interconnected quadratically in each sub-zone, and a wrong combination of the nodal values can give a false zero residuals. A special type of virtual field is introduced to avoid this issue, based on the Gaussian normal distribution function. The virtual position field is defined as,

$$ \delta x_1 = \begin{cases} X_1(Y^2 - 2.25)(|Y| \le 1.5) \\ 0(|Y| > 1.5) \end{cases}, \quad \delta x_2 = \begin{cases} \Phi_{i=1,2,3}(Y)(|Y| \le 1.5) \\ 0(|Y| > 1.5) \end{cases} \qquad (6) $$

where $\Phi_{i=1,2,3}$ is the cumulative distribution function,

$$ \Phi_{i=1,2,3}(Y) = \frac{1}{2} \left\{ 1 + \operatorname{erf}\left(\frac{Y - \mu_i}{\sigma_N \sqrt{2}}\right) \right\}. \qquad (7) $$

Then, the derivatives of the virtual position, or the virtual deformation gradient is given as

$$ \delta F = \frac{\partial \delta x}{\partial x} = \begin{pmatrix} Y^2 - 2.25 & 2X_1 Y \\ 0 & \varphi_{i=1,2,3}(Y) \end{pmatrix} \qquad (8) $$

where $\varphi_{i=1,2,3}$ is the Gaussian normal distribution function,

$$ \varphi_{i=1,2,3}(Y) = \frac{1}{\sigma_N \sqrt{2\pi}} \exp \left\{ -\frac{1}{2} \left( \frac{Y - \mu_i}{\sigma_N} \right)^2 \right\}. \qquad (9) $$

The Gaussian distribution can act as an amplifier of the strain (or stress) data within a target range in the sub-zones, when different $\sigma_N$ and $\mu_i$ are adopted. In this work, three independent virtual fields with values of $\sigma_N = 0.5$ and $\mu_i = 1, 0, -1$ for i=1,2,3 were used. Fig. 2 schematically illustrates the role of the Gaussian distribution type virtual fields. The $\sigma_N$ in Gaussian distribution determines the width of a larger value area, $2\sigma$. Considering that the length of the sub-zone is 3 mm, $\sigma_N = 0.5$ is used. Also, $\mu_i$ determines the peak location in the Gaussian distribution, and three different $\mu_i = 1, 0, -1$ are applied. Consequently, three independent virtual fields amplify the mechanical data of the top, middle, and bottom regions, selectively. Owing to these three Gaussian distribution-type virtual fields, three nodal-normalized constitutive parameters can be uniquely optimized despite that these parameters are interconnected with quadratic interpolation functions.
<table><caption>Table 1. Elastic constants and the Swift law parameters of the AA6061-T6 sheets.</caption><thead><tr><th colspan="2">Elastic constants</th><th colspan="3">Swift law parameters</th></tr><tr><th><i>E</i> [GPa]</th><th><i>ν</i></th><th><i>K</i> [MPa]</th><th><i>e</i><sub>0</sub></th><th><i>n</i></th></tr></thead><tbody><tr><td>70</td><td>0.33</td><td>508.1</td><td>0.03602</td><td>0.1577</td></tr></tbody></table>

Meanwhile, quadratic interpolation of the constitutive parameters makes it hard to uniquely identify when conventional types of virtual fields are employed. It is because of that the three nodal-normalized constitutive material parameters are interconnected quadratically in each sub-zone, and a wrong combination of the nodal values can give a false zero residuals. A special type of virtual field is introduced to avoid this issue, based on the Gaussian normal distribution function. The virtual position field is defined as,

$$ \delta x_1 = \begin{cases} X_1(Y^2 - 2.25)(|Y| \le 1.5) \\ 0(|Y| > 1.5) \end{cases}, \quad \delta x_2 = \begin{cases} \Phi_{i=1,2,3}(Y)(|Y| \le 1.5) \\ 0(|Y| > 1.5) \end{cases} \qquad (6) $$

where $\Phi_{i=1,2,3}$ is the cumulative distribution function,

$$ \Phi_{i=1,2,3}(Y) = \frac{1}{2} \left\{ 1 + \operatorname{erf}\left(\frac{Y - \mu_i}{\sigma_N \sqrt{2}}\right) \right\}. \qquad (7) $$

Then, the derivatives of the virtual position, or the virtual deformation gradient is given as

$$ \delta \mathbf{F} = \frac{\partial \delta x}{\partial x} = \begin{pmatrix} Y^2 - 2.25 & 2X_1 Y \\ 0 & \varphi_{i=1,2,3}(Y) \end{pmatrix} \qquad (8) $$

where $\varphi_{i=1,2,3}$ is the Gaussian normal distribution function,

$$ \varphi_{i=1,2,3}(Y) = \frac{1}{\sigma_N \sqrt{2\pi}} \exp \left\{ -\frac{1}{2} \left( \frac{Y - \mu_i}{\sigma_N} \right)^2 \right\}. \qquad (9) $$

The Gaussian distribution can act as an amplifier of the strain (or stress) data within a target range in the sub-zones, when different $\sigma_N$ and $\mu_i$ are adopted. In this work, three independent virtual fields with values of $\sigma_N = 0.5$ and $\mu_i = 1, 0, -1$ for $i=1,2,3$ were used. Fig. 2 schematically illustrates the role of the Gaussian distribution type virtual fields. The $\sigma_N$ in Gaussian distribution determines the width of a larger value area, $2\sigma$. Considering that the length of the sub-zone is 3 mm, $\sigma_N = 0.5$ is used. Also, $\mu_i$ determines the peak location in the Gaussian distribution, and three different $\mu_i = 1, 0, -1$ are applied. Consequently, three independent virtual fields amplify the mechanical data of the top, middle, and bottom regions, selectively. Owing to these three Gaussian distribution-type virtual fields, three nodal-normalized constitutive parameters can be uniquely optimized despite that these parameters are interconnected with quadratic interpolation functions.
<figure><img src="image_1680.png" /><figcaption>Fig. 2. The Gaussian normal distribution type virtual fields for amplifying data in different regions in the sub-zone.</figcaption></figure>

## Materials and experiments

The friction stir welded AA6061-T6 sheets were investigated. The thickness of the base material sheets was 3 mm. Rectangular workpieces of 300 mm x 150 mm were joined into 300 mm size foursquare sheets. The FSW conditions were as follows. The pin of the FSW tool was truncated cone shape, with a smaller bottom diameter and larger top diameter of 4 and 6 mm. The diameter of the FSW tool shoulder was 12.72 mm. The tool travel speed was 500 mm/min. The FSW was conducted with two conditions, with tool rotation speeds of 1200 and 1600 RPM.

The tensile specimens for applying the FE-VFM were made by wire electric discharge machining (W-EDM) following the geometries in Fig. 1(a). The FSW process remains a trace of the tool path like a plumbing groove. So, the surface of the tensile specimen was

trimmed by the W-EDM, and the thickness of the specimen was reduced to 2.3 mm. For the validation of the measured local strain hardening distributions, a different configuration of tensile tests was conducted, which loaded along the weld line direction, as shown in Fig. 3. These test results were compared with the FEM predictions conducted with the inverse measurements.

All the tensile tests were conducted with the same condition. Instron 8801 100 kN UTM was used, and the test speed was 1 mm/min, which is within quasi-static conditions. Deformations were measured by the commercial DIC packages, the VIC-3D. In the validation tensile tests, displacements were measured by the virtual extensometer with a gauge length of 40 mm in DIC post-analysis, and the force-displacement curves measured by this were compared with the validation FEM predictions.

<figure><img src="image_1681.png" /><figcaption>Fig. 3. Schematic illustration of the validation tensile test specimens.</figcaption></figure>

## Result and discussion

Fig. 4(a) shows the virtual work residuals (VWR) for each sub-zone in STEP1, which were calculated using the base material properties in Table 1. The VWR shows considerably larger values from the -12 to 9 mm range for the 1200 RPM case, and the -12 to 12 mm range for the
<figure><img src="image_1679.png" /><figcaption>Fig. 2. The Gaussian normal distribution type virtual fields for amplifying data in different regions in the sub-zone.</figcaption></figure>

## Materials and experiments

The friction stir welded AA6061-T6 sheets were investigated. The thickness of the base material sheets was 3 mm. Rectangular workpieces of 300 mm x 150 mm were joined into 300 mm size foursquare sheets. The FSW conditions were as follows. The pin of the FSW tool was truncated cone shape, with a smaller bottom diameter and larger top diameter of 4 and 6 mm. The diameter of the FSW tool shoulder was 12.72 mm. The tool travel speed was 500 mm/min. The FSW was conducted with two conditions, with tool rotation speeds of 1200 and 1600 RPM.

The tensile specimens for applying the FE-VFM were made by wire electric discharge machining (W-EDM) following the geometries in Fig. 1(a). The FSW process remains a trace of the tool path like a plumbing groove. So, the surface of the tensile specimen was

trimmed by the W-EDM, and the thickness of the specimen was reduced to 2.3 mm. For the validation of the measured local strain hardening distributions, a different configuration of tensile tests was conducted, which loaded along the weld line direction, as shown in Fig. 3. These test results were compared with the FEM predictions conducted with the inverse measurements.

All the tensile tests were conducted with the same condition. Instron 8801 100 kN UTM was used, and the test speed was 1 mm/min, which is within quasi-static conditions. Deformations were measured by the commercial DIC packages, the VIC-3D. In the validation tensile tests, displacements were measured by the virtual extensometer with a gauge length of 40 mm in DIC post-analysis, and the force-displacement curves measured by this were compared with the validation FEM predictions.

<figure><img src="image_1680.png" /><figcaption>Fig. 3. Schematic illustration of the validation tensile test specimens.</figcaption></figure>

## Result and discussion

Fig. 4(a) shows the virtual work residuals (VWR) for each sub-zone in STEP1, which were calculated using the base material properties in Table 1. The VWR shows considerably larger values from the -12 to 9 mm range for the 1200 RPM case, and the -12 to 12 mm range for the
1600 mm case. These regions are the WAZ, and the outside regions are the base materials unaffected by the FSW process. This result indicates that the WAZ is smaller for the 1200 RPM case compared to the 1600 RPM case, on the retreating side (RS). The RS is a side where the direction of rotation and tool translation direction are opposite, and heat generated during the FSW is smaller compared to the advancing side (AS). For the 1200 RPM, less frictional heat would generated compared to the 1600 RPM, and the smaller WAZ can be formed.

<figure><img src="image_1682.png" /><figcaption>Fig. 4. Inverse identification results. (a) Virtual work residuals for each sub-zone in the WAZ boundaries identification (in STEP1). (b) Distributions of the normal-normalized constitutive parameters in the Swift laws (in STEP2).</figcaption></figure>

Fig. 4(b) shows distributions of the constitutive parameters in the Swift laws. Note that solid lines are the quadratic interpolations of the normalized values by base material values, and the square dots denote nodal-normalized values which are the 1st direct inverse measurements in this work. It can be seen that the parameter n values within the -6 to 6 mm range exhibit almost twice larger than the base material value for both 1200 and 1600 RPM cases. The Swift law parameter n is highly related to the uniform elongation by the Considère condition. Also, this region is exactly matched with the diameter of the FSW tool shoulder. Hence, this result indicates that welded materials were softened by the FSW process. Also, rapid drops of the n values are found just near and outside of this range, which corresponds to the thermo-mechanically affected zone (TMAZ). A lower n value indicates a shorter uniform elongation in the tensile tests, and the n value was only around 40% compared to the base material properties. The TMAZ is typically considered to be a weak location in the FSWed sheets, hence, the result of sudden valley-shape distributions of n value in this area can be interpreted that the lower uniform elongation was measured in these relatively narrow TMAZ regions.

Fig. 5 shows the force-displacement curves in the validation tensile tests loaded in weld line directions. The FEM analyses were conducted by using the distributions of the local strain hardening in Fig. 4(a). It can be noticed that the experimental values (empty circle dots) and FEM predictions (solid lines) show very good matches. This confirms that the suggested FE-VFM method can measure local strain hardening distributions of the welded sheets with good accuracy.
1600 mm case. These regions are the WAZ, and the outside regions are the base materials unaffected by the FSW process. This result indicates that the WAZ is smaller for the 1200 RPM case compared to the 1600 RPM case, on the retreating side (RS). The RS is a side where the direction of rotation and tool translation direction are opposite, and heat generated during the FSW is smaller compared to the advancing side (AS). For the 1200 RPM, less frictional heat would generated compared to the 1600 RPM, and the smaller WAZ can be formed.

<figure><img src="image_1681.png" /><figcaption>Fig. 4. Inverse identification results. (a) Virtual work residuals for each sub-zone in the WAZ boundaries identification (in STEP1). (b) Distributions of the normal-normalized constitutive parameters in the Swift laws (in STEP2).</figcaption></figure>

Fig. 4(b) shows distributions of the constitutive parameters in the Swift laws. Note that solid lines are the quadratic interpolations of the normalized values by base material values, and the square dots denote nodal-normalized values which are the 1st direct inverse measurements in this work. It can be seen that the parameter n values within the -6 to 6 mm range exhibit almost twice larger than the base material value for both 1200 and 1600 RPM cases. The Swift law parameter n is highly related to the uniform elongation by the Considère condition. Also, this region is exactly matched with the diameter of the FSW tool shoulder. Hence, this result indicates that welded materials were softened by the FSW process. Also, rapid drops of the n values are found just near and outside of this range, which corresponds to the thermo-mechanically affected zone (TMAZ). A lower n value indicates a shorter uniform elongation in the tensile tests, and the n value was only around 40% compared to the base material properties. The TMAZ is typically considered to be a weak location in the FSWed sheets, hence, the result of sudden valley-shape distributions of n value in this area can be interpreted that the lower uniform elongation was measured in these relatively narrow TMAZ regions.

Fig. 5 shows the force-displacement curves in the validation tensile tests loaded in weld line directions. The FEM analyses were conducted by using the distributions of the local strain hardening in Fig. 4(a). It can be noticed that the experimental values (empty circle dots) and FEM predictions (solid lines) show very good matches. This confirms that the suggested FE-VFM method can measure local strain hardening distributions of the welded sheets with good accuracy.
<figure><img src="image_1683.png" /><figcaption>Fig. 5. Tensile tests and FEM analyses for the validation of the inverse measurements.</figcaption></figure>

**Summary**

Inverse identification of the local strain hardening distributions for the FSWed AA6061-T6 sheets was conducted by the FE-VFM. In the FE-VFM, new numerical schemes were introduced, including the two-step identification procedure, the sub-zone independent identification with quadratic interpolation of the constitutive parameters, and the new Guass distribution function type virtual fields. Measured distributions of local strain hardening well capture a material softening due to the FSW process, and the weaker material properties in the TMAZ. Also, experimental force-displacement curves and the FEM predictions using the inverse measurements showed good matches in the validation tensile tests, which confirms the capability of the proposed inverse identification method.

**Acknowledgment**

This work was supported by the Fundamental Research Program of the Korea Institute of Materials Science (KIMS, PNK9800) and the National Research Council of Science and Technology (NST) grant by the Korea government (MSIT) (CRC23011-210).

**References**

[1] X. He, F. Gu, A. Ball, A review of numerical analysis of friction stir welding, Prog. Mat. Sci. 65 (2014) 1-66. https://doi.org/10.1016/j.pmatsci.2014.03.003

[2] D. Myung, W. Noh, J.-H. Kim, J. Kong, S.-T. Hong, M.-G. Lee, Probing the Mechanism of Friction Stir Welding with ALE Based Finite Element Simulations and Its Application to Strength Prediction of Welded Aluminum, Metals Mat. Int. 27 (2021) 650-666. https://doi.org/10.1007/s12540-020-00901-8

[3] G. Güleryüz, Relationship between FSW parameters and hardness of the ferritic steel joints: Modeling and optimization, Vacuum 178 (2020) 109449. https://doi.org/10.1016/j.vacuum.2020.109449

[4] J.H. Kim, F. Barlat, C. Kim, K. Chung, Themo-mechanical and microstructural modeling of friction stir welding of 6111-T4 aluminum alloys, Metals Mat. Int. 15 (2009) 125-132. https://doi.org/10.1007/s12540-009-0125-5

[5] D. Rao, J. Heerens, G. Alves Pinheiro, J.F. dos Santos, N. Huber, On characterisation of local stress-strain properties in friction stir welded aluminium AA5083 sheets using micro-tensile specimen testing and instrumented indentation technique, Mat. Sci. Eng. A 527 (2010) 5018-5025. https://doi.org/10.1016/j.msea.2010.04.047
<figure><img src="image_1682.png" /><figcaption>Fig. 5. Tensile tests and FEM analyses for the validation of the inverse measurements.</figcaption></figure>

**Summary**

Inverse identification of the local strain hardening distributions for the FSWed AA6061-T6 sheets was conducted by the FE-VFM. In the FE-VFM, new numerical schemes were introduced, including the two-step identification procedure, the sub-zone independent identification with quadratic interpolation of the constitutive parameters, and the new Guass distribution function type virtual fields. Measured distributions of local strain hardening well capture a material softening due to the FSW process, and the weaker material properties in the TMAZ. Also, experimental force-displacement curves and the FEM predictions using the inverse measurements showed good matches in the validation tensile tests, which confirms the capability of the proposed inverse identification method.

**Acknowledgment**

This work was supported by the Fundamental Research Program of the Korea Institute of Materials Science (KIMS, PNK9800) and the National Research Council of Science and Technology (NST) grant by the Korea government (MSIT) (CRC23011-210).

**References**

[1] X. He, F. Gu, A. Ball, A review of numerical analysis of friction stir welding, Prog. Mat. Sci. 65 (2014) 1-66. https://doi.org/10.1016/j.pmatsci.2014.03.003

[2] D. Myung, W. Noh, J.-H. Kim, J. Kong, S.-T. Hong, M.-G. Lee, Probing the Mechanism of Friction Stir Welding with ALE Based Finite Element Simulations and Its Application to Strength Prediction of Welded Aluminum, Metals Mat. Int. 27 (2021) 650-666. https://doi.org/10.1007/s12540-020-00901-8

[3] G. Güleryüz, Relationship between FSW parameters and hardness of the ferritic steel joints: Modeling and optimization, Vacuum 178 (2020) 109449. https://doi.org/10.1016/j.vacuum.2020.109449

[4] J.H. Kim, F. Barlat, C. Kim, K. Chung, Themo-mechanical and microstructural modeling of friction stir welding of 6111-T4 aluminum alloys, Metals Mat. Int. 15 (2009) 125-132. https://doi.org/10.1007/s12540-009-0125-5

[5] D. Rao, J. Heerens, G. Alves Pinheiro, J.F. dos Santos, N. Huber, On characterisation of local stress-strain properties in friction stir welded aluminium AA5083 sheets using micro-tensile specimen testing and instrumented indentation technique, Mat. Sci. Eng. A 527 (2010) 5018-5025. https://doi.org/10.1016/j.msea.2010.04.047
[6] D. Bernard, D.G. Hattingh, W.E. Goosen, M.N. James, High Speed Friction Stir Welding of 5182-H111 Alloy: Temperature and Microstructural Insights into Deformation Mechanisms, Metals Mat. Int. 27 (2021) 2821-2836. https://doi.org/10.1007/s12540-020-00622-y

[7] G. Cao, M. Ren, Y. Zhang, W. Peng, T. Li, A partitioning method for friction stir welded joint of AA2219 based on tensile test, Metals 10 (2020) 65. https://doi.org/10.3390/met10010065

[8] M.A. Sutton, J.H. Yan, S. Avril, F. Pierron, S.M. Adeeb, Identification of heterogeneous constitutive parameters in a welded specimen: Uniform stress and virtual fields methods for material property estimation, Exp. Mech. 48 (2008) 451-464. https://doi.org/10.1007/s11340-008-9132-6

[9] F. Tucci, A. Andrade-Campos, S. Thuillier, P. Carlone, Calibration of the Elasto-Plastic Properties of Friction Stir Welded Blanks in Aluminum Alloy AA6082, Key Eng. Mater. 926 (2022) 2183-2192. https://doi.org/10.4028/p-3a8e45

[10] G. Le Louëdec, F. Pierron, M.A. Sutton, A.P. Reynolds, Identification of the Local Elasto-Plastic Behavior of FSW Welds Using the Virtual Fields Method, Exp. Mech. 53 (2013) 849-859. https://doi.org/10.1007/s11340-012-9679-0

[11] G. Le Louëdec, F. Pierron, M.A. Sutton, C. Siviour, A.P. Reynolds, Identification of the Dynamic Properties of Al 5456 FSW Welds Using the Virtual Fields Method, J. Dynamic Behavior Mater. 1 (2015) 176-190. https://doi.org/10.1007/s40870-015-0014-6

[12] C. Kim, J.H. Kim, M.-G. Lee, Identification of Inhomogeneous Plastic Constitutive Models of Friction Stir Welded Aluminum Alloy Sheets Using Virtual Fields Method, Conf. Proc. Soci. Exp. Mech. (2020) 157-161. https://doi.org/10.1007/978-3-030-30098-2_24

[13] H. Jiang, Z. Lei, R. Bai, W. Wu, et al., Identifying elasto-plastic damage coupling model of laser-welded aluminum alloy by virtual field method and digital image correlation, Opt. Laser Tech. 129 (2020) 106268. https://doi.org/10.1016/j.optlastec.2020.106268

[14] C. Kim, D. Myung, D. Kim, M.-G. Lee, Inhomogeneous flow stresses in FSW jointed aluminum alloy sheets inversely identified by FE-VFM, Int. J Mech. Sci. 245 (2023) 108097. https://doi.org/10.1016/j.ijmecsci.2022.108097

[15] C. Kim, M.-G. Lee, Finite element-based virtual fields method with pseudo-real deformation fields for identifying constitutive parameters, Int. J of Sol. Struct. 233 (2021) 111204. https://doi.org/10.1016/j.ijsolstr.2021.111204

[16] C. Kim, J.H. Kim, M.-G. Lee, A virtual fields method for identifying anisotropic elastic constants of fiber reinforced composites using a single tension test: Theory and validation, Composites Part B: Eng. 200 (2020) 108338. https://doi.org/10.1016/j.compositesb.2020.108338
[6] D. Bernard, D.G. Hattingh, W.E. Goosen, M.N. James, High Speed Friction Stir Welding of 5182-H111 Alloy: Temperature and Microstructural Insights into Deformation Mechanisms, Metals Mat. Int. 27 (2021) 2821-2836. https://doi.org/10.1007/s12540-020-00622-y

[7] G. Cao, M. Ren, Y. Zhang, W. Peng, T. Li, A partitioning method for friction stir welded joint of AA2219 based on tensile test, Metals 10 (2020) 65. https://doi.org/10.3390/met10010065

[8] M.A. Sutton, J.H. Yan, S. Avril, F. Pierron, S.M. Adeeb, Identification of heterogeneous constitutive parameters in a welded specimen: Uniform stress and virtual fields methods for material property estimation, Exp. Mech. 48 (2008) 451-464. https://doi.org/10.1007/s11340-008-9132-6

[9] F. Tucci, A. Andrade-Campos, S. Thuillier, P. Carlone, Calibration of the Elasto-Plastic Properties of Friction Stir Welded Blanks in Aluminum Alloy AA6082, Key Eng. Mater. 926 (2022) 2183-2192. https://doi.org/10.4028/p-3a8e45

[10] G. Le Louëdec, F. Pierron, M.A. Sutton, A.P. Reynolds, Identification of the Local Elasto-Plastic Behavior of FSW Welds Using the Virtual Fields Method, Exp. Mech. 53 (2013) 849-859. https://doi.org/10.1007/s11340-012-9679-0

[11] G. Le Louëdec, F. Pierron, M.A. Sutton, C. Siviour, A.P. Reynolds, Identification of the Dynamic Properties of Al 5456 FSW Welds Using the Virtual Fields Method, J. Dynamic Behavior Mater. 1 (2015) 176-190. https://doi.org/10.1007/s40870-015-0014-6

[12] C. Kim, J.H. Kim, M.-G. Lee, Identification of Inhomogeneous Plastic Constitutive Models of Friction Stir Welded Aluminum Alloy Sheets Using Virtual Fields Method, Conf. Proc. Soci. Exp. Mech. (2020) 157-161. https://doi.org/10.1007/978-3-030-30098-2_24

[13] H. Jiang, Z. Lei, R. Bai, W. Wu, et al., Identifying elasto-plastic damage coupling model of laser-welded aluminum alloy by virtual field method and digital image correlation, Opt. Laser Tech. 129 (2020) 106268. https://doi.org/10.1016/j.optlastec.2020.106268

[14] C. Kim, D. Myung, D. Kim, M.-G. Lee, Inhomogeneous flow stresses in FSW jointed aluminum alloy sheets inversely identified by FE-VFM, Int. J Mech. Sci. 245 (2023) 108097. https://doi.org/10.1016/j.ijmecsci.2022.108097

[15] C. Kim, M.-G. Lee, Finite element-based virtual fields method with pseudo-real deformation fields for identifying constitutive parameters, Int. J of Sol. Struct. 233 (2021) 111204. https://doi.org/10.1016/j.ijsolstr.2021.111204

[16] C. Kim, J.H. Kim, M.-G. Lee, A virtual fields method for identifying anisotropic elastic constants of fiber reinforced composites using a single tension test: Theory and validation, Composites Part B: Eng. 200 (2020) 108338. https://doi.org/10.1016/j.compositesb.2020.108338
# Identifiability analysis of material identification using nonlinear VFM

STARMAN Bojan¹,²,a, ZHANG Yi²,b, LAVA Pascal³,c, HALILOVIČ Miroslav¹,d, and
COPPIETERS Sam²,e*

¹University of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, 1000 Ljubljana, Slovenia

²KU Leuven, Department of Materials Engineering, Ghent Campus, 9000 Ghent, Belgium

³MatchID NV, Ghent, Belgium

abojan.starman@fs.uni-lj.si, byi.zhang@kuleuven.be, cpascal.lava@matchid.eu,

dmiroslav.halilovic@fs.uni-lj.si, esam.coppieters@kuleuven.be

**Keywords:** Inverse Identification, Digital Image Correlation, Virtual Field Method, Anisotropy

**Abstract.** The success of inverse material model identification depends on the interaction between the adopted material model, the design of the heterogeneous specimens, the quality of the full-field measurements and the employed inverse identification method. Although inverse identification with full fields usually uses either FEMU or nonlinear VFM algorithms, a range of specimen designs and heterogeneity indicators have been proposed to assess the quality of the measured field and specimen design. While many studies investigate the effects of strain field heterogeneity on material model identification, few of them address the comprehensive interaction of all the above features and investigate their interactions during inverse identification through identifiability analysis. In this study, we analyze the identifiability of the parameters of the YLD2000-2d model used to describe the plastic anisotropy of steel sheet DC04 using a perforated biaxial specimen with the nonlinear VFM method. For this purpose, we performed a virtual DIC experiment with known material parameters by simulating the test in ABAQUS/Standard, generating synthetic images and reconstructing the strains via stereo DIC. Before inverse identification with a nonlinear sensitivity-based VFM, we analyzed the sensitivity of the virtual work to parameter changes and performed an identifiability analysis.

**Introduction**

The reliability of numerical analyses in structural mechanics depends largely on the accurate modeling of materials. For metallic materials, the models are usually formulated in the form of complex algebraic differential equations, especially when a comprehensive representation of advanced material properties is desired. In addition to the formulation itself, these models include a number of parameters that must be determined by appropriate material tests.

Particularly when characterizing sheet metal, the free parameters of the models are traditionally determined using uniaxial standard tests. Usually, only certain parameters can be determined with analytical expressions derived from statically determined test configurations. Consequently, these tests are designed to measure homogeneous deformations determined by electrical strain gages or extensometers, which easily correlate with the global force measured by the load cell of the testing machine. This approach is commonly referred to as Materials testing 1.0 (MT1.0) [1].

Alternatively, the advent of optical technologies has opened new possibilities for characterizing material behavior. The use of 3D scanners and cameras enables the precise measurement of heterogeneous body displacements and deformations under load. The ability to accurately track and quantify displacements and strains enables engineers to assess material behavior and make informed design decisions.

<figure><img src="image_1684.png" /></figure>
# Identifiability analysis of material identification using nonlinear VFM

STARMAN Bojan¹,²,a, ZHANG Yi²,b, LAVA Pascal³,c, HALILOVIČ Miroslav¹,d, and
COPPIETERS Sam²,e*

¹University of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, 1000 Ljubljana, Slovenia

²KU Leuven, Department of Materials Engineering, Ghent Campus, 9000 Ghent, Belgium

³MatchID NV, Ghent, Belgium

abojan.starman@fs.uni-lj.si, byi.zhang@kuleuven.be, cpascal.lava@matchid.eu,

dmiroslav.halilovic@fs.uni-lj.si, esam.coppieters@kuleuven.be

**Keywords:** Inverse Identification, Digital Image Correlation, Virtual Field Method, Anisotropy

**Abstract.** The success of inverse material model identification depends on the interaction between the adopted material model, the design of the heterogeneous specimens, the quality of the full-field measurements and the employed inverse identification method. Although inverse identification with full fields usually uses either FEMU or nonlinear VFM algorithms, a range of specimen designs and heterogeneity indicators have been proposed to assess the quality of the measured field and specimen design. While many studies investigate the effects of strain field heterogeneity on material model identification, few of them address the comprehensive interaction of all the above features and investigate their interactions during inverse identification through identifiability analysis. In this study, we analyze the identifiability of the parameters of the YLD2000-2d model used to describe the plastic anisotropy of steel sheet DC04 using a perforated biaxial specimen with the nonlinear VFM method. For this purpose, we performed a virtual DIC experiment with known material parameters by simulating the test in ABAQUS/Standard, generating synthetic images and reconstructing the strains via stereo DIC. Before inverse identification with a nonlinear sensitivity-based VFM, we analyzed the sensitivity of the virtual work to parameter changes and performed an identifiability analysis.

**Introduction**

The reliability of numerical analyses in structural mechanics depends largely on the accurate modeling of materials. For metallic materials, the models are usually formulated in the form of complex algebraic differential equations, especially when a comprehensive representation of advanced material properties is desired. In addition to the formulation itself, these models include a number of parameters that must be determined by appropriate material tests.

Particularly when characterizing sheet metal, the free parameters of the models are traditionally determined using uniaxial standard tests. Usually, only certain parameters can be determined with analytical expressions derived from statically determined test configurations. Consequently, these tests are designed to measure homogeneous deformations determined by electrical strain gages or extensometers, which easily correlate with the global force measured by the load cell of the testing machine. This approach is commonly referred to as Materials testing 1.0 (MT1.0) [1].

Alternatively, the advent of optical technologies has opened new possibilities for characterizing material behavior. The use of 3D scanners and cameras enables the precise measurement of heterogeneous body displacements and deformations under load. The ability to accurately track and quantify displacements and strains enables engineers to assess material behavior and make informed design decisions.

<figure><img src="image_1683.png" /></figure>
Digital image correlation (DIC) [2] is one of the most commonly used techniques for measuring strain fields. In this approach, sequentially captured images are processed to extract valuable information that facilitates the study of both local behavior and global responses in structures. In addition, DIC can theoretically operate in real time and provides the ability to monitor structural responses under dynamic loading conditions. This feature makes DIC a powerful tool for various applications ranging from material testing to structural analysis and product development.

These novel technologies unleashed a potential for development of novel material calibration strategies, where a richness of full-field measurement of displacements, strains and eventually stresses over a complete Region of interest (ROI) can be analyzed and employed for material calibration. This triggered also a reconsideration of test specimen designs which are needed to be tailored for novel optical measurement technologies to fully exploit the potential, resulting in a new paradigm in material testing using complex tests, full-field measurements and inverse identification to identify mechanical constitutive parameters, known as Material testing 2.0 (MT2.0) [3].

There are a number of new identification methods, including the Virtual Field Method (VFM) [4], [5] the Equilibrium Gap Method (EGM) [6] and the Constitutive Equation Gap Method (CEGM) [7]. In addition, methods such as the Reciprocity Gap Method (RGM) [8], Finite Element Model Updating (FEMU) [9], [10], Integrated Mechanical Image Correlation (I-MIC) [11], Integrated Digital Image Correlation (I-DIC) [12] and various combined theoretical-experimental approaches [13-14] have been developed. While the novel identification strategies aim to reduce the calibration effort compared to conventional material testing, they mainly utilize the heterogeneity of deformation to determine several material model parameters simultaneously. In the most desirable and ideal case, the concept of MT2.0 would allow the identification of all material parameters from a single heterogeneous test, which remains an open question to date.

Originally, these methods were applied to problems of elasticity, but have evolved over time. Currently, the most commonly used techniques for solving plasticity problems are FEMU and VFM. FEMU, known for its user-friendly approach, uses an optimization algorithm to minimize the discrepancy between the measured and calculated results of finite element analysis (FEA). This method has proven successful in identifying parameters for various plasticity models, including Hill48 [15] and YLD2000-2d [16]. On the other hand, VFM, which is based on the principles of virtual work, was originally developed for the determination of elastic properties [17]. In later applications, however, the range of application was successfully extended to plasticity problems [18-22].

Both FEMU and VFM suffer from certain drawbacks that prevent them from gaining wider acceptance in the industrial engineering community. Inverse identification with FEMU is extremely computationally intensive for plasticity models because the optimization algorithm requires the results of FEM for all material points in a region of interest (ROI) in each iteration, as well as repeating all FEM analyses for perturbed values of all model parameters in each iteration. Each FEM analysis also involves equilibrium iterations, and for each iteration of the FEM analysis, the complete stress field needs to be computed. VFM is much more computationally efficient for elasticity, but becomes computationally intensive for plasticity, where the stress field must also be reconstructed at each iteration.

The second problem hindering wider practical application arises from the conditioning of the optimization problem, which may arise solely from the user's understanding of the underlying identification/optimization problem [23]. Consequently, the optimization process can be ill-posed and lead to poor convergence, a non-unique solution or, in extreme cases, to no solution at all. An illustrative example of such an ill-posed optimization scenario is the attempt to identify isotropic hardening behavior from the linear elastic material response or to determine the hardening exponent for small plastic strains. In both cases, no information is contained in the measured strain
Digital image correlation (DIC) [2] is one of the most commonly used techniques for measuring strain fields. In this approach, sequentially captured images are processed to extract valuable information that facilitates the study of both local behavior and global responses in structures. In addition, DIC can theoretically operate in real time and provides the ability to monitor structural responses under dynamic loading conditions. This feature makes DIC a powerful tool for various applications ranging from material testing to structural analysis and product development.

These novel technologies unleashed a potential for development of novel material calibration strategies, where a richness of full-field measurement of displacements, strains and eventually stresses over a complete Region of interest (ROI) can be analyzed and employed for material calibration. This triggered also a reconsideration of test specimen designs which are needed to be tailored for novel optical measurement technologies to fully exploit the potential, resulting in a new paradigm in material testing using complex tests, full-field measurements and inverse identification to identify mechanical constitutive parameters, known as Material testing 2.0 (MT2.0) [3].

There are a number of new identification methods, including the Virtual Field Method (VFM) [4], [5] the Equilibrium Gap Method (EGM) [6] and the Constitutive Equation Gap Method (CEGM) [7]. In addition, methods such as the Reciprocity Gap Method (RGM) [8], Finite Element Model Updating (FEMU) [9], [10], Integrated Mechanical Image Correlation (I-MIC) [11], Integrated Digital Image Correlation (I-DIC) [12] and various combined theoretical-experimental approaches [13-14] have been developed. While the novel identification strategies aim to reduce the calibration effort compared to conventional material testing, they mainly utilize the heterogeneity of deformation to determine several material model parameters simultaneously. In the most desirable and ideal case, the concept of MT2.0 would allow the identification of all material parameters from a single heterogeneous test, which remains an open question to date.

Originally, these methods were applied to problems of elasticity, but have evolved over time. Currently, the most commonly used techniques for solving plasticity problems are FEMU and VFM. FEMU, known for its user-friendly approach, uses an optimization algorithm to minimize the discrepancy between the measured and calculated results of finite element analysis (FEA). This method has proven successful in identifying parameters for various plasticity models, including Hill48 [15] and YLD2000-2d [16]. On the other hand, VFM, which is based on the principles of virtual work, was originally developed for the determination of elastic properties [17]. In later applications, however, the range of application was successfully extended to plasticity problems [18-22].

Both FEMU and VFM suffer from certain drawbacks that prevent them from gaining wider acceptance in the industrial engineering community. Inverse identification with FEMU is extremely computationally intensive for plasticity models because the optimization algorithm requires the results of FEM for all material points in a region of interest (ROI) in each iteration, as well as repeating all FEM analyses for perturbed values of all model parameters in each iteration. Each FEM analysis also involves equilibrium iterations, and for each iteration of the FEM analysis, the complete stress field needs to be computed. VFM is much more computationally efficient for elasticity, but becomes computationally intensive for plasticity, where the stress field must also be reconstructed at each iteration.

The second problem hindering wider practical application arises from the conditioning of the optimization problem, which may arise solely from the user's understanding of the underlying identification/optimization problem [23]. Consequently, the optimization process can be ill-posed and lead to poor convergence, a non-unique solution or, in extreme cases, to no solution at all. An illustrative example of such an ill-posed optimization scenario is the attempt to identify isotropic hardening behavior from the linear elastic material response or to determine the hardening exponent for small plastic strains. In both cases, no information is contained in the measured strain
fields and there is no physical correlation between the parameters sought and the measured response.

While the two problems described above can be associated with a low sensitivity of the resulting strain or stress field to parameter changes, the more difficult case occurs when parameters are correlated. Such an example would also lead to an ill-posed problem where no unique solution can be found and the optimization process does not effectively converge to the desired solution. An illustrative example is the problem of identifying isotropic strain hardening and plastic anisotropy when all parameters describing hardening and anisotropy are included in the optimization. In this particular case, the same full-field response can be obtained by proportionally increasing the parameters for the hardening curve or the anisotropy.

In this paper, we present an approach to perform an identifiability analysis of the nonlinear Virtual Field Method (VFM) specifically applied to the identification of plastic anisotropy using the perforated biaxial cruciform specimen developed by Coppieters et al. [24]. The first step is to investigate the sensitivity of the stress fields associated with the virtual work and thus determine the rank of the parameter sensitivity. In the final step, we calculate the correlation matrix and find that only a certain set of parameters can be reliably identified for the specific material model and heterogeneous test configuration.

## Methodology

In order to investigate the identifiability of the parameters of a material model when they are identified using nonlinear sensitivity based VFM [25], we have relied on a methodology of virtual experimentation, where a numerical simulation is performed to obtain synthetically deformed images, which are afterwards processed with DIC. The identifiability of the material parameters—when the identification is subjected to nonlinear VFM—is then investigated.

For this particular purpose, the entire methodology was carried out in the following steps, which are shown in Fig. 1:

- Finite element simulation of a perforated biaxial cruciform sheet metal test specimen. The shape of the specimen was taken from [24], [26], while the material of the sheet was assumed to correspond to the actual material behavior of a 1.2 mm thick DC04 steel sheet. The material behavior is assumed to be plastically anisotropic, using the YLD2000-2d yield function to describe plastic anisotropy and Swift's hardening law for isotropic hardening behavior. The material is fully characterized in Coppieters et al. [27].

- Generation of synthetic images of the deformed specimen shape. Based on the calculated displacement field and the corresponding nodal displacements, a deformed speckle images can be generated from a generic optimal speckle pattern. The procedure for generating images is described in detail by Lava et al. [28].

- DIC post-processing. After the deformed images have been generated, the data can be post-processed with a DIC code that generates the “experimental” displacement field. Since the virtual experiment is generated using a known material model, the virtual data can be used to evaluate the identification quality of the nonlinear sensitivity based VFM. It should be noted that the chosen virtual experimentation approach enables to mimic the entire measurement chain, including image noise and effects caused by camera calibration.

- Analysis of identifiability. Prior to the inverse identification procedure using nonlinear VFM, an identifiability analysis should be performed. The identifiability analysis includes a calculation of the internal virtual work sensitivity matrix, which is further analyzed for the sensitivity rank and the collinearity of the parameters. The procedure closely follows the work of Zhang et al. [28] and [29], who applied the identifiability analysis to FEMU and deep notch experiment.
fields and there is no physical correlation between the parameters sought and the measured response.

While the two problems described above can be associated with a low sensitivity of the resulting strain or stress field to parameter changes, the more difficult case occurs when parameters are correlated. Such an example would also lead to an ill-posed problem where no unique solution can be found and the optimization process does not effectively converge to the desired solution. An illustrative example is the problem of identifying isotropic strain hardening and plastic anisotropy when all parameters describing hardening and anisotropy are included in the optimization. In this particular case, the same full-field response can be obtained by proportionally increasing the parameters for the hardening curve or the anisotropy.

In this paper, we present an approach to perform an identifiability analysis of the nonlinear Virtual Field Method (VFM) specifically applied to the identification of plastic anisotropy using the perforated biaxial cruciform specimen developed by Coppieters et al. [24]. The first step is to investigate the sensitivity of the stress fields associated with the virtual work and thus determine the rank of the parameter sensitivity. In the final step, we calculate the correlation matrix and find that only a certain set of parameters can be reliably identified for the specific material model and heterogeneous test configuration.

## Methodology

In order to investigate the identifiability of the parameters of a material model when they are identified using nonlinear sensitivity based VFM [25], we have relied on a methodology of virtual experimentation, where a numerical simulation is performed to obtain synthetically deformed images, which are afterwards processed with DIC. The identifiability of the material parameters—when the identification is subjected to nonlinear VFM—is then investigated.

For this particular purpose, the entire methodology was carried out in the following steps, which are shown in Fig. 1:

- Finite element simulation of a perforated biaxial cruciform sheet metal test specimen. The shape of the specimen was taken from [24], [26], while the material of the sheet was assumed to correspond to the actual material behavior of a 1.2 mm thick DC04 steel sheet. The material behavior is assumed to be plastically anisotropic, using the YLD2000-2d yield function to describe plastic anisotropy and Swift's hardening law for isotropic hardening behavior. The material is fully characterized in Coppieters et al. [27].

- Generation of synthetic images of the deformed specimen shape. Based on the calculated displacement field and the corresponding nodal displacements, a deformed speckle images can be generated from a generic optimal speckle pattern. The procedure for generating images is described in detail by Lava et al. [28].

- DIC post-processing. After the deformed images have been generated, the data can be post-processed with a DIC code that generates the “experimental” displacement field. Since the virtual experiment is generated using a known material model, the virtual data can be used to evaluate the identification quality of the nonlinear sensitivity based VFM. It should be noted that the chosen virtual experimentation approach enables to mimic the entire measurement chain, including image noise and effects caused by camera calibration.

- Analysis of identifiability. Prior to the inverse identification procedure using nonlinear VFM, an identifiability analysis should be performed. The identifiability analysis includes a calculation of the internal virtual work sensitivity matrix, which is further analyzed for the sensitivity rank and the collinearity of the parameters. The procedure closely follows the work of Zhang et al. [28] and [29], who applied the identifiability analysis to FEMU and deep notch experiment.
<figure><img src="image_1685.png" /><figcaption>Fig. 1. The concept of virtual experimentation. The idea is to simulate the material response of the heterogeneous sample with known material parameters (a) and to re-identify the parameters (e) by synthetic image deformation (b) and DIC processing of the deformed images (c). Before inverse identification with the nonlinear VFM, the user should know which parameters are identifiable for the material model used based on the employed heterogeneous test (d).</figcaption></figure>

## Results

The identifiability analysis was conducted following the methodology outlined above. To assess this, the sensitivity of the internal virtual work to the YLD2000-2d parameters was evaluated at each time frame and DIC acquisition point. Utilizing the extracted data, a sensitivity matrix was constructed, and from this matrix, the sensitivity rank of the YLD2000-2d parameters was determined. The results of the sensitivity rank for each parameter are depicted in a bar chart in Fig. 2(a). Notably, parameter $\alpha_7$, associated with shear behavior in YLD2000-2d, significantly influences the internal virtual work. This parameter is conventionally determined through the standard uniaxial tensile test in the diagonal direction to the sheet's rolling direction (45°).

The dominance of shear behavior is evident in Fig. 2(b), where the contour plot displays the maximum principal stress component. It reveals that stress is most pronounced at the fillets, oriented at a 45° angle to the horizontal (rolling) direction, resembling behavior observed in the uniaxial tensile test in the diagonal direction. Finally, from the perspective of the performance of the optimization algorithm, the value of $\alpha_7$ would effectively converge to the reference value, but a slightly different value of $\alpha_7$ would greatly affect the values of other parameters and normally, convergence of such optimization problem is poor. For such behavior a better strategy is to fix or determine $\alpha_7$ from uniaxial tensile test in the diagonal direction. Poor convergence of the optimization problem is usually associated with poor conditioning of the system matrix, which is usually characterized by the conditioning number defined as ratio between the largest and the smallest principal value of the system's matrix. If the conditioning number is high, poor convergence is to be expected.
<figure><img src="image_1684.png" /><figcaption>Fig. 1. The concept of virtual experimentation. The idea is to simulate the material response of the heterogeneous sample with known material parameters (a) and to re-identify the parameters (e) by synthetic image deformation (b) and DIC processing of the deformed images (c). Before inverse identification with the nonlinear VFM, the user should know which parameters are identifiable for the material model used based on the employed heterogeneous test (d).</figcaption></figure>

## Results

The identifiability analysis was conducted following the methodology outlined above. To assess this, the sensitivity of the internal virtual work to the YLD2000-2d parameters was evaluated at each time frame and DIC acquisition point. Utilizing the extracted data, a sensitivity matrix was constructed, and from this matrix, the sensitivity rank of the YLD2000-2d parameters was determined. The results of the sensitivity rank for each parameter are depicted in a bar chart in Fig. 2(a). Notably, parameter $\alpha_7$, associated with shear behavior in YLD2000-2d, significantly influences the internal virtual work. This parameter is conventionally determined through the standard uniaxial tensile test in the diagonal direction to the sheet's rolling direction (45°).

The dominance of shear behavior is evident in Fig. 2(b), where the contour plot displays the maximum principal stress component. It reveals that stress is most pronounced at the fillets, oriented at a 45° angle to the horizontal (rolling) direction, resembling behavior observed in the uniaxial tensile test in the diagonal direction. Finally, from the perspective of the performance of the optimization algorithm, the value of $\alpha_7$ would effectively converge to the reference value, but a slightly different value of $\alpha_7$ would greatly affect the values of other parameters and normally, convergence of such optimization problem is poor. For such behavior a better strategy is to fix or determine $\alpha_7$ from uniaxial tensile test in the diagonal direction. Poor convergence of the optimization problem is usually associated with poor conditioning of the system matrix, which is usually characterized by the conditioning number defined as ratio between the largest and the smallest principal value of the system's matrix. If the conditioning number is high, poor convergence is to be expected.
<figure><img src="image_1686.png" /><figcaption>Fig. 2. Analysis of identifiability by parameter sensitivity rank. (a) When attempting to identify all parameters of YLD2000-2d from a perforated biaxial tensile test, it is clear that the parameter $\alpha_7$ exerts the greatest influence on the internal virtual work. (b) This behavior is expected since the specimen at the fillets, where the maximum principal stress is highest, closely mirrors the behavior observed in the uniaxial tensile test.</figcaption></figure>

Based on the conclusion above, and to fully exploit the potential of heterogeneous testing, it is reasonable to eliminate the parameter $\alpha_7$ from the identification procedure and determine its value from another test. The rationale behind such decision is eliminating the largest eigenvalues from the system's matrix and thus improving the conditioning of the optimization. Based on the assumption that the parameter $\alpha_7$ is fixed, we re-evaluated the sensitivity matrix and parameters' sensitivity rank. The results are presented in Fig. 3(a).

<figure><img src="image_1687.png" /><figcaption>Fig. 3. Analysis of identifiability with parameter $\alpha_7$ fixed. (a) Parameters $\alpha_1$, $\alpha_2$, $\alpha_5$ and $\alpha_8$ have the highest sensitivity rank and are potentially identifiable. (b) This may further be analyzed by collinearity analysis, where the interaction of parameters is analyzed via correlation matrix.</figcaption></figure>
<figure><img src="image_1685.png" /><figcaption>Fig. 2. Analysis of identifiability by parameter sensitivity rank. (a) When attempting to identify all parameters of YLD2000-2d from a perforated biaxial tensile test, it is clear that the parameter $\alpha_7$ exerts the greatest influence on the internal virtual work. (b) This behavior is expected since the specimen at the fillets, where the maximum principal stress is highest, closely mirrors the behavior observed in the uniaxial tensile test.</figcaption></figure>

Based on the conclusion above, and to fully exploit the potential of heterogeneous testing, it is reasonable to eliminate the parameter $\alpha_7$ from the identification procedure and determine its value from another test. The rationale behind such decision is eliminating the largest eigenvalues from the system's matrix and thus improving the conditioning of the optimization. Based on the assumption that the parameter $\alpha_7$ is fixed, we re-evaluated the sensitivity matrix and parameters' sensitivity rank. The results are presented in Fig. 3(a).

<figure><img src="image_1686.png" /><figcaption>Fig. 3. Analysis of identifiability with parameter $\alpha_7$ fixed. (a) Parameters $\alpha_1$, $\alpha_2$, $\alpha_5$ and $\alpha_8$ have the highest sensitivity rank and are potentially identifiable. (b) This may further be analyzed by collinearity analysis, where the interaction of parameters is analyzed via correlation matrix.</figcaption></figure>
In this scenario, it is evident that the parameters $α_1, α_2, α_5$ and $α_8$ have the highest sensitivity rank, suggesting that the test inherently contains sufficient information for their identification. The inclusion of other parameters with a lower sensitivity rank would worsen the conditioning number of the system matrix and impair the optimization. However, it is crucial to consider the possibility that these parameters behave as linearly dependent during the inverse identification process. To investigate the possible interaction between $α_1, α_2, α_5$ and $α_8$, a correlation matrix was computed and is presented in Fig. 3(b).

The correlation matrix reveals pairwise interactions between the parameters, with off-diagonal values exceeding ±0.8 potentially causing ill-conditioning of the identification. In our particular case, focusing on the pairwise interactions between $α_1, α_2, α_5$ and $α_8$, the largest observed value is 0.36, suggesting that these parameters are not significantly correlated, allowing for their unique identification from the test.

It is important to note that the presented analysis relies on synthetic data, and the identification quality of some parameters may be compromised by noisy data. Addressing this concern involves evaluating the parameters' confidence intervals, and this aspect is a subject for future studies [31].

## Summary

To summarize, we have presented the identifiability study of the YLD2000-2d plastic anisotropy model, whose parameters are to be identified from a perforated biaxial tensile test using sensitivity-based nonlinear VFM.

- For this purpose, we performed a virtual experiment in which a numerical simulation is carried out to obtain synthetically deformed images that are subsequently processed with DIC. In a further step, the identifiability of the material parameters is investigated when the identification is subjected to nonlinear VFM.

- We found that when all parameters are subjected to inverse identification, only the parameter $α_7$ can be reliably identified, as the test closely mirrors the behavior observed in the uniaxial tensile test at 45° measured from the rolling direction of the sheet.

- To fully exploit the potential of heterogeneous testing, the parameter $α_7$ was fixed in a further step since a slightly different value of $α_7$ would greatly affect the values of other parameters and normally, convergence of such optimization problem is poor. Based on fixed value of $α_7$ we re-evaluated the identifiability of the other parameters and found that parameters $α_1, α_2, α_5$ and $α_8$ can be uniquely identified from the test, whereas all other parameters should be also fixed. Finally, further studies should be conducted to investigate the quality of the identification of these parameters when applied to noisy data.

## Acknowledgement

M. Halilovič and B. Starman would like to acknowledge the Slovenian Research Agency for its financial support of research program P2-0263. S. Coppieters, B. Starman and Y. Zhang gratefully acknowledge the support from the Research Fund for Coal and Steel under grant agreement No 888153 (EU-RFCS 2019 project No. 888153 | vForm-xSteels).

## Disclaimer

The results reflect only the authors' view, and the European Commission is not responsible for any use that may be made of the information it contains.

## References

[1] F. Pierron, "Material Testing 2.0: A brief review," *Strain*, vol. 59, no. 3, p. e12434, 2023. https://doi.org/10.1111/str.12434

[2] H. Schreier, J.-J. Orteu, and M. A. Sutton, *Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications*. Boston, MA: Springer US, 2009. https://doi.org/10.1007/978-0-387-78747-3
In this scenario, it is evident that the parameters $α_1, α_2, α_5$ and $α_8$ have the highest sensitivity rank, suggesting that the test inherently contains sufficient information for their identification. The inclusion of other parameters with a lower sensitivity rank would worsen the conditioning number of the system matrix and impair the optimization. However, it is crucial to consider the possibility that these parameters behave as linearly dependent during the inverse identification process. To investigate the possible interaction between $α_1, α_2, α_5$ and $α_8$, a correlation matrix was computed and is presented in Fig. 3(b).

The correlation matrix reveals pairwise interactions between the parameters, with off-diagonal values exceeding ±0.8 potentially causing ill-conditioning of the identification. In our particular case, focusing on the pairwise interactions between $α_1, α_2, α_5$ and $α_8$, the largest observed value is 0.36, suggesting that these parameters are not significantly correlated, allowing for their unique identification from the test.

It is important to note that the presented analysis relies on synthetic data, and the identification quality of some parameters may be compromised by noisy data. Addressing this concern involves evaluating the parameters' confidence intervals, and this aspect is a subject for future studies [31].

## Summary

To summarize, we have presented the identifiability study of the YLD2000-2d plastic anisotropy model, whose parameters are to be identified from a perforated biaxial tensile test using sensitivity-based nonlinear VFM.

- For this purpose, we performed a virtual experiment in which a numerical simulation is carried out to obtain synthetically deformed images that are subsequently processed with DIC. In a further step, the identifiability of the material parameters is investigated when the identification is subjected to nonlinear VFM.

- We found that when all parameters are subjected to inverse identification, only the parameter $α_7$ can be reliably identified, as the test closely mirrors the behavior observed in the uniaxial tensile test at 45° measured from the rolling direction of the sheet.

- To fully exploit the potential of heterogeneous testing, the parameter $α_7$ was fixed in a further step since a slightly different value of $α_7$ would greatly affect the values of other parameters and normally, convergence of such optimization problem is poor. Based on fixed value of $α_7$ we re-evaluated the identifiability of the other parameters and found that parameters $α_1, α_2, α_5$ and $α_8$ can be uniquely identified from the test, whereas all other parameters should be also fixed. Finally, further studies should be conducted to investigate the quality of the identification of these parameters when applied to noisy data.

## Acknowledgement

M. Halilovič and B. Starman would like to acknowledge the Slovenian Research Agency for its financial support of research program P2-0263. S. Coppieters, B. Starman and Y. Zhang gratefully acknowledge the support from the Research Fund for Coal and Steel under grant agreement No 888153 (EU-RFCS 2019 project No. 888153 | vForm-xSteels).

## Disclaimer

The results reflect only the authors' view, and the European Commission is not responsible for any use that may be made of the information it contains.

## References

[1] F. Pierron, "Material Testing 2.0: A brief review," *Strain*, vol. 59, no. 3, p. e12434, 2023. https://doi.org/10.1111/str.12434

[2] H. Schreier, J.-J. Orteu, and M. A. Sutton, *Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications*. Boston, MA: Springer US, 2009. https://doi.org/10.1007/978-0-387-78747-3
[3] F. Pierron and M. Grédiac, “Towards Material Testing 2.0. A review of test design for identification of constitutive parameters from full-field measurements,” *Strain*, vol. n/a, no. n/a, p. e12370, 2020. https://doi.org/10.1111/str.12370

[4] M. Grédiac, F. Pierron, S. Avril, E. Toussaint, and M. Rossi, “Virtual Fields Method,” in *Full-Field Measurements and Identification in Solid Mechanics*, John Wiley & Sons, Ltd, 2012, pp. 301–330. doi: 10.1002/9781118578469.ch11

[5] F. Pierron, S. Zhavoronok, and M. Grédiac, “Identification of the through-thickness properties of thick laminated tubes using the virtual fields method,” *Int. J. Solids Struct.*, vol. 37, no. 32, pp. 4437–4453, Aug. 2000. https://doi.org/10.1016/S0020-7683(99)00149-3

[6] D. Claire, F. Hild, and S. Roux, “A finite element formulation to identify damage fields: the equilibrium gap method,” *Int. J. Numer. Methods Eng.*, vol. 61, no. 2, pp. 189–208, 2004. https://doi.org/10.1002/nme.1057

[7] E. Pagnacco, A. Moreau, and D. Lemosse, “Inverse strategies for the identification of elastic and viscoelastic material parameters using full-field measurements,” *Mater. Sci. Eng. A*, vol. 452–453, pp. 737–745, Apr. 2007. https://doi.org/10.1016/j.msea.2006.10.122

[8] H. D. Bui, A. Constantinescu, and H. Maigre, “Numerical identification of linear cracks in 2D elastodynamics using the instantaneous reciprocity gap,” *Inverse Probl.*, vol. 20, no. 4, p. 993, May 2004. https://doi.org/10.1088/0266-5611/20/4/001

[9] F. Mathieu, H. Leclerc, F. Hild, and S. Roux, “Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC,” *Exp. Mech.*, vol. 55, no. 1, pp. 105–119, Jan. 2015. https://doi.org/10.1007/s11340-014-9888-9

[10] J. Kajberg and G. Lindkvist, “Characterisation of materials subjected to large strains by inverse modelling based on in-plane displacement fields,” *Int. J. Solids Struct.*, vol. 41, no. 13, pp. 3439–3459, Jun. 2004. https://doi.org/10.1016/j.ijsolstr.2004.02.021

[11] J. Réthoré, Muhibullah, T. Elguedj, M. Coret, P. Chaudet, and A. Combescure, “Robust identification of elasto-plastic constitutive law parameters from digital images using 3D kinematics,” *Int. J. Solids Struct.*, vol. 50, no. 1, pp. 73–85, Jan. 2013. https://doi.org/10.1016/j.ijsolstr.2012.09.002

[12] A. P. Ruybalid, J. P. M. Hoefnagels, O. van der Sluis, and M. G. D. Geers, “Comparison of the identification performance of conventional FEM updating and integrated DIC,” *Int. J. Numer. Methods Eng.*, vol. 106, no. 4, pp. 298–320, 2016. https://doi.org/10.1002/nme.5127

[13] S. Coppieters and T. Kuwabara, “Identification of Post-Necking Hardening Phenomena in Ductile Sheet Metal,” *Exp. Mech.*, vol. 54, no. 8, pp. 1355–1371, Oct. 2014. https://doi.org/10.1007/s11340-014-9900-4

[14] S. Coppieters, S. Cooreman, H. Sol, P. Van Houtte, and D. Debruyne, “Identification of the post-necking hardening behaviour of sheet metal by comparison of the internal and external work in the necking zone,” *J. Mater. Process. Technol.*, vol. 211, no. 3, pp. 545–552, Mar. 2011. https://doi.org/10.1016/j.jmatprotec.2010.11.015

[15] D. Lecompte, S. Cooreman, S. Coppieters, J. Vantomme, H. Sol, and D. Debruyne, “Parameter identification for anisotropic plasticity model using digital image correlation,” *Eur. J. Comput. Mech.*, vol. 18, no. 3–4, pp. 393–418, Jan. 2009. https://doi.org/10.3166/ejcm.18.393-418

[16] A. Maček, B. Starman, N. Mole, and M. Halilovič, “Calibration of Advanced Yield Criteria Using Uniaxial and Heterogeneous Tensile Test Data,” *Metals*, vol. 10, no. 4, p. 542, Apr. 2020. https://doi.org/10.3390/met10040542

[17] S. Avril and F. Pierron, “General framework for the identification of constitutive parameters from full-field measurements in linear elasticity,” *Int. J. Solids Struct.*, vol. 44, no. 14, pp. 4978–5002, Jul. 2007. https://doi.org/10.1016/j.ijsolstr.2006.12.018
[3] F. Pierron and M. Grédiac, “Towards Material Testing 2.0. A review of test design for identification of constitutive parameters from full-field measurements,” *Strain*, vol. n/a, no. n/a, p. e12370, 2020. https://doi.org/10.1111/str.12370

[4] M. Grédiac, F. Pierron, S. Avril, E. Toussaint, and M. Rossi, “Virtual Fields Method,” in *Full-Field Measurements and Identification in Solid Mechanics*, John Wiley & Sons, Ltd, 2012, pp. 301–330. doi: 10.1002/9781118578469.ch11

[5] F. Pierron, S. Zhavoronok, and M. Grédiac, “Identification of the through-thickness properties of thick laminated tubes using the virtual fields method,” *Int. J. Solids Struct.*, vol. 37, no. 32, pp. 4437–4453, Aug. 2000. https://doi.org/10.1016/S0020-7683(99)00149-3

[6] D. Claire, F. Hild, and S. Roux, “A finite element formulation to identify damage fields: the equilibrium gap method,” *Int. J. Numer. Methods Eng.*, vol. 61, no. 2, pp. 189–208, 2004. https://doi.org/10.1002/nme.1057

[7] E. Pagnacco, A. Moreau, and D. Lemosse, “Inverse strategies for the identification of elastic and viscoelastic material parameters using full-field measurements,” *Mater. Sci. Eng. A*, vol. 452–453, pp. 737–745, Apr. 2007. https://doi.org/10.1016/j.msea.2006.10.122

[8] H. D. Bui, A. Constantinescu, and H. Maigre, “Numerical identification of linear cracks in 2D elastodynamics using the instantaneous reciprocity gap,” *Inverse Probl.*, vol. 20, no. 4, p. 993, May 2004. https://doi.org/10.1088/0266-5611/20/4/001

[9] F. Mathieu, H. Leclerc, F. Hild, and S. Roux, “Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC,” *Exp. Mech.*, vol. 55, no. 1, pp. 105–119, Jan. 2015. https://doi.org/10.1007/s11340-014-9888-9

[10] J. Kajberg and G. Lindkvist, “Characterisation of materials subjected to large strains by inverse modelling based on in-plane displacement fields,” *Int. J. Solids Struct.*, vol. 41, no. 13, pp. 3439–3459, Jun. 2004. https://doi.org/10.1016/j.ijsolstr.2004.02.021

[11] J. Réthoré, Muhibullah, T. Elguedj, M. Coret, P. Chaudet, and A. Combescure, “Robust identification of elasto-plastic constitutive law parameters from digital images using 3D kinematics,” *Int. J. Solids Struct.*, vol. 50, no. 1, pp. 73–85, Jan. 2013. https://doi.org/10.1016/j.ijsolstr.2012.09.002

[12] A. P. Ruybalid, J. P. M. Hoefnagels, O. van der Sluis, and M. G. D. Geers, “Comparison of the identification performance of conventional FEM updating and integrated DIC,” *Int. J. Numer. Methods Eng.*, vol. 106, no. 4, pp. 298–320, 2016. https://doi.org/10.1002/nme.5127

[13] S. Coppieters and T. Kuwabara, “Identification of Post-Necking Hardening Phenomena in Ductile Sheet Metal,” *Exp. Mech.*, vol. 54, no. 8, pp. 1355–1371, Oct. 2014. https://doi.org/10.1007/s11340-014-9900-4

[14] S. Coppieters, S. Cooreman, H. Sol, P. Van Houtte, and D. Debruyne, “Identification of the post-necking hardening behaviour of sheet metal by comparison of the internal and external work in the necking zone,” *J. Mater. Process. Technol.*, vol. 211, no. 3, pp. 545–552, Mar. 2011. https://doi.org/10.1016/j.jmatprotec.2010.11.015

[15] D. Lecompte, S. Cooreman, S. Coppieters, J. Vantomme, H. Sol, and D. Debruyne, “Parameter identification for anisotropic plasticity model using digital image correlation,” *Eur. J. Comput. Mech.*, vol. 18, no. 3–4, pp. 393–418, Jan. 2009. https://doi.org/10.3166/ejcm.18.393-418

[16] A. Maček, B. Starman, N. Mole, and M. Halilovič, “Calibration of Advanced Yield Criteria Using Uniaxial and Heterogeneous Tensile Test Data,” *Metals*, vol. 10, no. 4, p. 542, Apr. 2020. https://doi.org/10.3390/met10040542

[17] S. Avril and F. Pierron, “General framework for the identification of constitutive parameters from full-field measurements in linear elasticity,” *Int. J. Solids Struct.*, vol. 44, no. 14, pp. 4978–5002, Jul. 2007. https://doi.org/10.1016/j.ijsolstr.2006.12.018
[18] M. Grédiac and F. Pierron, "Applying the Virtual Fields Method to the identification of elasto-plastic constitutive parameters," *Int. J. Plast.*, vol. 22, no. 4, pp. 602-627, Apr. 2006. https://doi.org/10.1016/j.ijplas.2005.04.007

[19] Y. Pannier, S. Avril, R. Rotinat, and F. Pierron, "Identification of Elasto-Plastic Constitutive Parameters from Statically Undetermined Tests Using the Virtual Fields Method," *Exp. Mech.*, vol. 46, no. 6, pp. 735-755, Dec. 2006. https://doi.org/10.1007/s11340-006-9822-x

[20] J.-H. Kim, F. Barlat, F. Pierron, and M.-G. Lee, "Determination of Anisotropic Plastic Constitutive Parameters Using the Virtual Fields Method," *Exp. Mech.*, vol. 54, no. 7, pp. 1189-1204, Sep. 2014. https://doi.org/10.1007/s11340-014-9879-x

[21] M. Rossi, F. Pierron, and M. Štamborská, "Application of the virtual fields method to large strain anisotropic plasticity," *Int. J. Solids Struct.*, vol. 97-98, pp. 322-335, Oct. 2016. https://doi.org/10.1016/j.ijsolstr.2016.07.015

[22] M. Rossi and F. Pierron, "Identification of plastic constitutive parameters at large deformations from three dimensional displacement fields," *Comput. Mech.*, vol. 49, no. 1, pp. 53-71, Jan. 2012. https://doi.org/10.1007/s00466-011-0627-0

[23] R. Brun, P. Reichert, and H. R. Künsch, "Practical identifiability analysis of large environmental simulation models," *Water Resour. Res.*, vol. 37, no. 4, pp. 1015-1030, 2001. https://doi.org/10.1029/2000WR900350

[24] S. Coppieters, T. Hakoyama, D. Debruyne, S. Takahashi, and T. Kuwabara, "Inverse Yield Locus Identification of Sheet Metal Using a Complex Cruciform in Biaxial Tension and Digital Image Correlation," *Proceedings*, vol. 2, no. 8, Art. no. 8, 2018. https://doi.org/10.3390/ICEM18-05208

[25] A. Marek, F. M. Davis, and F. Pierron, "Sensitivity-based virtual fields for the non-linear virtual fields method," *Comput. Mech.*, vol. 60, no. 3, pp. 409-431, Sep. 2017. https://doi.org/10.1007/s00466-017-1411-6

[26] S. Cooreman, "Identification of the plastic material behaviour through full-field displacement measurements and inverse methods," VUB, 2008.

[27] S. Coppieters et al., "On the synergy between physical and virtual sheet metal testing: calibration of anisotropic yield functions using a microstructure-based plasticity model," *Int. J. Mater. Form.*, vol. 12, no. 5, pp. 741-759, Sep. 2019. https://doi.org/10.1007/s12289-018-1444-1

[28] P. Lava, S. Cooreman, S. Coppieters, M. De Strycker, and D. Debruyne, "Assessment of measuring errors in DIC using deformation fields generated by plastic FEA," *Opt. Lasers Eng.*, vol. 47, no. 7, pp. 747-753, Jul. 2009. https://doi.org/10.1016/j.optlaseng.2009.03.007

[29] Y. Zhang, A. Yamanaka, S. Cooreman, T. Kuwabara, and S. Coppieters, "Inverse identification of plastic anisotropy through multiple non-conventional mechanical experiments," *Int. J. Solids Struct.*, vol. 285, p. 112534, Dec. 2023. https://doi.org/10.1016/j.ijsolstr.2023.112534

[30] Y. Zhang, A. Van Bael, A. Andrade-Campos, and S. Coppieters, "Parameter identifiability analysis: Mitigating the non-uniqueness issue in the inverse identification of an anisotropic yield function," *Int. J. Solids Struct.*, vol. 243, p. 111543, May 2022. https://doi.org/10.1016/j.ijsolstr.2022.111543

[31] A. Maček, B. Starman, S. Coppieters, J. Urevc, and M. Halilovič, "Confidence intervals of inversely identified material model parameters: A novel two-stage error propagation model based on stereo DIC system uncertainty," *Opt. Lasers Eng.*, vol. 174, p. 107958, Mar. 2024. https://doi.org/10.1016/j.optlaseng.2023.107958
[18] M. Grédiac and F. Pierron, "Applying the Virtual Fields Method to the identification of elasto-plastic constitutive parameters," *Int. J. Plast.*, vol. 22, no. 4, pp. 602-627, Apr. 2006. https://doi.org/10.1016/j.ijplas.2005.04.007

[19] Y. Pannier, S. Avril, R. Rotinat, and F. Pierron, "Identification of Elasto-Plastic Constitutive Parameters from Statically Undetermined Tests Using the Virtual Fields Method," *Exp. Mech.*, vol. 46, no. 6, pp. 735-755, Dec. 2006. https://doi.org/10.1007/s11340-006-9822-x

[20] J.-H. Kim, F. Barlat, F. Pierron, and M.-G. Lee, "Determination of Anisotropic Plastic Constitutive Parameters Using the Virtual Fields Method," *Exp. Mech.*, vol. 54, no. 7, pp. 1189-1204, Sep. 2014. https://doi.org/10.1007/s11340-014-9879-x

[21] M. Rossi, F. Pierron, and M. Štamborská, "Application of the virtual fields method to large strain anisotropic plasticity," *Int. J. Solids Struct.*, vol. 97-98, pp. 322-335, Oct. 2016. https://doi.org/10.1016/j.ijsolstr.2016.07.015

[22] M. Rossi and F. Pierron, "Identification of plastic constitutive parameters at large deformations from three dimensional displacement fields," *Comput. Mech.*, vol. 49, no. 1, pp. 53-71, Jan. 2012. https://doi.org/10.1007/s00466-011-0627-0

[23] R. Brun, P. Reichert, and H. R. Künsch, "Practical identifiability analysis of large environmental simulation models," *Water Resour. Res.*, vol. 37, no. 4, pp. 1015-1030, 2001. https://doi.org/10.1029/2000WR900350

[24] S. Coppieters, T. Hakoyama, D. Debruyne, S. Takahashi, and T. Kuwabara, "Inverse Yield Locus Identification of Sheet Metal Using a Complex Cruciform in Biaxial Tension and Digital Image Correlation," *Proceedings*, vol. 2, no. 8, Art. no. 8, 2018. https://doi.org/10.3390/ICEM18-05208

[25] A. Marek, F. M. Davis, and F. Pierron, "Sensitivity-based virtual fields for the non-linear virtual fields method," *Comput. Mech.*, vol. 60, no. 3, pp. 409-431, Sep. 2017. https://doi.org/10.1007/s00466-017-1411-6

[26] S. Cooreman, "Identification of the plastic material behaviour through full-field displacement measurements and inverse methods," VUB, 2008.

[27] S. Coppieters et al., "On the synergy between physical and virtual sheet metal testing: calibration of anisotropic yield functions using a microstructure-based plasticity model," *Int. J. Mater. Form.*, vol. 12, no. 5, pp. 741-759, Sep. 2019. https://doi.org/10.1007/s12289-018-1444-1

[28] P. Lava, S. Cooreman, S. Coppieters, M. De Strycker, and D. Debruyne, "Assessment of measuring errors in DIC using deformation fields generated by plastic FEA," *Opt. Lasers Eng.*, vol. 47, no. 7, pp. 747-753, Jul. 2009. https://doi.org/10.1016/j.optlaseng.2009.03.007

[29] Y. Zhang, A. Yamanaka, S. Cooreman, T. Kuwabara, and S. Coppieters, "Inverse identification of plastic anisotropy through multiple non-conventional mechanical experiments," *Int. J. Solids Struct.*, vol. 285, p. 112534, Dec. 2023. https://doi.org/10.1016/j.ijsolstr.2023.112534

[30] Y. Zhang, A. Van Bael, A. Andrade-Campos, and S. Coppieters, "Parameter identifiability analysis: Mitigating the non-uniqueness issue in the inverse identification of an anisotropic yield function," *Int. J. Solids Struct.*, vol. 243, p. 111543, May 2022. https://doi.org/10.1016/j.ijsolstr.2022.111543

[31] A. Maček, B. Starman, S. Coppieters, J. Urevc, and M. Halilovič, "Confidence intervals of inversely identified material model parameters: A novel two-stage error propagation model based on stereo DIC system uncertainty," *Opt. Lasers Eng.*, vol. 174, p. 107958, Mar. 2024. https://doi.org/10.1016/j.optlaseng.2023.107958
# Topology optimization for the design of additively manufactured hot stamping tools

YU Daoming¹,², a*, RACHIK Mohamed¹,ᵇ, BRUN Gilles²,ᶜ, BLAISE Alexandre²,ᵈ and SARRE Benjamin²,ᵉ

¹Laboratoire Roberval, FRE UTC-CNRS 2012, Alliance Sorbonne université, Université de Technologie de Compiègne, 60203 Compiègne, France

²ArcelorMittal Global R&D Montataire, 60160 Montataire, France

ᵃdaoming.yu@utc.fr, daoming.yu@arcelormittal.com, ᵇmohamed.rachik@utc.fr,

ᶜgilles.brun@arcelormittal.com, ᵈalexandre.blaise@arcelormittal.com,

ᵉbenjamin.sarre@arcelormittal.com

**Keywords:** Topology Optimization, Additive Manufacturing, Hot-Stamping, Simulation

**Abstract.** The hot stamping process makes it possible to obtain complex-shaped parts with high mechanical properties. The hot stamping production rate is driven by the cooling system inside the tool. Improving the efficiency of the cooling system is a key factor as it reduces part production time and thus lower the cost. But an optimized cooling system requires complex internal geometries. That is why the development of additive manufacturing (AM) opens up new and good prospects for designing and making hot stamping tools with high cooling efficiency. This research proposes a two-step topology optimization procedure for hot stamping tools' design. In the first step, the fluid/thermal topology optimization is used to find the best design for the cooling system. As the printing time and cost for the additively manufactured tool depend on the amount of material used, the second step focuses on removing unnecessary materials in hot stamping tools.

## Introduction

Hot stamping is a widely used process in the manufacturing of automotive parts. A thin sheet of Press Hardened Steel (PHS) is heated in a furnace at a high temperature (around 900 °C) and then formed and quenched in the stamping tool. Quenching is a key step for part productivity. Effective cooling system is necessary to minimize the quenching time [1]. Usually, the cooling channels are drilled straight in the tools. On one hand, it is the easiest way to build the cooling system. On the other hand, this design is not the most efficient one. The aim of this work is to propose a general methodology to achieve the most efficient cooling design while the complexity of the resulting design and the difficulties to manufacture it are overcome by using additive manufacturing.

Modeling hot stamping is a multiphysics problem as it involves the use of mechanics, thermal engineering, and metallurgy. To optimize the cooling channels, solid mechanics, heat transfer and fluid mechanics must be considered. Thus, to reduce the difficulties, the optimization problem can be split into two sub-optimization problems: thermal-mechanical coupled optimization and fluid-thermal coupled optimization.

## Thermal-mechanical coupled topology optimization

Thanks to the advantages of additive manufacturing, the cooling channels can reach closer areas of the tools' surface to increase cooling efficiency and reduce the production time. However, printing the whole tool may highly increase the tool's building cost, as the conventional hot stamping tools are usually cut from solid material. One direct way to reduce the printing cost is to use the hybrid printing method. Only the working part of the tool is printed on a conventional machined substrate. As the printing cost is directly related to the printed volume, the tool's printing cost can thus be reduced. On the contrary, the hybrid printing method can be difficult to apply,

<figure><img src="image_1688.png" /></figure>
# Topology optimization for the design of additively manufactured hot stamping tools

YU Daoming<fnref n="1" />,<fnref n="2" />, a*, RACHIK Mohamed<fnref n="1" />,ᵇ, BRUN Gilles<fnref n="2" />,ᶜ, BLAISE Alexandre<fnref n="2" />,ᵈ and SARRE Benjamin<fnref n="2" />,ᵉ

**Keywords:** Topology Optimization, Additive Manufacturing, Hot-Stamping, Simulation

**Abstract.** The hot stamping process makes it possible to obtain complex-shaped parts with high mechanical properties. The hot stamping production rate is driven by the cooling system inside the tool. Improving the efficiency of the cooling system is a key factor as it reduces part production time and thus lower the cost. But an optimized cooling system requires complex internal geometries. That is why the development of additive manufacturing (AM) opens up new and good prospects for designing and making hot stamping tools with high cooling efficiency. This research proposes a two-step topology optimization procedure for hot stamping tools' design. In the first step, the fluid/thermal topology optimization is used to find the best design for the cooling system. As the printing time and cost for the additively manufactured tool depend on the amount of material used, the second step focuses on removing unnecessary materials in hot stamping tools.

## Introduction

Hot stamping is a widely used process in the manufacturing of automotive parts. A thin sheet of Press Hardened Steel (PHS) is heated in a furnace at a high temperature (around 900 °C) and then formed and quenched in the stamping tool. Quenching is a key step for part productivity. Effective cooling system is necessary to minimize the quenching time [1]. Usually, the cooling channels are drilled straight in the tools. On one hand, it is the easiest way to build the cooling system. On the other hand, this design is not the most efficient one. The aim of this work is to propose a general methodology to achieve the most efficient cooling design while the complexity of the resulting design and the difficulties to manufacture it are overcome by using additive manufacturing.

Modeling hot stamping is a multiphysics problem as it involves the use of mechanics, thermal engineering, and metallurgy. To optimize the cooling channels, solid mechanics, heat transfer and fluid mechanics must be considered. Thus, to reduce the difficulties, the optimization problem can be split into two sub-optimization problems: thermal-mechanical coupled optimization and fluid-thermal coupled optimization.

## Thermal-mechanical coupled topology optimization

Thanks to the advantages of additive manufacturing, the cooling channels can reach closer areas of the tools' surface to increase cooling efficiency and reduce the production time. However, printing the whole tool may highly increase the tool's building cost, as the conventional hot stamping tools are usually cut from solid material. One direct way to reduce the printing cost is to use the hybrid printing method. Only the working part of the tool is printed on a conventional machined substrate. As the printing cost is directly related to the printed volume, the tool's printing cost can thus be reduced. On the contrary, the hybrid printing method can be difficult to apply,

<figure><img src="image_1687.png" /></figure>

<footnote n="1">¹Laboratoire Roberval, FRE UTC-CNRS 2012, Alliance Sorbonne université, Université de Technologie de Compiègne, 60203 Compiègne, France</footnote>

<footnote n="2">²ArcelorMittal Global R&D Montataire, 60160 Montataire, France</footnote>

<footnote>ᵃdaoming.yu@utc.fr, daoming.yu@arcelormittal.com, ᵇmohamed.rachik@utc.fr,</footnote>

<footnote>ᶜgilles.brun@arcelormittal.com, ᵈalexandre.blaise@arcelormittal.com,</footnote>

<footnote>ᵉbenjamin.sarre@arcelormittal.com</footnote>
depending on the tool's geometry, the cooling design, and the used printing technology. For example, the Selective Laser Melting (SLM) technology needs a flat interface between the printed part and the substrate. All upper side of this interface will be printed, irrespective of necessity. Thus, material saving should be taken into account during the tool's design.

Topology optimization is a good tool to answer the printing costs problem. Pure mechanical topology optimization is now very common and widely used in the design of various parts. It can find the best material (solid or void) distribution and use only the necessary solid material to lighten the parts. By reducing the printed hot stamping tool, the printing cost is then reduced as well. Many commercial CAE software already include this function. Topology optimization has already been applied on the cold stamping tools using LS-TaSC with the LS-DYNA solver [2]. The optimized tools have been printed and used successfully, showing good production performance, particularly in terms of fatigue and wear resistance.

In a hot stamping tool, the volume of material and its distribution influence not only the mechanical performance but also the heat transfer from the blank to the cooling channels. Thus, the cooling performance of the tool should be considered in the optimization algorithm. Top3d [3] is used in this study for the thermal-mechanical coupled topology optimization. The idea is to study the influence of the amount of material used on both mechanical and thermal performance of the tool.

A standard topology optimization for linear elastic stress problem consists in minimizing the mechanical compliance of the system which is written as:

$$min: J_m^M = F^T U = U^T K_m U$$

$$s.t.: \frac{v_s}{v^*} \le f \qquad (1)$$

where $J_m^M$ is the mechanical compliance of the system, $U$ and $F$ are respectively the nodal displacement and the nodal force vectors, $K_m$ is the stiffness matrix, $v_s$ and $v^*$ are respectively the used material volume and design domain volume. $f$ is the volume fraction.

Heat conduction optimization problem is very similar to the mechanical optimization problem. Instead of minimizing the system's mechanical compliance, the thermal compliance is minimized, by replacing the nodal displacement vector and the global stiffness matrix by the nodal temperature vector and the conductivity matrix respectively:

$$min: J_t^M = T^T K_t T$$

$$s.t.: \frac{v_s}{v^*} \le f \qquad (2)$$

where $J_t^M$ is the so-called system's thermal compliance, $T$ is the nodal temperature vector, $K_t$ is the conductivity matrix.

To combine the mechanical and the thermal optimization, the objective function needs to consider the tool's mechanical and thermal performance at the same time. As the values of $J_m^M$ and $J_t^M$ are not of the same order of magnitude, a weighting coefficient is usually used to combine the two objective functions. Choosing the appropriate value of such a weighting coefficient can be a tricky task. To overcome this difficulty, the coupled optimization problem can be formulated as:

$$min: J = J_m^M * J_t^M$$

$$s.t.: \frac{v_s}{v^*} \le f \qquad (3)$$
depending on the tool's geometry, the cooling design, and the used printing technology. For example, the Selective Laser Melting (SLM) technology needs a flat interface between the printed part and the substrate. All upper side of this interface will be printed, irrespective of necessity. Thus, material saving should be taken into account during the tool's design.

Topology optimization is a good tool to answer the printing costs problem. Pure mechanical topology optimization is now very common and widely used in the design of various parts. It can find the best material (solid or void) distribution and use only the necessary solid material to lighten the parts. By reducing the printed hot stamping tool, the printing cost is then reduced as well. Many commercial CAE software already include this function. Topology optimization has already been applied on the cold stamping tools using LS-TaSC with the LS-DYNA solver [2]. The optimized tools have been printed and used successfully, showing good production performance, particularly in terms of fatigue and wear resistance.

In a hot stamping tool, the volume of material and its distribution influence not only the mechanical performance but also the heat transfer from the blank to the cooling channels. Thus, the cooling performance of the tool should be considered in the optimization algorithm. Top3d [3] is used in this study for the thermal-mechanical coupled topology optimization. The idea is to study the influence of the amount of material used on both mechanical and thermal performance of the tool.

A standard topology optimization for linear elastic stress problem consists in minimizing the mechanical compliance of the system which is written as:

$$min: J_m^M = F^T U = U^T K_m U$$

$$s.t.: \frac{v_s}{v^*} \le f \qquad (1)$$

where $J_m^M$ is the mechanical compliance of the system, $U$ and $F$ are respectively the nodal displacement and the nodal force vectors, $K_m$ is the stiffness matrix, $v_s$ and $v^*$ are respectively the used material volume and design domain volume. $f$ is the volume fraction.

Heat conduction optimization problem is very similar to the mechanical optimization problem. Instead of minimizing the system's mechanical compliance, the thermal compliance is minimized, by replacing the nodal displacement vector and the global stiffness matrix by the nodal temperature vector and the conductivity matrix respectively:

$$min: J_t^M = T^T K_t T$$

$$s.t.: \frac{v_s}{v^*} \le f \qquad (2)$$

where $J_t^M$ is the so-called system's thermal compliance, $T$ is the nodal temperature vector, $K_t$ is the conductivity matrix.

To combine the mechanical and the thermal optimization, the objective function needs to consider the tool's mechanical and thermal performance at the same time. As the values of $J_m^M$ and $J_t^M$ are not of the same order of magnitude, a weighting coefficient is usually used to combine the two objective functions. Choosing the appropriate value of such a weighting coefficient can be a tricky task. To overcome this difficulty, the coupled optimization problem can be formulated as:

$$min: J = J_m^M * J_t^M$$

$$s.t.: \frac{v_s}{v^*} \le f \qquad (3)$$
where $J^M$ is the coupled objective function which is the product of $J_m^M$ and $J_t^M$. This can be interpreted as: the algorithm seeks the distribution of the material to ensure good mechanical and thermal performance at the same time.

## Thermal-mechanical optimization of a simplified tool

Common hot-stamped automotive parts generally have an omega-shaped cross-section, such as the A or B pillars, as the omega shape offers a compromise between part's mechanical strength and formability. The previously described optimization procedure is then used to investigate a simplified Omega tool. The simplified Omega tool geometry, the mechanical boundary conditions and the thermal boundary conditions are described in Fig. 1. On the tool working surface, the thermomechanical loading experienced by the tool during the forming process is represented, in a simplified manner, by a pressure and a heat flux density.

<figure><img src="image_1689.png" /><figcaption>Fig. 1. Left) Simplified Omega punch used for optimization; Right) Boundary conditions for optimization.</figcaption></figure>

The displacement is set to zero on the fixed face of the tool and the temperature is set to zero on the same face. The tool has 6 cooling channels 10 mm in diameter and for simplicity, the temperature is set to zero at the surface of the channels. The optimization algorithm focuses on modifying only the solid material within the design domain (grey zone in Fig. 1 (Right)), while the position and geometry of the cooling channels remain unchanged because fluid mechanics is not taken into account.

To demonstrate the capability of the optimization procedure, five optimal configurations were determined by varying the solid volume fraction parameter $f$ in Eq. 3, modifying the amount of solid material in the tool. The five optimal configurations are shown in Fig. 2. The vertical and horizontal structures observed in the tool result from the optimization algorithm. This ensures the mechanical strength of the tool. In addition, the material between the working surface of the tool and the cooling channels is kept to ensure heat transfer between the latter. These optimal designs are compared with the original tool in terms of thermomechanical behavior using finite element models with Abaqus software.
where $J^M$ is the coupled objective function which is the product of $J_m^M$ and $J_t^M$. This can be interpreted as: the algorithm seeks the distribution of the material to ensure good mechanical and thermal performance at the same time.

## Thermal-mechanical optimization of a simplified tool

Common hot-stamped automotive parts generally have an omega-shaped cross-section, such as the A or B pillars, as the omega shape offers a compromise between part's mechanical strength and formability. The previously described optimization procedure is then used to investigate a simplified Omega tool. The simplified Omega tool geometry, the mechanical boundary conditions and the thermal boundary conditions are described in Fig. 1. On the tool working surface, the thermomechanical loading experienced by the tool during the forming process is represented, in a simplified manner, by a pressure and a heat flux density.

<figure><img src="image_1688.png" /><figcaption>Fig. 1. Left) Simplified Omega punch used for optimization; Right) Boundary conditions for optimization.</figcaption></figure>

The displacement is set to zero on the fixed face of the tool and the temperature is set to zero on the same face. The tool has 6 cooling channels 10 mm in diameter and for simplicity, the temperature is set to zero at the surface of the channels. The optimization algorithm focuses on modifying only the solid material within the design domain (grey zone in Fig. 1 (Right)), while the position and geometry of the cooling channels remain unchanged because fluid mechanics is not taken into account.

To demonstrate the capability of the optimization procedure, five optimal configurations were determined by varying the solid volume fraction parameter $f$ in Eq. 3, modifying the amount of solid material in the tool. The five optimal configurations are shown in Fig. 2. The vertical and horizontal structures observed in the tool result from the optimization algorithm. This ensures the mechanical strength of the tool. In addition, the material between the working surface of the tool and the cooling channels is kept to ensure heat transfer between the latter. These optimal designs are compared with the original tool in terms of thermomechanical behavior using finite element models with Abaqus software.
<figure><img src="image_1690.png" /><figcaption>Fig. 1 Different optimal designs obtained by thermal-mechanical coupled topology optimization.</figcaption></figure>

To validate the tool cooling performances, the optimized tool shapes are subjected to a cyclic thermal loading and the heat transfer is analyzed for each configuration. For each thermal cycle, a blank with an initial temperature of 750 °C (to take into account the transfer of blank from the furnace to the press) is brought into contact with the tool for 8 seconds. This corresponds to the quenching operation during the hot stamping process. Afterwards, the contact condition is deactivated for 22 s (zero heat transfer to represent the tool opening operation). The material properties, summarized in Table 1, are supposed as temperature independent. A total of 30 cycles were simulated to ensure stabilization of the temperature field.

<figure><img src="image_1691.png" /><figcaption>Fig. 2 Blank temperature after the 30th quenching.</figcaption></figure>

In hot stamping process, the quenching time is determined by the temperature of the hot spot on the work piece when the tool is opened. Thus, the temperature of the blank after the 30th strokes is shown in the Fig. 3. It can be observed that the hot spots are located at point C and D due to the geometry of the tool.

The difference between the optimized tool's hot spot temperature and that of the initial tool is depicted in Fig. 4 for each of the five investigated configurations. These results show that the hot spot temperature increases when the quantity of solid material is reduced. But it should be noted that when 60 % of the material is used in the tool, the temperature only increases by 2 °C. For a hot spot temperature of around 160 °C, it is obvious that this slight difference of 2 °C will not impact productivity.
<figure><img src="image_1689.png" /><figcaption>Fig. 1 Different optimal designs obtained by thermal-mechanical coupled topology optimization.</figcaption></figure>

To validate the tool cooling performances, the optimized tool shapes are subjected to a cyclic thermal loading and the heat transfer is analyzed for each configuration. For each thermal cycle, a blank with an initial temperature of 750 °C (to take into account the transfer of blank from the furnace to the press) is brought into contact with the tool for 8 seconds. This corresponds to the quenching operation during the hot stamping process. Afterwards, the contact condition is deactivated for 22 s (zero heat transfer to represent the tool opening operation). The material properties, summarized in Table 1, are supposed as temperature independent. A total of 30 cycles were simulated to ensure stabilization of the temperature field.

<figure><img src="image_1690.png" /><figcaption>Fig. 2 Blank temperature after the 30th quenching.</figcaption></figure>

In hot stamping process, the quenching time is determined by the temperature of the hot spot on the work piece when the tool is opened. Thus, the temperature of the blank after the 30th strokes is shown in the Fig. 3. It can be observed that the hot spots are located at point C and D due to the geometry of the tool.

The difference between the optimized tool's hot spot temperature and that of the initial tool is depicted in Fig. 4 for each of the five investigated configurations. These results show that the hot spot temperature increases when the quantity of solid material is reduced. But it should be noted that when 60 % of the material is used in the tool, the temperature only increases by 2 °C. For a hot spot temperature of around 160 °C, it is obvious that this slight difference of 2 °C will not impact productivity.
<table><caption>Table 1. Simulation parameters.</caption><thead><tr><th></th><th>Tool material properties (H13)</th><th>Blank material properties (Usibor®1500)</th></tr></thead><tbody><tr><td>Thermal conductivity (W/m²·K)</td><td>25</td><td>20</td></tr><tr><td>Density (kg/m³)</td><td colspan="2">7800</td></tr><tr><td>Specific heat capacity (J/kg·K)</td><td colspan="2">460</td></tr><tr><td rowspan="2">Contact property (thermal transfer coefficient) (W/m²·K)</td><td colspan="2">On horizontal surfaces: 4730</td></tr><tr><td colspan="2">On vertical surfaces: 4140</td></tr><tr><td>Film condition on bottom surface and cooling channels' surfaces</td><td colspan="2">3340W/m²·K to 25°C</td></tr><tr><td>Thermal expansion coefficient (J/K)</td><td colspan="2">1.26 10<sup>-5</sup></td></tr></tbody></table>

<figure><img src="image_1692.png" /><figcaption>Fig. 3 Hot spot's temperature difference.</figcaption></figure>

In addition to cooling efficiency, mechanical resistance is the other requirement for a hot stamping tool. To verify that the removal of solid materials from the tool does not excessively compromise its load capacity, a stress analysis was carried out on the optimized designs of the tool. For these stress analyses, the boundary conditions are those used for topology optimization. Regarding loading, two load cases were considered. The first load case consists of a static pressure applied on the tool's working surface. The obtained von Mises stress field and the displacement field are depicted in Fig. 5 for both the original tool and the optimized tool with 50 % material. Despite the maximum stress increasing from 119.1 MPa to 177 MPa, it remains significantly below the yielding stress of H13 steel for example (a typical steel used for hot stamping tools). The maximum displacement increases from 0.01264 mm to 0.03407 mm which is again very small.

<figure><img src="image_1693.png" /><figcaption>Fig. 4 von Mises stress field and the displacement field under static loads.</figcaption></figure>

The second load case combines the previously described one (pressure on the tool's working surface) with the cyclic thermal loading (thermal expansion). During the last quenching cycle, the maximum stress is reached after a contact time between the flank and the tool of 0.8 seconds. The obtained von Mises stress fields are shown in Fig. 6, and the stress evolution along the surface is plotted in Fig. 7. The maximum stress is located at points C and D due to the local high temperature. It should be noted that thermal stress is strongly linked to the temperature field. Since the tool are subjected to the same thermal loads and the cooling performance has been well optimized, the stress distributions on the contact surface are very similar. Compared to the original tool, the optimized tool with 50 % material has approximately the same maximum stress. Thus, the optimized tool has the same load bearing capacity as the initial tool with a material saving of 50 %.
<table><caption>Table 1. Simulation parameters.</caption><thead><tr><th></th><th>Tool material properties (H13)</th><th>Blank material properties (Usibor®1500)</th></tr></thead><tbody><tr><td>Thermal conductivity (W/m²·K)</td><td>25</td><td>20</td></tr><tr><td>Density (kg/m³)</td><td colspan="2">7800</td></tr><tr><td>Specific heat capacity (J/kg·K)</td><td colspan="2">460</td></tr><tr><td>Contact property (thermal transfer coefficient) (W/m²·K)</td><td colspan="2">On horizontal surfaces: 4730<br>On vertical surfaces: 4140</td></tr><tr><td>Film condition on bottom surface and cooling channels' surfaces</td><td colspan="2">3340W/m²·K to 25°C</td></tr><tr><td>Thermal expansion coefficient (J/K)</td><td colspan="2">1.26 10⁻⁵</td></tr></tbody></table>

<figure><img src="image_1691.png" /><figcaption>Fig. 3 Hot spot's temperature difference.</figcaption></figure>

In addition to cooling efficiency, mechanical resistance is the other requirement for a hot stamping tool. To verify that the removal of solid materials from the tool does not excessively compromise its load capacity, a stress analysis was carried out on the optimized designs of the tool. For these stress analyses, the boundary conditions are those used for topology optimization. Regarding loading, two load cases were considered. The first load case consists of a static pressure applied on the tool's working surface. The obtained von Mises stress field and the displacement field are depicted in Fig. 5 for both the original tool and the optimized tool with 50 % material. Despite the maximum stress increasing from 119.1 MPa to 177 MPa, it remains significantly below the yielding stress of H13 steel for example (a typical steel used for hot stamping tools). The maximum displacement increases from 0.01264 mm to 0.03407 mm which is again very small.

<figure><img src="image_1692.png" /><figcaption>Fig. 4 von Mises stress field and the displacement field under static loads.</figcaption></figure>

The second load case combines the previously described one (pressure on the tool's working surface) with the cyclic thermal loading (thermal expansion). During the last quenching cycle, the maximum stress is reached after a contact time between the flank and the tool of 0.8 seconds. The obtained von Mises stress fields are shown in Fig. 6, and the stress evolution along the surface is plotted in Fig. 7. The maximum stress is located at points C and D due to the local high temperature. It should be noted that thermal stress is strongly linked to the temperature field. Since the tool are subjected to the same thermal loads and the cooling performance has been well optimized, the stress distributions on the contact surface are very similar. Compared to the original tool, the optimized tool with 50 % material has approximately the same maximum stress. Thus, the optimized tool has the same load bearing capacity as the initial tool with a material saving of 50 %.
<figure><img src="image_1694.png" /><figcaption>Fig. 5 von Mises stress field, under thermal expansion Left) Original tool; Right) Tool with 50 % material.</figcaption></figure>

<figure><img src="image_1695.png" /><figcaption>Fig. 6. Surface von Mises stress at t=0.8 s.</figcaption></figure>

## Fluid-thermal coupled topology optimization

As presented in the previous sections, optimizing the material used in the tool may reduce the printing cost without changing the tool's cooling performance and its load bearing capacity. However, to improve the cooling efficiency of the tool, the geometry and position of the cooling channels need to be optimized and this involves taking into account the fluid flow in the cooling channels. To achieve this, a fluid-thermal coupled topology optimization algorithm is developed in the software FreeFEM [4]. Most topological optimization algorithms in fluid mechanics focus on laminar flow with the Navier-Stokes (NS) equations. [5]. When the flow velocity becomes high, the turbulence appears in the flow. The turbulence problem in fluid mechanics is highly nonlinear and unstable, which makes fluid topology optimization a very difficult task. Research reported in the literature shows that fluid topology optimization algorithms have been developed and successfully used for flows with Reynolds numbers up to 200. In the hot stamping tools used for productions, the fluid flow in the cooling channels generally has a Reynolds number of between 4000 and 11000 [6]. Such high flow velocities are sought to increase the cooling efficiency of the tool, as turbulence promotes heat transfer from the tool to the cooling fluid. To overcome this difficulty, the NS equation is simplified as Stokes equation. Without the convective term in the NS equation, the turbulence is neglected, and the flow is considered as fully laminar. This simplification may lead to errors in the calculation of the temperature field in the fluid domain. However, the temperature is more sensitive to the position of the channels. Thus, this method should be able to create continuous channels at hot areas to dissipate heat.

In optimizing heat transfer systems, both hydraulic performance and thermal performance are concurrently considered. The pressure drop of the system is usually minimized to prevent unreliable solutions. It can be written as [7]:

$$J_f^F = - \int_{\Gamma_{Inlet} \cup \Gamma_{Outlet}} p \vec{u} \cdot \vec{n} \, d\Gamma \quad (4)$$

where **p** is the flow pressure and $\vec{u}$ is the velocity vector.
<figure><img src="image_1693.png" /><figcaption>Fig. 5 von Mises stress field, under thermal expansion Left) Original tool; Right) Tool with 50 % material.</figcaption></figure>

<figure><img src="image_1694.png" /><figcaption>Fig. 6. Surface von Mises stress at t=0.8 s.</figcaption></figure>

Fluid-thermal coupled topology optimization

As presented in the previous sections, optimizing the material used in the tool may reduce the printing cost without changing the tool's cooling performance and its load bearing capacity. However, to improve the cooling efficiency of the tool, the geometry and position of the cooling channels need to be optimized and this involves taking into account the fluid flow in the cooling channels. To achieve this, a fluid-thermal coupled topology optimization algorithm is developed in the software FreeFEM [4]. Most topological optimization algorithms in fluid mechanics focus on laminar flow with the Navier-Stokes (NS) equations. [5]. When the flow velocity becomes high, the turbulence appears in the flow. The turbulence problem in fluid mechanics is highly nonlinear and unstable, which makes fluid topology optimization a very difficult task. Research reported in the literature shows that fluid topology optimization algorithms have been developed and successfully used for flows with Reynolds numbers up to 200. In the hot stamping tools used for productions, the fluid flow in the cooling channels generally has a Reynolds number of between 4000 and 11000 [6]. Such high flow velocities are sought to increase the cooling efficiency of the tool, as turbulence promotes heat transfer from the tool to the cooling fluid. To overcome this difficulty, the NS equation is simplified as Stokes equation. Without the convective term in the NS equation, the turbulence is neglected, and the flow is considered as fully laminar. This simplification may lead to errors in the calculation of the temperature field in the fluid domain. However, the temperature is more sensitive to the position of the channels. Thus, this method should be able to create continuous channels at hot areas to dissipate heat.

In optimizing heat transfer systems, both hydraulic performance and thermal performance are concurrently considered. The pressure drop of the system is usually minimized to prevent unreliable solutions. It can be written as [7]:

$$J_f^F = - \int_{\Gamma_{Inlet} \cup \Gamma_{Outlet}} p \vec{u} \cdot \vec{n} \, d\Gamma \quad (4)$$

where **p** is the flow pressure and $\vec{u}$ is the velocity vector.
As explained earlier, the cooling performance of a hot stamping tool is quantified by the workpiece hot spot temperature. By eliminating the effects of poor contact between the tool and the workpiece, the hot spot on the tool and the hot spot on the workpiece are identical. Consequently, the maximum temperature on the working surface of the tool will be used. When the optimization goal is to minimize the hot spot temperature on the tool working surface, the objective function can be written as [8]:

$$J_t^F = \left( \int_{\Gamma_{Contact}} T^\beta d\Gamma \right)^{\frac{1}{\beta}} \qquad (5)$$

where T is the temperature field, β is a positive large number and ΓContact is the tool working surface.

However, the calculation of extremely large values of Tβ may lead to computational issues. To overcome this difficulty, the following thermal objective function is used.

$$J_t^F = \int_{\Gamma_{Contact}} \left(\frac{T}{40}\right)^{20} d\Gamma \qquad (6)$$

The values are chosen based on numerical tests. As the temperature on the surface is always higher than 40 °C, the term of T/40 is larger than 1. By powered by 20, the cooling performance at the hot spot takes more importance during the optimization. Combining Eq. 4 and Eq.6, the coupled fluid-thermal topology optimization problem can be reformulated as:

$$\begin{array}{l} \min: \quad J^F = J_f^F + \tau_{obj} * J_t^F \\ \text{s.t.}: \quad \frac{V_f}{V^*} \le f \end{array} \qquad (7)$$

where τobj is a weighting factor for optimization function and Vf is the fluid volume in the domain.

Applying the topology optimization to a large-scale 3D problem such as the hot stamping tool involves many computation difficulties, in particular the efficiency in terms of computational time. To achieve this, the domain decomposition method is used for parallel computation [9]. Highly efficient solvers are also implemented for solving fluid equations, thanks to the interface between FreeFEM and PETSc [13]. Mesh adaptivity is also used to ensure the solution quality while keeping reasonable computational time. These mentioned methods are implemented into optimization algorithm to ensure that the computation time remains acceptable. To address the challenge posed by heat transport problem with high Péclet numbers, the Streamline Upwind Petrov-Galerkin (SUPG) type stabilization method is employed along with adaptive mesh near the fluid solid interface [10]. The techniques described above are combined with the adjoint method to achieve a highly efficient optimization procedure.

## Fluid-thermal topology optimization of a simplified tool

The developed fluid-thermal topology optimization procedure is tested on a simplified geometry, Fig. 8A, which is representative of a hot stamping tool. Water is used as the coolant for cooling channels. The value of volume fraction parameter f in the Eq. 7 is chosen such that the fluid volume is identical to the one used in a conventional cooling channel shown in the Fig. 9A.
As explained earlier, the cooling performance of a hot stamping tool is quantified by the workpiece hot spot temperature. By eliminating the effects of poor contact between the tool and the workpiece, the hot spot on the tool and the hot spot on the workpiece are identical. Consequently, the maximum temperature on the working surface of the tool will be used. When the optimization goal is to minimize the hot spot temperature on the tool working surface, the objective function can be written as [8]:

$$J_t^F = \left( \int_{\Gamma_{Contact}} T^\beta d\Gamma \right)^{\frac{1}{\beta}} \qquad (5)$$

where T is the temperature field, β is a positive large number and Γcontact is the tool working surface.

However, the calculation of extremely large values of Tβ may lead to computational issues. To overcome this difficulty, the following thermal objective function is used.

$$J_t^F = \int_{\Gamma_{Contact}} \left(\frac{T}{40}\right)^{20} d\Gamma \qquad (6)$$

The values are chosen based on numerical tests. As the temperature on the surface is always higher than 40 °C, the term of T/40 is larger than 1. By powered by 20, the cooling performance at the hot spot takes more importance during the optimization. Combining Eq. 4 and Eq.6, the coupled fluid-thermal topology optimization problem can be reformulated as:

$$\begin{array}{l} \min: \quad J^F = J_f^F + \tau_{obj} * J_t^F \\ \text{s.t.}: \quad \frac{V_f}{V^*} \le f \end{array} \qquad (7)$$

where τobj is a weighting factor for optimization function and Vf is the fluid volume in the domain.

Applying the topology optimization to a large-scale 3D problem such as the hot stamping tool involves many computation difficulties, in particular the efficiency in terms of computational time. To achieve this, the domain decomposition method is used for parallel computation [9]. Highly efficient solvers are also implemented for solving fluid equations, thanks to the interface between FreeFEM and PETSc [13]. Mesh adaptivity is also used to ensure the solution quality while keeping reasonable computational time. These mentioned methods are implemented into optimization algorithm to ensure that the computation time remains acceptable. To address the challenge posed by heat transport problem with high Péclet numbers, the Streamline Upwind Petrov-Galerkin (SUPG) type stabilization method is employed along with adaptive mesh near the fluid solid interface [10]. The techniques described above are combined with the adjoint method to achieve a highly efficient optimization procedure.

## Fluid-thermal topology optimization of a simplified tool

The developed fluid-thermal topology optimization procedure is tested on a simplified geometry, Fig. 8A, which is representative of a hot stamping tool. Water is used as the coolant for cooling channels. The value of volume fraction parameter f in the Eq. 7 is chosen such that the fluid volume is identical to the one used in a conventional cooling channel shown in the Fig. 9A.
<figure><img src="image_1696.png" /><figcaption>Fig. 7. A) Simplified geometry for optimization B) Design domain for optimization.</figcaption></figure>

As the optimization algorithm is for steady laminar flow, the cyclic thermal load is simplified as a constant heat flux 53625 W/m², applied on the contact surface. This heat flux represents the averaged heat flux resulting from a 0.75 mm thickness blank being cooled from 750 °C to 200 °C within a cycle time of 30 seconds. To prevent water from going directly to the surface, creating an unrealizable solution, a solid offset on surfaces is set as a constraint for optimization (Fig. 8B). The inlet temperature is fixed as 20 °C. The flow velocity at inlet is supposed as a fully developed laminar flow using Hagen-Poiseuille equation:

$$ \vec{u}_{\text{Inlet}} \cdot \vec{n} = u_{\text{max}} \left(1 - \frac{r^2}{R^2}\right) \qquad (8) $$

where $u_{max}$ is the maximum velocity at the center of inlet, chosen as 1 m/s. r is the distance to the inlet center and R is the radius of the inlet which is 4 mm.

<figure><img src="image_1697.png" /><figcaption>Fig. 8 A) Conventional cooling channel B) Rectangular cooling channel.</figcaption></figure>

Choosing different values for the weight parameter $\tau_{obj}$ in Eq.7 may lead to different design solutions. When $\tau_{obj}$ is too high, the algorithm would neglect the hydraulic optimization leading to unrealizable solution with non-connected channels. On the other hand, when $\tau_{obj}$ is too low, the cooling performance of the tool would be degraded. In this work, an interesting solution is found by initializing $\tau_{obj}$ at 1000 and gradually decreasing it to 10 during the iterations. A twisted shape channel is obtained to better mix the water inside the channel. The solution is then exported and reconstructed in CAD software, shown in the Fig. 10.

Obviously, the optimization algorithm directs the water flow toward the contact surface to extract the heat. Moreover, as the $\tau_{obj}$ is high enough, a twisted channel geometry is found. As the cold water is heated up while flowing through the channel, the topology optimization algorithm has well captured the temperature gradient in the channel. The twisted geometry promotes the mixing of heated water with cold water to increase cooling efficiency. Some similar geometries are artificially created following the same idea [11].
<figure><img src="image_1695.png" /><figcaption>Fig. 7. A) Simplified geometry for optimization B) Design domain for optimization.</figcaption></figure>

As the optimization algorithm is for steady laminar flow, the cyclic thermal load is simplified as a constant heat flux 53625 W/m², applied on the contact surface. This heat flux represents the averaged heat flux resulting from a 0.75 mm thickness blank being cooled from 750 °C to 200 °C within a cycle time of 30 seconds. To prevent water from going directly to the surface, creating an unrealizable solution, a solid offset on surfaces is set as a constraint for optimization (Fig. 8B). The inlet temperature is fixed as 20 °C. The flow velocity at inlet is supposed as a fully developed laminar flow using Hagen-Poiseuille equation:

$$ \vec{u}_{\text{Inlet}} \cdot \vec{n} = u_{\text{max}} \cdot \left(1 - \frac{r^2}{R^2}\right) \qquad (8) $$

where $u_{max}$ is the maximum velocity at the center of inlet, chosen as 1 m/s. r is the distance to the inlet center and R is the radius of the inlet which is 4 mm.

<figure><img src="image_1696.png" /><figcaption>Fig. 8 A) Conventional cooling channel B) Rectangular cooling channel.</figcaption></figure>

Choosing different values for the weight parameter $\tau_{obj}$ in Eq.7 may lead to different design solutions. When $\tau_{obj}$ is too high, the algorithm would neglect the hydraulic optimization leading to unrealizable solution with non-connected channels. On the other hand, when $\tau_{obj}$ is too low, the cooling performance of the tool would be degraded. In this work, an interesting solution is found by initializing $\tau_{obj}$ at 1000 and gradually decreasing it to 10 during the iterations. A twisted shape channel is obtained to better mix the water inside the channel. The solution is then exported and reconstructed in CAD software, shown in the Fig. 10.

Obviously, the optimization algorithm directs the water flow toward the contact surface to extract the heat. Moreover, as the $\tau_{obj}$ is high enough, a twisted channel geometry is found. As the cold water is heated up while flowing through the channel, the topology optimization algorithm has well captured the temperature gradient in the channel. The twisted geometry promotes the mixing of heated water with cold water to increase cooling efficiency. Some similar geometries are artificially created following the same idea [11].
The optimized cooling channel has been compared with the conventional cooling channel (Fig. 9A) though a CFD model developed in COMSOL Multiphysics. As the optimized channel has a larger width compared to the conventional channel, a rectangular cooling channel (Fig. 9B) is included in the comparison. All these three cooling channels are designed with the same fluid volume. The steady $k-\epsilon$ turbulent model is chosen to simulate the real situation [12]. The domain is meshed using tetrahedral elements with 8 layers of prismatic element in the fluid domain. After a mesh convergence study, the global element size in the fluid is set to 1 mm.

<figure><img src="image_1698.png" /><figcaption>Fig. 9 Optimal solution obtained by fluid-thermal coupled topology optimization.</figcaption></figure>

Fig. 11 shows a comparison between three designs in terms of temperature evolution along the contact surface centerline. It can be observed that the conventional channel (circular shape) has an overall higher temperature along the surface centerline. The temperature with rectangular channel increases almost linearly. The highest temperatures are all located near the end of the surface due to the geometry. Though, the temperature with optimized channel is better controlled. The maximum temperature obtained with the optimized channel is more than 5 °C lower than those obtained with the other designs. Considering this is a simplified case, the gain in temperature may be much higher while facing to a more complex real tool.

<figure><img src="image_1699.png" /><figcaption>Fig. 10 Temperature comparison with three cooling channels.</figcaption></figure>

## Summary

In this study, topology optimization has shown its high suitability for designing hot stamping tools. In thermal-mechanical coupled optimization, around 40% of the solid material in a simplified Omega tool, without changing the tool's cooling performance. Even with a 50% material reduction in material, the reduction in tool load bearing capacity is not critical. This may lead to the price reduction and the building time saving for hot stamping tools. In fluid-thermal coupled optimization, the algorithm shows well its ability for cooling system designs. It does not only send water to the contact surface, which can be potential used for conformal cooling design, but also creates the twisted cooling channel. Even if the flow is supposed as fully laminar during optimization, the heated water mixed well with cold water to achieve a better cooling performance. In the case of a real industrial tool, the design process would start with fluid-thermal optimization for the cooling system design, followed by applying thermal-mechanical optimization for the solid
The optimized cooling channel has been compared with the conventional cooling channel (Fig. 9A) though a CFD model developed in COMSOL Multiphysics. As the optimized channel has a larger width compared to the conventional channel, a rectangular cooling channel (Fig. 9B) is included in the comparison. All these three cooling channels are designed with the same fluid volume. The steady $k-\epsilon$ turbulent model is chosen to simulate the real situation [12]. The domain is meshed using tetrahedral elements with 8 layers of prismatic element in the fluid domain. After a mesh convergence study, the global element size in the fluid is set to 1 mm.

<figure><img src="image_1697.png" /><figcaption>Fig. 9 Optimal solution obtained by fluid-thermal coupled topology optimization.</figcaption></figure>

Fig. 11 shows a comparison between three designs in terms of temperature evolution along the contact surface centerline. It can be observed that the conventional channel (circular shape) has an overall higher temperature along the surface centerline. The temperature with rectangular channel increases almost linearly. The highest temperatures are all located near the end of the surface due to the geometry. Though, the temperature with optimized channel is better controlled. The maximum temperature obtained with the optimized channel is more than 5 °C lower than those obtained with the other designs. Considering this is a simplified case, the gain in temperature may be much higher while facing to a more complex real tool.

<figure><img src="image_1698.png" /><figcaption>Fig. 10 Temperature comparison with three cooling channels.</figcaption></figure>

## Summary

In this study, topology optimization has shown its high suitability for designing hot stamping tools. In thermal-mechanical coupled optimization, around 40% of the solid material in a simplified Omega tool, without changing the tool's cooling performance. Even with a 50% material reduction in material, the reduction in tool load bearing capacity is not critical. This may lead to the price reduction and the building time saving for hot stamping tools. In fluid-thermal coupled optimization, the algorithm shows well its ability for cooling system designs. It does not only send water to the contact surface, which can be potential used for conformal cooling design, but also creates the twisted cooling channel. Even if the flow is supposed as fully laminar during optimization, the heated water mixed well with cold water to achieve a better cooling performance. In the case of a real industrial tool, the design process would start with fluid-thermal optimization for the cooling system design, followed by applying thermal-mechanical optimization for the solid
volume reduction. Further studies will continue applying topology optimization algorithms on real hot stamping tools to answer the following questions:

- Printability is not considered in the optimization algorithms

- Mechanical resistance is not considered during the cooling system design

- Using Stokes flow in optimization may not represent the real flow distribution

## Acknowledgment

This work has been performed within the framework of FuseMetal joint laboratory (ArcelorMittal France, UTC, CNRS). This laboratory has received funding from the European Union (FEDER Région Hauts-de-France).

## References

[1] D. Chantzis, X. Liu, D.J. Politis et al. Review on additive manufacturing of tooling for hot stamping. Int J Adv Manuf Technol 109 (2020) 87-107. https://doi.org/10.1007/s00170-020-05622-1

[2] N. Asnafi, J. Rajalampi, D. Aspenberg et al. Production Tools Made by Additive Manufacturing Through Laser-based Powder Bed Fusion. Berg Huettenmaenn Monatsh 165 (2020) 125-136. https://doi.org/10.1007/s00501-020-00961-8

[3] K. Liu, A. Tovar, An efficient 3D topology optimization code written in Matlab. Struct Multidisc Optim 50 (2014) 1175-1196. https://doi.org/10.1007/s00158-014-1107-x

[4] F. Hecht. New development in FreeFem++. Journal of Numerical Mathematics. 20.10.1515/jnum-2012-0013.

[5] J. Alexandersen, C.S. Andreasen. A Review of Topology Optimisation for Fluid-Based Problems. Fluids 5(1):29 (2020). https://doi.org/10.3390/fluids5010029

[6] Z. Shan, M. Zhang, C. Jiang et al. Basic study on die cooling system of hot stamping process. ATDM 2010, Beijing, 2010, pp. 1-4, doi: 10.1049/cp.2010.1248.

[7] Y. Sato, K. Yaji, K. Izui et al. An Optimum Design Method for a Thermal-Fluid Device Incorporating Multiobjective Topology Optimization With an Adaptive Weighting Scheme. ASME. J. Mech. Des. 140(3) (2018) 031402. https://doi.org/10.1115/1.4038209

[8] M. Yu, X. Wang, J. Gu et al. A synergic topology optimization approach on distribution of cooling channels and diverse-intensity heat sources for liquid-cooled heat sink. Struct Multidisc Optim 65 (2022) 48. https://doi.org/10.1007/s00158-021-03113-9

[9] P. Jolivet, V. Dolean, F. Hecht et al. High performance domain decomposition methods on massively parallel architectures with freefem++. Journal of Numerical Mathematics, 20(3-4) (2012) 287-302. https://doi.org/10.1515/jnum-2012-0015

[10] S. Sun, P. Liebersbach, X. Qian. 3D Topology Optimization of Heat Sinks for Liquid Cooling. Applied Thermal Engineering. 178. 115540. 10.1016/j.applthermaleng.2020.115540.

[11] A. Saysroy, S. Eiamsa-ard. Enhancing convective heat transfer in laminar and turbulent flow regions using multi-channel twisted tape inserts. International Journal of Thermal Sciences. 121. 55-74. 10.1016/j.ijthermalsci.2017.07.002.

[12] P. Hu, B. He, L. Ying. Numerical investigation on cooling performance of hot stamping tool with various channel designs. Applied Thermal Engineering. 96. 10.1016/j.applthermaleng.2015.10.154.

[13] S. Balay, S. Abhyankar, MF. Adams, et al. PETSc Web page. https://www.mcs.anl.gov/petsc; 2019.
volume reduction. Further studies will continue applying topology optimization algorithms on real hot stamping tools to answer the following questions:

- Printability is not considered in the optimization algorithms

- Mechanical resistance is not considered during the cooling system design

- Using Stokes flow in optimization may not represent the real flow distribution

## Acknowledgment

This work has been performed within the framework of FuseMetal joint laboratory (ArcelorMittal France, UTC, CNRS). This laboratory has received funding from the European Union (FEDER Région Hauts-de-France).

## References

[1] D. Chantzis, X. Liu, D.J. Politis et al. Review on additive manufacturing of tooling for hot stamping. Int J Adv Manuf Technol 109 (2020) 87-107. https://doi.org/10.1007/s00170-020-05622-1

[2] N. Asnafi, J. Rajalampi, D. Aspenberg et al. Production Tools Made by Additive Manufacturing Through Laser-based Powder Bed Fusion. Berg Huettenmaenn Monatsh 165 (2020) 125-136. https://doi.org/10.1007/s00501-020-00961-8

[3] K. Liu, A. Tovar, An efficient 3D topology optimization code written in Matlab. Struct Multidisc Optim 50 (2014) 1175-1196. https://doi.org/10.1007/s00158-014-1107-x

[4] F. Hecht. New development in FreeFem++. Journal of Numerical Mathematics. 20.10.1515/jnum-2012-0013.

[5] J. Alexandersen, C.S. Andreasen. A Review of Topology Optimisation for Fluid-Based Problems. Fluids 5(1):29 (2020). https://doi.org/10.3390/fluids5010029

[6] Z. Shan, M. Zhang, C. Jiang et al. Basic study on die cooling system of hot stamping process. ATDM 2010, Beijing, 2010, pp. 1-4, doi: 10.1049/cp.2010.1248.

[7] Y. Sato, K. Yaji, K. Izui et al. An Optimum Design Method for a Thermal-Fluid Device Incorporating Multiobjective Topology Optimization With an Adaptive Weighting Scheme. ASME. J. Mech. Des. 140(3) (2018) 031402. https://doi.org/10.1115/1.4038209

[8] M. Yu, X. Wang, J. Gu et al. A synergic topology optimization approach on distribution of cooling channels and diverse-intensity heat sources for liquid-cooled heat sink. Struct Multidisc Optim 65 (2022) 48. https://doi.org/10.1007/s00158-021-03113-9

[9] P. Jolivet, V. Dolean, F. Hecht et al. High performance domain decomposition methods on massively parallel architectures with freefem++. Journal of Numerical Mathematics, 20(3-4) (2012) 287-302. https://doi.org/10.1515/jnum-2012-0015

[10] S. Sun, P. Liebersbach, X. Qian. 3D Topology Optimization of Heat Sinks for Liquid Cooling. Applied Thermal Engineering. 178. 115540. 10.1016/j.applthermaleng.2020.115540.

[11] A. Saysroy, S. Eiamsa-ard. Enhancing convective heat transfer in laminar and turbulent flow regions using multi-channel twisted tape inserts. International Journal of Thermal Sciences. 121. 55-74. 10.1016/j.ijthermalsci.2017.07.002.

[12] P. Hu, B. He, L. Ying. Numerical investigation on cooling performance of hot stamping tool with various channel designs. Applied Thermal Engineering. 96. 10.1016/j.applthermaleng.2015.10.154.

[13] S. Balay, S. Abhyankar, MF. Adams, et al. PETSc Web page. https://www.mcs.anl.gov/petsc; 2019.
# Surrogate modeling for multi-objective optimization in the high-precision production of LiDAR glass optics

VU Anh Tuan¹,ᵃ*, PARIA Hamidreza¹,ᵇ, GRUNWALD Tim¹,ᶜ and BERGS Thomas¹,²,ᵈ

¹Fraunhofer Institute for Production Technology IPT, Aachen, Germany

²Manufacturing Technology Institute (MTI), RWTH Aachen University, Aachen, Germany

ᵃanh.tuan.vu@ipt.fraunhofer.de, ᵇhamidreza.paria@ipt.fraunhofer.de,

ᶜtim.grunwald@ipt.fraunhofer.de, ᵈthomas.bergs@ipt.fraunhofer.de

**Keywords:** Nonisothermal Glass Molding, Glass Optics, LiDAR, Surrogate Modeling, Bayesian Optimization, Industry 4.0

**Abstract.** This study addresses the ever-increasing demands on glass optics for LiDAR systems in autonomous vehicles, highlighting the pivotal role of the recently developed Nonisothermal Glass Molding (NGM) in enabling the mass production of complex and precise glass optics. While NGM promises a significant advancement in cost- and energy-efficient solutions, achieving the requisite shape and form accuracy for high-precision optics remains a persistent challenge. The research focuses on expediting the development phase, presenting a methodology that strategically utilizes a sparse dataset for determining optimized molding parameters with a minimized number of experimental trials. Importantly, our method highlights the exceptional ability of a robust surrogate model to precisely predict the accuracy outputs of glass optics, strongly influenced by numerous input molding parameters of the NGM process. This significance emphasizes the surrogate model, which emerges as a promising alternative to inefficient traditional methods, such as time-consuming experiments or computation-intensive simulations, particularly in the realm of high-precision production for LiDAR glass optics. In contributing to optics manufacturing advancements, this study also aligns with contemporary trends in digitalization and Industry 4.0 within modern optics production, thereby fostering innovation in the automotive industry.

**Introduction**

The rapid evolution of autonomous vehicle technology has ushered in a new era of transportation, promising unprecedented safety and efficiency. At the forefront of this technological revolution is Light Detection and Ranging (LiDAR), a key sensing technology that enables vehicles to navigate and perceive their surroundings with remarkable precision. As LiDAR systems play a critical role in providing real-time, three-dimensional mapping of the vehicle's environment, the efficacy of the optical components is paramount to ensure accurate and rapid data acquisition. In the quest to enhance LiDAR systems, the focus on optical components has become pivotal, serving as a cornerstone for achieving heightened performance and reliability. Glass optics, known for outstanding mechanical strength, durability, and exceptionally high transparency, stands out as a compelling solution to meet the rigorous demands of LiDAR systems. Nevertheless, amidst the recent growth in LiDAR technology, glass optics manufacturers are facing rising demands, including increasing geometrical complexity, high form accuracy, surface roughness, exceptional optical functionalities, reduced energy consumption, and cost-effectiveness from the dynamic LiDAR market.

Addressing these challenges, the “Glass4AutoFuture” project at Fraunhofer IPT, Germany aimed to develop a cost- and energy-efficient manufacturing technology enabling the high-volume production of optical glass components essential for future car concepts. Nonisothermal Glass

<figure><img src="image_1700.png" /></figure>
# Surrogate modeling for multi-objective optimization in the high-precision production of LiDAR glass optics

VU Anh Tuan<fnref n="1" />,ᵃ*, PARIA Hamidreza<fnref n="1" />,ᵇ, GRUNWALD Tim<fnref n="1" />,ᶜ and BERGS Thomas<fnref n="1" />,<fnref n="1" />,<fnref n="2" />,ᵈ

**Keywords:** Nonisothermal Glass Molding, Glass Optics, LiDAR, Surrogate Modeling, Bayesian Optimization, Industry 4.0

**Abstract.** This study addresses the ever-increasing demands on glass optics for LiDAR systems in autonomous vehicles, highlighting the pivotal role of the recently developed Nonisothermal Glass Molding (NGM) in enabling the mass production of complex and precise glass optics. While NGM promises a significant advancement in cost- and energy-efficient solutions, achieving the requisite shape and form accuracy for high-precision optics remains a persistent challenge. The research focuses on expediting the development phase, presenting a methodology that strategically utilizes a sparse dataset for determining optimized molding parameters with a minimized number of experimental trials. Importantly, our method highlights the exceptional ability of a robust surrogate model to precisely predict the accuracy outputs of glass optics, strongly influenced by numerous input molding parameters of the NGM process. This significance emphasizes the surrogate model, which emerges as a promising alternative to inefficient traditional methods, such as time-consuming experiments or computation-intensive simulations, particularly in the realm of high-precision production for LiDAR glass optics. In contributing to optics manufacturing advancements, this study also aligns with contemporary trends in digitalization and Industry 4.0 within modern optics production, thereby fostering innovation in the automotive industry.

**Introduction**

The rapid evolution of autonomous vehicle technology has ushered in a new era of transportation, promising unprecedented safety and efficiency. At the forefront of this technological revolution is Light Detection and Ranging (LiDAR), a key sensing technology that enables vehicles to navigate and perceive their surroundings with remarkable precision. As LiDAR systems play a critical role in providing real-time, three-dimensional mapping of the vehicle's environment, the efficacy of the optical components is paramount to ensure accurate and rapid data acquisition. In the quest to enhance LiDAR systems, the focus on optical components has become pivotal, serving as a cornerstone for achieving heightened performance and reliability. Glass optics, known for outstanding mechanical strength, durability, and exceptionally high transparency, stands out as a compelling solution to meet the rigorous demands of LiDAR systems. Nevertheless, amidst the recent growth in LiDAR technology, glass optics manufacturers are facing rising demands, including increasing geometrical complexity, high form accuracy, surface roughness, exceptional optical functionalities, reduced energy consumption, and cost-effectiveness from the dynamic LiDAR market.

Addressing these challenges, the “Glass4AutoFuture” project at Fraunhofer IPT, Germany aimed to develop a cost- and energy-efficient manufacturing technology enabling the high-volume production of optical glass components essential for future car concepts. Nonisothermal Glass

<figure><img src="image_1699.png" /></figure>

<footnote n="1">¹Fraunhofer Institute for Production Technology IPT, Aachen, Germany</footnote>

<footnote n="2">²Manufacturing Technology Institute (MTI), RWTH Aachen University, Aachen, Germany</footnote>

<footnote>ᵃanh.tuan.vu@ipt.fraunhofer.de, ᵇhamidreza.paria@ipt.fraunhofer.de,</footnote>

<footnote>ᶜtim.grunwald@ipt.fraunhofer.de, ᵈthomas.bergs@ipt.fraunhofer.de</footnote>
Molding (NGM), a replication-based process, has emerged as a viable technology for manufacturing complex and precision glass optics on a mass production scale [1,2]. In NGM, a glass preform is compressed between two mold halves at high temperature, and after cooling, optics with the desired shape are achieved without additional post-processing steps like grinding and polishing [3]. With the ability to produce multiple optics in parallel and eliminating the need for postprocessing, NGM holds the potential for high-volume production [4].

Nevertheless, the primary challenge that remains with NGM technology is to achieve the desired shape and form accuracy required for high precision LiDAR optics. These characteristics are influenced by multiple factors, including highly nonlinear deformation behaviors arising from nonequilibrium viscoelasticity of glass, the heat transfers between glass and molding environments, and numerous forming parameters. Understanding the thermal-mechanical behaviors of glass during the forming process and optimizing the forming parameters are crucial to satisfy the requirements of the glass optics [5]. These tasks are far from trivial, given the multitude of parameters influencing the accuracy of the optics, making the process development through an experimental trial-and-error method inefficient.

This study presents a method for determining optimized molding parameters with a minimal number of experimental trials. The proposed solution involves the development of a robust surrogate model that enables the multi-objective optimization of glass molding process parameters using a sparse dataset. The procedure consists of two steps. First, a process simulation model was established, accounting for the highly nonlinear viscoelastic deformation behavior of glass and all heat transfer phenomena occurring in the forming process. The simulation model serves as a tool to generate data for developing the surrogate model. The subsequent step involves exploring different surrogate models, acquisition functions, and selection strategies to minimize iterations while converging to the global solution. The results demonstrate that the surrogate model accurately approximates the outputs features of the product quality, offering reliable estimations that can substitute the computationally expensive simulations or time-consuming experiments.

## Optics demonstration and experimental setup

To demonstrate the methodology proposed in this study, an optical lens was purposefully designed for an optics demonstration. Fig. 1 presents the lens design and underlies the accuracy requirements for the demonstration. The design features a plano-convex lens commonly utilized for laser beam collimation in numerous LiDAR laser scanner systems. The precision of the optical lens is crucial for ensuring the high performance of the LiDAR system, as any single photon loss in the laser beam is unaffordable. Consequently, the form accuracies of the plane and spherical surfaces of the lens emerge as critical considerations. For this reason, form deviation of the spherical surface and center thickness were specifically identified as the primary criteria for optimizing the NGM process.

<table><caption>Fig. 1. Design and specifications of the optical demonstration.</caption><tr><td>Material</td><td>B270i</td></tr><tr><td>Refractive index n<sub>d</sub></td><td>1.52308 +/- 0.001</td></tr><tr><td>Abbe number</td><td>58.57 +/- 0.0818</td></tr><tr><td>Birefringence</td><td>10 nm/cm</td></tr></table><table><tr><td colspan="2">Specifications according to DIN ISO 10110</td></tr><tr><td>Form accuracy</td><td>3/ 0.005</td></tr><tr><td>Center thickness</td><td>25.7 +/-0.01</td></tr><tr><td>Bubbles and other inclusions</td><td>1/3x0.16</td></tr><tr><td>Inhomogeneity and Striae</td><td>2/0;3</td></tr><tr><td>Surface imperfections</td><td>5/3x0.16</td></tr></table>
Molding (NGM), a replication-based process, has emerged as a viable technology for manufacturing complex and precision glass optics on a mass production scale [1,2]. In NGM, a glass preform is compressed between two mold halves at high temperature, and after cooling, optics with the desired shape are achieved without additional post-processing steps like grinding and polishing [3]. With the ability to produce multiple optics in parallel and eliminating the need for postprocessing, NGM holds the potential for high-volume production [4].

Nevertheless, the primary challenge that remains with NGM technology is to achieve the desired shape and form accuracy required for high precision LiDAR optics. These characteristics are influenced by multiple factors, including highly nonlinear deformation behaviors arising from nonequilibrium viscoelasticity of glass, the heat transfers between glass and molding environments, and numerous forming parameters. Understanding the thermal-mechanical behaviors of glass during the forming process and optimizing the forming parameters are crucial to satisfy the requirements of the glass optics [5]. These tasks are far from trivial, given the multitude of parameters influencing the accuracy of the optics, making the process development through an experimental trial-and-error method inefficient.

This study presents a method for determining optimized molding parameters with a minimal number of experimental trials. The proposed solution involves the development of a robust surrogate model that enables the multi-objective optimization of glass molding process parameters using a sparse dataset. The procedure consists of two steps. First, a process simulation model was established, accounting for the highly nonlinear viscoelastic deformation behavior of glass and all heat transfer phenomena occurring in the forming process. The simulation model serves as a tool to generate data for developing the surrogate model. The subsequent step involves exploring different surrogate models, acquisition functions, and selection strategies to minimize iterations while converging to the global solution. The results demonstrate that the surrogate model accurately approximates the outputs features of the product quality, offering reliable estimations that can substitute the computationally expensive simulations or time-consuming experiments.

## Optics demonstration and experimental setup

To demonstrate the methodology proposed in this study, an optical lens was purposefully designed for an optics demonstration. Fig. 1 presents the lens design and underlies the accuracy requirements for the demonstration. The design features a plano-convex lens commonly utilized for laser beam collimation in numerous LiDAR laser scanner systems. The precision of the optical lens is crucial for ensuring the high performance of the LiDAR system, as any single photon loss in the laser beam is unaffordable. Consequently, the form accuracies of the plane and spherical surfaces of the lens emerge as critical considerations. For this reason, form deviation of the spherical surface and center thickness were specifically identified as the primary criteria for optimizing the NGM process.

<table><caption>Fig. 1. Design and specifications of the optical demonstration.</caption><tr><td>Material</td><td>B270i</td></tr><tr><td>Refractive index n<sub>d</sub></td><td>1.52308 +/- 0.001</td></tr><tr><td>Abbe number</td><td>58.57 +/- 0.0818</td></tr><tr><td>Birefringence</td><td>10 nm/cm</td></tr></table><table><tr><td colspan="2">Specifications according to DIN ISO 10110</td></tr><tr><td>Form accuracy</td><td>3/ 0.005</td></tr><tr><td>Center thickness</td><td>25.7 +/-0.01</td></tr><tr><td>Bubbles and other inclusions</td><td>1/3x0.16</td></tr><tr><td>Inhomogeneity and Striae</td><td>2/0;3</td></tr><tr><td>Surface imperfections</td><td>5/3x0.16</td></tr></table>
Fig. 2 outlines the process chain, experimental setup, and the process diagram indicating the parameters designed for producing the chosen lens demonstration. In the initial step, multiple glass preforms were heated in an external oven, while the two mold halves were preheated at a constant temperature, referred to as the *mold temperature* ($T_M$). Once the preforms reached a predefined *glass temperature* ($T_G$), each preform was automatically transferred to the molding position by a robot arm. At this stage, the glass preform was compressed in the space between two mold halves at a constant *pressing velocity* ($v$), until a predefined *pressing force* ($F$), set to prevent glass breakage, was attained. The force was maintained for a certain period – the *pressing time* ($t_H$) – to allow for complete deformation of the glass into the mold cavity. Following the molding step, the molded glass optics were transferred to an annealing oven, where multiple molded glass lenses underwent sequential cooling steps, comprising three *cooling times* ($C_{t,i}$, $i=1-3$). Each cooling time was carefully defined to allow the release of stress remaining after the molding step, ultimately achieving the desired shape of the molded lens after cooling at room temperature.

<figure><img src="image_1701.png" /><figcaption>Fig. 2. Description of NGM process and the process parameters required for optimization.</figcaption></figure>

## Process modeling

As depicted in Fig. 2, the molding process involves a set of 8 parameters, comprising of the set temperature of glass for molding $T_G$, predefined mold temperature $T_M$, molding velocity $v$, pressing force $F$, molding time $t_H$, and three phases of cooling in an annealing furnace $C_{t,i}$ where $i = 1-3$. The product quality is highly sensitive to these parameters, where adjustments to individual parameters can cause significant alterations in the final shape and form accuracy of the molded lens. The sensitivity arises from the complex, nonlinear thermoviscoelastic behaviors exhibited by glass during the molding step [6–8], coupled with nonlinear shrinkage behaviors resulting from structural relaxation during the cooling step [9–11]. Consequently, the conventional trial-and-error method for adjusting process parameters is both demanding and labor-intensive.

To address these challenges, this study introduces a numerical modeling approach for the NGM process. The overall goal is to understand the correlation between the process parameters (inputs)
Fig. 2 outlines the process chain, experimental setup, and the process diagram indicating the parameters designed for producing the chosen lens demonstration. In the initial step, multiple glass preforms were heated in an external oven, while the two mold halves were preheated at a constant temperature, referred to as the *mold temperature* ($T_M$). Once the preforms reached a predefined *glass temperature* ($T_G$), each preform was automatically transferred to the molding position by a robot arm. At this stage, the glass preform was compressed in the space between two mold halves at a constant *pressing velocity* ($v$), until a predefined *pressing force* ($F$), set to prevent glass breakage, was attained. The force was maintained for a certain period – the *pressing time* ($t_H$) – to allow for complete deformation of the glass into the mold cavity. Following the molding step, the molded glass optics were transferred to an annealing oven, where multiple molded glass lenses underwent sequential cooling steps, comprising three *cooling times* ($C_{t,i}$, $i=1-3$). Each cooling time was carefully defined to allow the release of stress remaining after the molding step, ultimately achieving the desired shape of the molded lens after cooling at room temperature.

<figure><img src="image_1700.png" /><figcaption>Fig. 2. Description of NGM process and the process parameters required for optimization.</figcaption></figure>

## Process modeling

As depicted in Fig. 2, the molding process involves a set of 8 parameters, comprising of the set temperature of glass for molding $T_G$, predefined mold temperature $T_M$, molding velocity $v$, pressing force $F$, molding time $t_H$, and three phases of cooling in an annealing furnace $C_{t,i}$ where $i = 1-3$. The product quality is highly sensitive to these parameters, where adjustments to individual parameters can cause significant alterations in the final shape and form accuracy of the molded lens. The sensitivity arises from the complex, nonlinear thermoviscoelastic behaviors exhibited by glass during the molding step [6–8], coupled with nonlinear shrinkage behaviors resulting from structural relaxation during the cooling step [9–11]. Consequently, the conventional trial-and-error method for adjusting process parameters is both demanding and labor-intensive.

To address these challenges, this study introduces a numerical modeling approach for the NGM process. The overall goal is to understand the correlation between the process parameters (inputs)
and the accuracy criteria (outputs), facilitating a multi-objective optimization method that minimizes the required dataset, given the computational expense of process simulation.

Fig. 3 presents the simulation model constructed to capture the thermal- and mechanical responses throughout the molding process. A coupled thermal-mechanical model was, therefore, necessary. The thermal model calculates the temperature distributions within the molding components, accounting for their thermal interactions with the molding environment. All relevant heat transfer phenomena at the boundaries, including convection and radiation, were considered. Convective and radiative heat transfer coefficients were defined through the utilization of infrared thermographic techniques, with the methodology and inverse heat transfer solution detailed in [12]. In addition, the heat transfer across the glass-mold interface, governed by the thermal contact conductance depending on the pressure *p*, interfacial temperature *T*, and surface roughness *Ra*, was adopted from prior works [13,14].

<figure><img src="image_1702.png" /><figcaption>Fig. 3. (a) Setup and molding components, and (b) coupled thermal-mechanical model of the NGM process [12].</figcaption></figure>

The mechanical model addresses the material responses of glass, including the stress relaxation and structural relaxation phenomena during molding and subsequent cooling. In the nonisothermal molding process, there are constantly heat exchanges at the glass surfaces with the colder molds and the molding environment, leading to a transition in the glass state from equilibrium to nonequilibrium [15]. Capturing the relaxation nature of the nonequilibrium glass necessitates a robust material model incorporating three variables: time, temperature, and thermal history [15]. A material model, developed in previous work [9], was adapted here to address the stress relaxation behaviors over the entire range of molding temperature, including thermo-viscoelasticity during the molding step and the shrinkage behaviors of glass due to the structural relaxation during the annealing step. Furthermore, the mechanical interaction at the glass-mold interface was modeled by a Coulomb's friction model. The friction coefficients, dependent on temperature *T*, pressure *p*, and sliding velocity *v*, were adopted from the previous work [16].

## Methodology

Conventional optimization methods commonly rely on 1st or 2nd order gradients of objective function or employ bio-inspired optimization algorithms such as evolutionary algorithms, which necessitate a large number of objective function evaluations to identify the global solution. In contrast, Bayesian Optimization (BO) stands out as a derivative-free machine learning method that endeavors to comprehend the intricate black box function of the system with limited and sparse dataset availability. Notably, BO exhibits efficiency in handling objective functions containing diverse systematic noises.
and the accuracy criteria (outputs), facilitating a multi-objective optimization method that minimizes the required dataset, given the computational expense of process simulation.

Fig. 3 presents the simulation model constructed to capture the thermal- and mechanical responses throughout the molding process. A coupled thermal-mechanical model was, therefore, necessary. The thermal model calculates the temperature distributions within the molding components, accounting for their thermal interactions with the molding environment. All relevant heat transfer phenomena at the boundaries, including convection and radiation, were considered. Convective and radiative heat transfer coefficients were defined through the utilization of infrared thermographic techniques, with the methodology and inverse heat transfer solution detailed in [12]. In addition, the heat transfer across the glass-mold interface, governed by the thermal contact conductance depending on the pressure *p*, interfacial temperature *T*, and surface roughness *Ra*, was adopted from prior works [13,14].

<figure><img src="image_1701.png" /><figcaption>Fig. 3. (a) Setup and molding components, and (b) coupled thermal-mechanical model of the NGM process [12].</figcaption></figure>

The mechanical model addresses the material responses of glass, including the stress relaxation and structural relaxation phenomena during molding and subsequent cooling. In the nonisothermal molding process, there are constantly heat exchanges at the glass surfaces with the colder molds and the molding environment, leading to a transition in the glass state from equilibrium to nonequilibrium [15]. Capturing the relaxation nature of the nonequilibrium glass necessitates a robust material model incorporating three variables: time, temperature, and thermal history [15]. A material model, developed in previous work [9], was adapted here to address the stress relaxation behaviors over the entire range of molding temperature, including thermo-viscoelasticity during the molding step and the shrinkage behaviors of glass due to the structural relaxation during the annealing step. Furthermore, the mechanical interaction at the glass-mold interface was modeled by a Coulomb's friction model. The friction coefficients, dependent on temperature *T*, pressure *p*, and sliding velocity *v*, were adopted from the previous work [16].

## Methodology

Conventional optimization methods commonly rely on 1st or 2nd order gradients of objective function or employ bio-inspired optimization algorithms such as evolutionary algorithms, which necessitate a large number of objective function evaluations to identify the global solution. In contrast, Bayesian Optimization (BO) stands out as a derivative-free machine learning method that endeavors to comprehend the intricate black box function of the system with limited and sparse dataset availability. Notably, BO exhibits efficiency in handling objective functions containing diverse systematic noises.
**Multi-Objective Bayesian Optimization (MOBO).** The versatility of BO becomes evident in multi-objective optimizations (MOO), where it diverges from single-objective optimization algorithms. Instead of converging to a single solution, MOO using BO converges to the Pareto optimal set of solutions within the design space of input features. This is especially pertinent when dealing with conflicting objective features, where enhancing one feature inevitably leads to the degradation of others. The set of corresponding output values for the Pareto optimal solutions is referred to as the *Pareto Front* (PF) set for a collection of conflicting objective parameters [17,18].

Building on these principles, the present study employed a MOBO framework, illustrated in Fig. 4. It encompasses four key steps: an initial sampling method, the construction of a surrogate model, the definition of an acquisition function, and the implementation of a selection strategy. This holistic approach allows MOBO to effectively navigate the complexities inherent in multi-objective optimization scenarios.

<figure><img src="image_1703.png" /><figcaption>Fig. 4. Flowchart of the multi-objective Bayesian optimization.</figcaption></figure>

**Surrogate model.** A surrogate model, a supervised machine learning approach, emulates complex and costly black box functions with high accuracy [19]. In this study, Gaussian Process (GP) was chosen among various supervised models like artificial neural networks (ANN) and tree-based methods due to its simplicity, fewer hyperparameters, and reduced sensitivity to the initial design of experiments [20]. GP proves advantageous in handling sparse datasets, avoiding overfitting issues that more complex models like ANN may encounter when trained on limited data. In this study, the GP surrogate model is trained independently for each objective parameter.

GP training: During the initial GP training on the black box function $f: \Omega \to \mathbb{R}^k$ with “k” elements of objectives, and design space of $\Omega \in \mathbb{R}^d$, a prior model was defined by mean function $m(\bar{x}): \Omega \to \mathbb{R}$, and Matern kernel (covariance) function $k(\bar{x}, \bar{x}^): \Omega \times \Omega \to \mathbb{R}$. This prior model was established without observing available data, forming a prior belief over the design space. For any given input variable $\bar{x}$, the prior of GP was expressed as:

$$f(\bar{x}) \sim N(m(\bar{x}), k(\bar{x}, \bar{x})). \qquad (1)$$
**Multi-Objective Bayesian Optimization (MOBO).** The versatility of BO becomes evident in multi-objective optimizations (MOO), where it diverges from single-objective optimization algorithms. Instead of converging to a single solution, MOO using BO converges to the Pareto optimal set of solutions within the design space of input features. This is especially pertinent when dealing with conflicting objective features, where enhancing one feature inevitably leads to the degradation of others. The set of corresponding output values for the Pareto optimal solutions is referred to as the *Pareto Front* (PF) set for a collection of conflicting objective parameters [17,18].

Building on these principles, the present study employed a MOBO framework, illustrated in Fig. 4. It encompasses four key steps: an initial sampling method, the construction of a surrogate model, the definition of an acquisition function, and the implementation of a selection strategy. This holistic approach allows MOBO to effectively navigate the complexities inherent in multi-objective optimization scenarios.

<figure><img src="image_1702.png" /><figcaption>Fig. 4. Flowchart of the multi-objective Bayesian optimization.</figcaption></figure>

**Surrogate model.** A surrogate model, a supervised machine learning approach, emulates complex and costly black box functions with high accuracy [19]. In this study, Gaussian Process (GP) was chosen among various supervised models like artificial neural networks (ANN) and tree-based methods due to its simplicity, fewer hyperparameters, and reduced sensitivity to the initial design of experiments [20]. GP proves advantageous in handling sparse datasets, avoiding overfitting issues that more complex models like ANN may encounter when trained on limited data. In this study, the GP surrogate model is trained independently for each objective parameter.

GP training: During the initial GP training on the black box function $f: \Omega \to \mathbb{R}^k$ with “k” elements of objectives, and design space of $\Omega \in \mathbb{R}^d$, a prior model was defined by mean function $m(\bar{x}): \Omega \to \mathbb{R}$, and Matern kernel (covariance) function $k(\bar{x}, \bar{x}^*): \Omega \times \Omega \to \mathbb{R}$. This prior model was established without observing available data, forming a prior belief over the design space. For any given input variable $\bar{x}$, the prior of GP was expressed as:

$$f(\bar{x}) \sim N(m(\bar{x}), k(\bar{x}, \bar{x}^*)) \quad (1)$$
Posterior GP training: Subsequently, a posterior of the GP was trained to refine the prior model using the available dataset, $D_n = \{(\bar{X}, \bar{Y})\} = \{(\bar{X}_1, \bar{Y}_1), (\bar{X}_2, \bar{Y}_2), ..., (\bar{X}_n, \bar{Y}_n)\}$ expressed as:

$$f(\bar{x}) \sim N(\mu(\bar{x}), \sigma(\bar{x})), \mu(\bar{x}) = m(\bar{x}) + kK^{-1}\bar{Y}, \sigma(\bar{x}) = k(\bar{x}, \bar{x}) - kK^{-1}k^T, \\ \text{where } k=k(\bar{x}, \bar{X}), \text{ and } K=k(\bar{X}, \bar{X}). \qquad (2)$$

Details for the surrogate model, prior, and posterior functions are elaborated in [17,21,22].

**Acquisition function:** Once the surrogate model $G_k$ was constructed to estimate each objective function $f_k$, $k \in N$ on available dataset, the acquisition function, serving as a cost-effective function to evaluate, was simply defined using the mean $\mu(\bar{x})$ and the uncertainty $\sigma(\bar{x})$ of the posterior function of the surrogate model:

$$u(\bar{x}|D_n) = u(\mu(\bar{x}), \sigma(\bar{x})). \qquad (3)$$

The primary objective of the acquisition function is to suggest the most promising single or batch of samples to be added to the dataset in the next iteration. This involves striking a balance between exploitation, focusing on the region with the highest mean value of the negative objective $\mu(\bar{x})$, and exploration, exploring areas with high objective function uncertainty $\sigma(\bar{x})$. The global maximum solution of the acquisition function, representing the trade-off between exploitation and exploration, guides the selection of new samples in the next iteration:

$$\bar{x}^{\text{new}} = \arg \max_{\bar{x} \in \Omega} u(\bar{x}|D_n). \qquad (4)$$

Examples of well-known acquisition functions include Expected Improvement (EI), Probability of Improvement (PI), Upper Confidence Bound (UCB), and Thompson sampling [18]. In this study, the Identity function, simply equal to mean $\mu(.)$ of $G_k$, is used as an acquisition function.

**Selection strategy:** After deriving the optimum points from the acquisition function, there are two methods for selecting sample(s) for the next iteration of the algorithm: sequential selection and batch selection. In sequential selection, one optimum point is chosen as the candidate for the next iteration in each cycle. While limiting the number of newly added samples to the dataset, it results in a high number of algorithm iterations and a low convergence rate. In contrast, the batch selection method suggests a batch of samples, improving the convergence rate of the algorithm. This strategy is particularly suitable when the black box function can be evaluated in parallel.

Three established methods for choosing the batch of samples for the upcoming iteration of MOBO are Kriging Believer (KB), Local Penalization (LP), and Believer-Penalizer (BP) [20–22].

## Results

The optimization process involves a set of 8 parameters. Table 1 provides an overview of the process parameters along with the selected range of physical constraints for each parameter. The table serves as a comprehensive reference, outlining the specific boundaries within which the MOBO framework operates to fine-tune the system and achieve optimal outcomes.

**Automated optimization procedure.** To automate the entire optimization procedure, the FEM model representing the entire NGM process, was treated as the black box function. This model, scripted in Python, seamlessly integrated with the MOBO algorithm. This integration allows the MOBO algorithm to generate any required simulation run through a simple command. The process parameters in Table 1 were defined as functional variables in Python, enabling the MOBO algorithm to efficiently modify and submit relevant simulations during its iterative runs. After each FEM simulation run's completion, essential output data was extracted from the .DAT file generated by ABAQUS and transmitted to MOBO as output parameter values.
Posterior GP training: Subsequently, a posterior of the GP was trained to refine the prior model using the available dataset, $D_n = \{(\bar{X}, \bar{Y})\} = \{(\bar{X}_1, \bar{Y}_1), (\bar{X}_2, \bar{Y}_2), ..., (\bar{X}_n, \bar{Y}_n)\}$ expressed as:

$$f(\bar{x}) \sim N(\mu(\bar{x}), \sigma(\bar{x})), \mu(\bar{x}) = m(\bar{x}) + kK^{-1}\bar{Y}, \sigma(\bar{x}) = k(\bar{x}, \bar{x}) - kK^{-1}k^T, \\ \text{where } k=k(\bar{x}, \bar{X}), \text{ and } K=k(\bar{X}, \bar{X}). \qquad (2)$$

Details for the surrogate model, prior, and posterior functions are elaborated in [17,21,22].

**Acquisition function:** Once the surrogate model $G_k$ was constructed to estimate each objective function $f_k, k \in N$ on available dataset, the acquisition function, serving as a cost-effective function to evaluate, was simply defined using the mean $\mu(\bar{x})$ and the uncertainty $\sigma(\bar{x})$ of the posterior function of the surrogate model:

$$u(\bar{x}|D_n) = u(\mu(\bar{x}), \sigma(\bar{x})). \qquad (3)$$

The primary objective of the acquisition function is to suggest the most promising single or batch of samples to be added to the dataset in the next iteration. This involves striking a balance between exploitation, focusing on the region with the highest mean value of the negative objective $\mu(\bar{x})$, and exploration, exploring areas with high objective function uncertainty $\sigma(\bar{x})$. The global maximum solution of the acquisition function, representing the trade-off between exploitation and exploration, guides the selection of new samples in the next iteration:

$$\bar{x}^{\text{new}} = \arg \max_{\bar{x} \in \Omega} u(\bar{x}|D_n). \qquad (4)$$

Examples of well-known acquisition functions include Expected Improvement (EI), Probability of Improvement (PI), Upper Confidence Bound (UCB), and Thompson sampling [18]. In this study, the Identity function, simply equal to mean $\mu(.)$ of $G_k$, is used as an acquisition function.

**Selection strategy:** After deriving the optimum points from the acquisition function, there are two methods for selecting sample(s) for the next iteration of the algorithm: sequential selection and batch selection. In sequential selection, one optimum point is chosen as the candidate for the next iteration in each cycle. While limiting the number of newly added samples to the dataset, it results in a high number of algorithm iterations and a low convergence rate. In contrast, the batch selection method suggests a batch of samples, improving the convergence rate of the algorithm. This strategy is particularly suitable when the black box function can be evaluated in parallel.

Three established methods for choosing the batch of samples for the upcoming iteration of MOBO are Kriging Believer (KB), Local Penalization (LP), and Believer-Penalizer (BP) [20–22].

## Results

The optimization process involves a set of 8 parameters. Table 1 provides an overview of the process parameters along with the selected range of physical constraints for each parameter. The table serves as a comprehensive reference, outlining the specific boundaries within which the MOBO framework operates to fine-tune the system and achieve optimal outcomes.

**Automated optimization procedure.** To automate the entire optimization procedure, the FEM model representing the entire NGM process, was treated as the black box function. This model, scripted in Python, seamlessly integrated with the MOBO algorithm. This integration allows the MOBO algorithm to generate any required simulation run through a simple command. The process parameters in Table 1 were defined as functional variables in Python, enabling the MOBO algorithm to efficiently modify and submit relevant simulations during its iterative runs. After each FEM simulation run's completion, essential output data was extracted from the .DAT file generated by ABAQUS and transmitted to MOBO as output parameter values.
<table><caption>Table 1. List of process parameters with their physical ranges defined for NGM process.</caption><thead><tr><th>Parameter</th><th>Denotation</th><th>Unit</th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>Glass temperature</td><td>T<sub>G</sub></td><td>[°C]</td><td>750.0</td><td>850.0</td></tr><tr><td>Molds temperature</td><td>T<sub>M</sub></td><td>[°C]</td><td>480.0</td><td>580.0</td></tr><tr><td>Molding velocity</td><td>v</td><td>[m/s]</td><td>1e-3</td><td>2e-2</td></tr><tr><td>Pressing force</td><td>F</td><td>[kN]</td><td>5.0</td><td>30.0</td></tr><tr><td>Pressing time</td><td>t<sub>H</sub></td><td>[s]</td><td>5.0</td><td>30.0</td></tr><tr><td>Cooling time 1</td><td>C<sub>t,1</sub></td><td>[s]</td><td>1800.0</td><td>18000.0</td></tr><tr><td>Cooling time 2</td><td>C<sub>t,2</sub></td><td>[s]</td><td>1680.0</td><td>168000.0</td></tr><tr><td>Cooling time 3</td><td>C<sub>t,3</sub></td><td>[s]</td><td>2742.0</td><td>13710.0</td></tr></tbody></table>

For MOBO implementation, the open-source Python library AutoOED was employed [21], and its hyperparameters were tuned following the methodology in [22]. The primary objective was to minimize two parameters of the black box function, representing the required accuracy criteria of the glass lenses. These objectives include: (i) optimizing the form accuracy of the molded lenses, quantified as the Root Mean Square Error (RMSE) of the convex surfaces between the desired shape y and that estimated by the FEM simulation ŷ, as per Eq. 5, and (ii) minimizing the absolute error between the center thickness of the lens and the target value given as 25.7 mm.

$$RMSE(\hat{y}, y) = \sqrt{\sum_{i=1}^{n} \frac{(\hat{y}_i - y_i)^2}{n}} \quad (5)$$

As presented in Fig. 4, the MOBO algorithm commenced with the generation of 20 quasi-random sample distributions along each input feature using Latin Hypercube Sampling (LHS) [23], denoted as $\bar{X}_0 = (X_1, X_2, ..., X_{20}) \in R^{20 \times 8}$. These samples constituted the initial design space. After running 20 simulations for each input sample, the objective parameters were evaluated, resulting in $\bar{Y}_0 = (Y_1, Y_2, ..., Y_{20}) \in R^{20 \times 2}$, referred to as the initial performance space. With both the initial design and performance spaces created, the MOBO algorithm was initiated, and through each iteration, 20 new samples were added to the dataset. These samples underwent FEM simulation to obtain the associated output parameters. This process continued until one of the predefined stopping criteria was attained. Two stopping criteria were defined: the first criterion was satisfied if the form accuracy error and center thickness error reach values below 5 µm and 0.1 mm, respectively. The second criterion was when the total number of iterations reaches 10, equivalent to a dataset of 200 samples.

The output results for each iteration were visually represented in a performance space plot, where the objective features were plotted against each other. Fig. 5(a) exhibits the outcomes of the LHS, while Fig. 5 (b) illustrates the results of the MOBO after 7 iterations. The upper plots showcase the full performance spaces, and the lower plots zoom in on the rectangular region delineated by a dash-dot line in the upper plot. The green region in the lower left corner of both lower plots represents the area where the first objective criterion was fulfilled. The grey zones indicate areas where only one of the objective features satisfied the criterion, while the orange zone signifies that none of the targets were achieved. Given that 20 samples were generated in each iteration, Fig. 5(b) only shows the best sample achieved after each iteration, considering both objectives. This selected set is termed the *Pareto Front* (PF) for that specific iteration.
<table><caption>Table 1. List of process parameters with their physical ranges defined for NGM process.</caption><thead><tr><th>Parameter</th><th>Denotation</th><th>Unit</th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>Glass temperature</td><td>T<sub>G</sub></td><td>[°C]</td><td>750.0</td><td>850.0</td></tr><tr><td>Molds temperature</td><td>T<sub>M</sub></td><td>[°C]</td><td>480.0</td><td>580.0</td></tr><tr><td>Molding velocity</td><td>v</td><td>[m/s]</td><td>1e-3</td><td>2e-2</td></tr><tr><td>Pressing force</td><td>F</td><td>[kN]</td><td>5.0</td><td>30.0</td></tr><tr><td>Pressing time</td><td>t<sub>H</sub></td><td>[s]</td><td>5.0</td><td>30.0</td></tr><tr><td>Cooling time 1</td><td>C<sub>t,1</sub></td><td>[s]</td><td>1800.0</td><td>18000.0</td></tr><tr><td>Cooling time 2</td><td>C<sub>t,2</sub></td><td>[s]</td><td>1680.0</td><td>168000.0</td></tr><tr><td>Cooling time 3</td><td>C<sub>t,3</sub></td><td>[s]</td><td>2742.0</td><td>13710.0</td></tr></tbody></table>

For MOBO implementation, the open-source Python library AutoOED was employed [21], and its hyperparameters were tuned following the methodology in [22]. The primary objective was to minimize two parameters of the black box function, representing the required accuracy criteria of the glass lenses. These objectives include: (i) optimizing the form accuracy of the molded lenses, quantified as the Root Mean Square Error (RMSE) of the convex surfaces between the desired shape y and that estimated by the FEM simulation ŷ, as per Eq. 5, and (ii) minimizing the absolute error between the center thickness of the lens and the target value given as 25.7 mm.

$$ \text{RMSE}(\hat{\mathbf{y}}, \mathbf{y}) = \sqrt{\sum_{i=1}^{n} \frac{(\hat{y}_i - y_i)^2}{n}}. \qquad (5) $$

As presented in Fig. 4, the MOBO algorithm commenced with the generation of 20 quasi-random sample distributions along each input feature using Latin Hypercube Sampling (LHS) [23], denoted as $\bar{X}_0 = (X_1, X_2, \dots, X_{20}) \in \mathbb{R}^{20 \times 8}$. These samples constituted the initial design space. After running 20 simulations for each input sample, the objective parameters were evaluated, resulting in $\bar{Y}_0 = (Y_1, Y_2, \dots, Y_{20}) \in \mathbb{R}^{20 \times 2}$, referred to as the initial performance space. With both the initial design and performance spaces created, the MOBO algorithm was initiated, and through each iteration, 20 new samples were added to the dataset. These samples underwent FEM simulation to obtain the associated output parameters. This process continued until one of the predefined stopping criteria was attained. Two stopping criteria were defined: the first criterion was satisfied if the form accuracy error and center thickness error reach values below 5 µm and 0.1 mm, respectively. The second criterion was when the total number of iterations reaches 10, equivalent to a dataset of 200 samples.

The output results for each iteration were visually represented in a performance space plot, where the objective features were plotted against each other. Fig. 5(a) exhibits the outcomes of the LHS, while Fig. 5 (b) illustrates the results of the MOBO after 7 iterations. The upper plots showcase the full performance spaces, and the lower plots zoom in on the rectangular region delineated by a dash-dot line in the upper plot. The green region in the lower left corner of both lower plots represents the area where the first objective criterion was fulfilled. The grey zones indicate areas where only one of the objective features satisfied the criterion, while the orange zone signifies that none of the targets were achieved. Given that 20 samples were generated in each iteration, Fig. 5(b) only shows the best sample achieved after each iteration, considering both objectives. This selected set is termed the *Pareto Front* (PF) for that specific iteration.
<figure><img src="image_1704.png" /><figcaption>Fig. 5. Results of the MOBO algorithm on FEM simulation of the NGM process.</figcaption></figure>

The inherent randomness of LHS is evident in Fig. 5(a), as data samples spread far from the defined target values, and only one best sample from the first generation is located in the grey zone. However, Fig. 5(b) demonstrates a significant improvement in objective values obtained by iteration 4, with the algorithm meeting the predefined criteria for the first time in iteration 5. Continuing for two more iterations, MOBO reveals that not only can global solutions be enhanced, but potential local minima in the design space can also be uncovered, thanks to its exploration and exploitation features. In iteration 7, 12 more samples lie in the grey zone, highlighting MOBO's capability not only in identifying global optimum points but also in intelligently designing experiments by detecting possible local minima, thus reducing the required number of data samples compared to conventional Design of Experiment (DoE) algorithms.

To assess the prediction accuracy of GP as the surrogate model in MOBO, the absolute error between the estimated and target values is given for each sample in Fig. 6. After 7 iterations, encompassing 160 samples including the initial LHS, the error plots for both output features nearly reach zero. The result reveals that with this sample size, the surrogate model accurately approximates the output features, offering reliable estimations that can substitute the computationally expensive FEM simulations. In addition, the GP prediction error significantly diminishes after generating 80 samples (4 iterations). This marks the first instance where MOBO detects the vicinity of the global optimum solution, as illustrated in Fig. 5. Despite slight fluctuations in the center thickness error plot during iteration 5, the algorithm effectively refines the global solution, emphasizing the robustness and efficacy of the MOBO methodology.
<figure><img src="image_1703.png" /><figcaption>Fig. 5. Results of the MOBO algorithm on FEM simulation of the NGM process.</figcaption></figure>

The inherent randomness of LHS is evident in Fig. 5(a), as data samples spread far from the defined target values, and only one best sample from the first generation is located in the grey zone. However, Fig. 5(b) demonstrates a significant improvement in objective values obtained by iteration 4, with the algorithm meeting the predefined criteria for the first time in iteration 5. Continuing for two more iterations, MOBO reveals that not only can global solutions be enhanced, but potential local minima in the design space can also be uncovered, thanks to its exploration and exploitation features. In iteration 7, 12 more samples lie in the grey zone, highlighting MOBO's capability not only in identifying global optimum points but also in intelligently designing experiments by detecting possible local minima, thus reducing the required number of data samples compared to conventional Design of Experiment (DoE) algorithms.

To assess the prediction accuracy of GP as the surrogate model in MOBO, the absolute error between the estimated and target values is given for each sample in Fig. 6. After 7 iterations, encompassing 160 samples including the initial LHS, the error plots for both output features nearly reach zero. The result reveals that with this sample size, the surrogate model accurately approximates the output features, offering reliable estimations that can substitute the computationally expensive FEM simulations. In addition, the GP prediction error significantly diminishes after generating 80 samples (4 iterations). This marks the first instance where MOBO detects the vicinity of the global optimum solution, as illustrated in Fig. 5. Despite slight fluctuations in the center thickness error plot during iteration 5, the algorithm effectively refines the global solution, emphasizing the robustness and efficacy of the MOBO methodology.
<figure><img src="image_1705.png" /><figcaption>Fig. 6. The absolute error between the outputs predicted by Gaussian process and those given by FEM simulation.</figcaption></figure>

**Summary**

Tackling the challenge of finding the global solution while minimizing operational time is a fundamental concern in optimizing processes with multiple variables, particularly in precision-demanding applications like glass molding. While experimental operations are time-consuming and labor-intensive, simulating glass molding processes requires huge computational resources. Our research enlightens the effectiveness of a MOBO framework as a surrogate model, offering an alternative approach to predict product accuracy based on multiple molding parameters. The methodology holds promising implications for industry applications. By conducting a limited number of trial experiments during the process development phase, industrial users can efficiently obtain optimized parameters aligned with specific product accuracy criteria, leading to accelerated phases of process development and production ramp-up. This study aligns with current trends in digitalization and Industry 4.0, contributing to innovation within the automotive industry.

**Acknowledgement**

This research was supported within the Fraunhofer and DFG transfer program under the project »Glass4AutoFuture«. Grant number: 032-600020.

**References**

[1] H. Kreilkamp, A.T. Vu, O. Dambon, N.F. Klocke, Non-Isothermal Glass Moulding of Complex Led Optics, in: S.K. Sundaram (Ed.), 77th Conference on Glass Problems, Wiley (2017) 141–149. https://doi.org/10.1002/9781119417507.ch13

[2] H. Kreilkamp, A.T. Vu, O. Dambon, F. Klocke, Replicative manufacturing of complex lighting optics by non-isothermal glass molding, in: D.H. Krevor, W.S. Beich, M.P. Schaub, A. Symmons (Eds.), Polymer Optics and Molded Glass Optics: Design, Fabrication, and Materials SPIE (2016) 99490B. https://doi.org/10.1117/12.2235848

[3] C. Strobl, P.A. Vogel, A.T. Vu, H. Mende, T. Grunwald, R.H. Schmitt, T. Bergs, Enabling Sustainability in Glass Optics Manufacturing by Wafer Scale Molding, KEM 926 (2022) 2371. https://doi.org/10.4028/p-hachrx

[4] A.T. Vu, H. Kreilkamp, O. Dambon, F. Klocke, Nonisothermal glass molding for the cost-efficient production of precision freeform optics, Opt. Eng 55 (2016) 71207. https://doi.org/10.1117/1.OE.55.7.071207

[5] A.T. Vu, H. Kreilkamp, B.J. Krishnamoorthi, O. Dambon, F. Klocke, A hybrid optimization approach in non-isothermal glass molding, in: Author(s) (2016) 40006.

[6] A.T. Vu, A.N. Vu, T. Grunwald, T. Bergs, Modeling of thermo-viscoelastic material behavior of glass over a wide temperature range in glass compression molding, J Am Ceram Soc 103 (2020) 2791–2807. https://doi.org/10.1111/jace.16963
<figure><img src="image_1704.png" /><figcaption>Fig. 6. The absolute error between the outputs predicted by Gaussian process and those given by FEM simulation.</figcaption></figure>

**Summary**

Tackling the challenge of finding the global solution while minimizing operational time is a fundamental concern in optimizing processes with multiple variables, particularly in precision-demanding applications like glass molding. While experimental operations are time-consuming and labor-intensive, simulating glass molding processes requires huge computational resources. Our research enlightens the effectiveness of a MOBO framework as a surrogate model, offering an alternative approach to predict product accuracy based on multiple molding parameters. The methodology holds promising implications for industry applications. By conducting a limited number of trial experiments during the process development phase, industrial users can efficiently obtain optimized parameters aligned with specific product accuracy criteria, leading to accelerated phases of process development and production ramp-up. This study aligns with current trends in digitalization and Industry 4.0, contributing to innovation within the automotive industry.

**Acknowledgement**

This research was supported within the Fraunhofer and DFG transfer program under the project »Glass4AutoFuture«. Grant number: 032-600020.

**References**

[1] H. Kreilkamp, A.T. Vu, O. Dambon, N.F. Klocke, Non-Isothermal Glass Moulding of Complex Led Optics, in: S.K. Sundaram (Ed.), 77th Conference on Glass Problems, Wiley (2017) 141–149. https://doi.org/10.1002/9781119417507.ch13

[2] H. Kreilkamp, A.T. Vu, O. Dambon, F. Klocke, Replicative manufacturing of complex lighting optics by non-isothermal glass molding, in: D.H. Krevor, W.S. Beich, M.P. Schaub, A. Symmons (Eds.), Polymer Optics and Molded Glass Optics: Design, Fabrication, and Materials SPIE (2016) 99490B. https://doi.org/10.1117/12.2235848

[3] C. Strobl, P.A. Vogel, A.T. Vu, H. Mende, T. Grunwald, R.H. Schmitt, T. Bergs, Enabling Sustainability in Glass Optics Manufacturing by Wafer Scale Molding, KEM 926 (2022) 2371. https://doi.org/10.4028/p-hachrx

[4] A.T. Vu, H. Kreilkamp, O. Dambon, F. Klocke, Nonisothermal glass molding for the cost-efficient production of precision freeform optics, Opt. Eng 55 (2016) 71207. https://doi.org/10.1117/1.OE.55.7.071207

[5] A.T. Vu, H. Kreilkamp, B.J. Krishnamoorthi, O. Dambon, F. Klocke, A hybrid optimization approach in non-isothermal glass molding, in: Author(s) (2016) 40006.

[6] A.T. Vu, A.N. Vu, T. Grunwald, T. Bergs, Modeling of thermo-viscoelastic material behavior of glass over a wide temperature range in glass compression molding, J Am Ceram Soc 103 (2020) 2791–2807. https://doi.org/10.1111/jace.16963
[7] G. Liu, A.T. Vu, O. Dambon, F. Klocke, Glass Material Modeling and its Molding Behavior, MRS Advances 2 (2017) 875-885. https://doi.org/10.1557/adv.2017.64

[8] F. Wang, Simulating the precision glass molding process. Zugl.: Aachen, Techn. Hochsch., Diss., 2013, 1st ed., Apprimus-Verl., Aachen, 2014.

[9] A.T. Vu, T. Grunwald, T. Bergs, Thermo-viscoelastic Modeling of Nonequilibrium Material Behavior of Glass in Nonisothermal Glass Molding, Procedia Manufacturing 47 (2020) 1561. https://doi.org/10.1016/j.promfg.2020.04.350

[10] T.D. Pallicity, A.T. Vu, K. Ramesh, P. Mahajan, G. Liu, O. Dambon, Birefringence measurement for validation of simulation of precision glass molding process, J Am Ceram Soc 100 (2017) 4680-4698. https://doi.org/10.1111/jace.15010

[11] F. Wang, Y. Chen, F. Klocke, G. Pongs, A.Y. Yi, Numerical Simulation Assisted Curve Compensation in Compression Molding of High Precision Aspherical Glass Lenses, Journal of Manufacturing Science and Engineering 131 (2009). https://doi.org/10.1115/1.3063652

[12] A.T. Vu, Modeling Relaxation Nature of Nonequilibrium Glass in Nonisothermal Glass Molding. Dissertation, 2023.

[13] A.T. Vu, A.N. Vu, G. Liu, T. Grunwald, O. Dambon, F. Klocke, T. Bergs, Experimental investigation of contact heat transfer coefficients in nonisothermal glass molding by infrared thermography, J Am Ceram Soc 102 (2019) 2116-2134. https://doi.org/10.1111/jace.16029

[14] A.T. Vu, T. Helmig, A.N. Vu, Y. Frekers, T. Grunwald, R. Kneer, T. Bergs, Numerical and experimental determinations of contact heat transfer coefficients in nonisothermal glass molding, J Am Ceram Soc 103 (2020) 1258-1269. https://doi.org/10.1111/jace.16756

[15] A.T. Vu, R.d.l.A. Avila Hernandez, T. Grunwald, T. Bergs, Modeling nonequilibrium thermoviscoelastic material behaviors of glass in nonisothermal glass molding, J Am Ceram Soc 105 (2022) 6799-6815. https://doi.org/10.1111/jace.18605

[16] A.T. Vu, T. Grunwald, T. Bergs, Friction in Glass Forming: Tribological Behaviors of Optical Glasses and uncoated Steel near Glass Transition Temperature, J Non-Cryst Solids (2024).

[17] B. Shahriari, K. Swersky, Z. Wang, R.P. Adams, N. de Freitas, Taking the Human Out of the Loop: A Review of Bayesian Optimization, Proc. IEEE 104 (2016) 148-175. https://doi.org/10.1109/JPROC.2015.2494218

[18] S. Greenhill, S. Rana, S. Gupta, P. Vellanki, S. Venkatesh, Bayesian Optimization for Adaptive Experimental Design: A Review, IEEE Access 8 (2020) 13937-13948. https://doi.org/10.1109/ACCESS.2020.2966228

[19] S. Ranftl, W. von der Linden, Bayesian Surrogate Analysis and Uncertainty Propagation, in: The 40th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, MDPI, Basel Switzerland (2021) 6. https://doi.org/10.3390/psf2021003006

[20] F. Mekki-Berrada, Z. Ren, T. Huang, W.K. Wong, F. Zheng, J. Xie, I.P.S. Tian et al., Two-step machine learning enables optimized nanoparticle synthesis, npj Comput Mater 7 (2021). https://doi.org/10.1038/s41524-021-00520-w

[21] Y. Tian, M.K. Luković, T. Erps, M. Foshey, W. Matusik, AutoOED: Automated Optimal Experiment Design Platform, arXiv, 2021.

[22] M.K. Luković, Y. Tian, W. Matusik, Diversity-Guided Multi-Objective Bayesian Optimization with Batch Evaluations, in: Proceedings of the 34th International Conference on Neural Information Processing Systems, Curran Associates Inc, Red Hook, NY, USA, 2020.

[23] P. Saves et al., SMT 2.0: A Surrogate Modeling Toolbox with a focus on hierarchical and mixed variables Gaussian processes, Advances in Engineering Software 188 (2024) 103571. https://doi.org/10.1016/j.advengsoft.2023.103571