<figure><img src="image_2108.png" /><figcaption>Figure 5: Comparison between the two method of strain measurement: by extensometer and by DIC for the FG microstructure.</figcaption></figure>

## Material model

Elasto-viscoplastic behavior: To model the behavior of the Ti-6Al-4V, an elasto-viscoplastic model is used. Strain-rate is partitioned between an elastic and a plastic strain-rate (Eq. 1). Elastic deformation is given according to Hooke's law (Eq. 2). The viscoplastic model chosen is a three- parameters hyperbolic sinus model (Eq. 3, Eq. 4). A two-parameter hardening law is used to describe the plasticity (Eq. 5). Finally, a von Mises yield surface mark the transition between elastic and viscoplastic deformation (Eq. 6).

$$ \bar{\epsilon} = \bar{\epsilon}_e + \bar{\epsilon}_p \qquad (1) $$

$$ \bar{\sigma} = \bar{C}\bar{\epsilon}_e \qquad (2) $$

$$ \bar{\epsilon}_p = \frac{3}{2} \dot{p} \frac{\bar{S}}{\sigma_{eq}} \qquad (3) $$

$$ \dot{p} = Asinh(\beta f)^n \qquad (4) $$

$$ R = Q(1 - e^{-bp}) \qquad (5) $$

$$ f = \sigma_{eq} - k - R \qquad (6) $$

## Model identification

Method description: the Young's modulus and elasticity limit are identified from the tensile test performed under small deformation condition. In order to identify the 3 viscous parameters, Eq.7 is written as a log-log equation. The relaxation test is used to describe the evolution of the viscous stress and the viscoplastic strain rate is calculated as a function of the stress rate and the Young's modulus (Eq. 8).

$$ \sigma_{vp} = \frac{1}{\beta} \arcsinh \left[ \left( \frac{1}{A} |\dot{\epsilon}_{vp}| \right)^{\frac{1}{n}} \right] \qquad (7) $$
<figure><img src="image_2099.png" /><figcaption>Figure 5: Comparison between the two method of strain measurement: by extensometer and by DIC for the FG microstructure.</figcaption></figure>

## Material model

Elasto-viscoplastic behavior: To model the behavior of the Ti-6Al-4V, an elasto-viscoplastic model is used. Strain-rate is partitioned between an elastic and a plastic strain-rate (Eq. 1). Elastic deformation is given according to Hooke's law (Eq. 2). The viscoplastic model chosen is a three- parameters hyperbolic sinus model (Eq. 3, Eq. 4). A two-parameter hardening law is used to describe the plasticity (Eq. 5). Finally, a von Mises yield surface mark the transition between elastic and viscoplastic deformation (Eq. 6).

$$ \bar{\varepsilon} = \bar{\varepsilon}_e + \bar{\varepsilon}_p \qquad (1) $$

$$ \bar{\sigma} = \bar{C} \bar{\varepsilon}_e \qquad (2) $$

$$ \bar{\varepsilon}_p = \frac{3}{2} \dot{p} \frac{\bar{S}}{\sigma_{eq}} \qquad (3) $$

$$ \dot{p} = Asinh(\beta f)^n \qquad (4) $$

$$ R = Q(1 - e^{-bp}) \qquad (5) $$

$$ f = \sigma_{eq} - k - R \qquad (6) $$

## Model identification

Method description: the Young's modulus and elasticity limit are identified from the tensile test performed under small deformation condition. In order to identify the 3 viscous parameters, Eq.7 is written as a log-log equation. The relaxation test is used to describe the evolution of the viscous stress and the viscoplastic strain rate is calculated as a function of the stress rate and the Young's modulus (Eq. 8).

$$ \sigma_{vp} = \frac{1}{\beta} \arcsinh \left[ \left( \frac{1}{A} |\dot{\varepsilon}_{vp}| \right)^{\frac{1}{n}} \right] \qquad (7) $$
$$ \dot{\epsilon}_{vp} = \frac{-\dot{\sigma}}{E} \qquad (8) $$

Once the viscoplastic parameters have been identified, the last parameters to determine are the two hardening parameters. A first identification is done using the tensile test performed under small deformation conditions, up to 0.02 deformation.

## Results

The evolutions of the 3 viscous parameters and the two hardening parameters with temperature is shown in Figure 6. Rollin et al. reported similar values and evolution for the viscous parameters [16]. The values of the hardening parameters b and Q (Figure 6b) show a decreasing evolution with temperature, the hardening is lower at higher temperature but also slower to reach its saturation value Q. The FG microstructure seems to present a higher hardening than the CG microstructure, at first contradicting the results of the relaxation tests. However, it should be considered that those hardening parameters identified give a faster saturation for the coarse-grained microstructure. No hardening is considered at 650 °C and above, and the parameters Q and b are set to 0. This identification does not yet include the traction to rupture test which, with higher strain levels, presents more hardening.

<figure><img src="image_2109.png" /><figcaption>Figure 6: Model parameters evolution with temperature for both microstructures</figcaption></figure>

After identification of the viscoplastic part, a simple isotropic hardening was fitted to complete the model. A comparison between the experimental data and the modelling of the tensile curves for the FG microstructure is presented in Figure 7. The predicted stresses at 0.02 strain is consistent with the experimental stresses. More discrepancy exists at the onset of hardening, with a more abrupt elasto-viscoplastic transition for the model which predict higher stress, particularly at the low temperatures. The maximum gap between experimental and model is quantified according to Eq. 9. The values fluctuate around a 10% difference, mostly at the strain value of the elasto-viscoplastic transition for each temperature.

$$ \sigma_{error} = \frac{|\sigma_{exp} - \sigma_{mod}|}{|\sigma_{exp}|} \qquad (9) $$
$$ \dot{\epsilon}_{vp} = \frac{-\dot{\sigma}}{E} \qquad (8) $$

Once the viscoplastic parameters have been identified, the last parameters to determine are the two hardening parameters. A first identification is done using the tensile test performed under small deformation conditions, up to 0.02 deformation.

## Results

The evolutions of the 3 viscous parameters and the two hardening parameters with temperature is shown in Figure 6. Rollin et al. reported similar values and evolution for the viscous parameters [16]. The values of the hardening parameters b and Q (Figure 6b) show a decreasing evolution with temperature, the hardening is lower at higher temperature but also slower to reach its saturation value Q. The FG microstructure seems to present a higher hardening than the CG microstructure, at first contradicting the results of the relaxation tests. However, it should be considered that those hardening parameters identified give a faster saturation for the coarse-grained microstructure. No hardening is considered at 650 °C and above, and the parameters Q and b are set to 0. This identification does not yet include the traction to rupture test which, with higher strain levels, presents more hardening.

<figure><img src="image_2100.png" /><figcaption>Figure 6: Model parameters evolution with temperature for both microstructures</figcaption></figure>

After identification of the viscoplastic part, a simple isotropic hardening was fitted to complete the model. A comparison between the experimental data and the modelling of the tensile curves for the FG microstructure is presented in Figure 7. The predicted stresses at 0.02 strain is consistent with the experimental stresses. More discrepancy exists at the onset of hardening, with a more abrupt elasto-viscoplastic transition for the model which predict higher stress, particularly at the low temperatures. The maximum gap between experimental and model is quantified according to Eq. 9. The values fluctuate around a 10% difference, mostly at the strain value of the elasto-viscoplastic transition for each temperature.

$$ \sigma_{error} = \frac{|\sigma_{exp} - \sigma_{mod}|}{|\sigma_{exp}|} \qquad (9) $$
<figure><img src="image_2110.png" /><figcaption>Figure 7: Comparison between the experimental loading curves and the model at each temperature for the FG microstructure</figcaption></figure>

**Conclusion**

This study investigated the mechanical behavior of two different microstructures of Ti-6Al-4V alloy at high temperature. The DIC method allowed to acquire information until very large strain (at least up to a true strain of 2). The results are in good agreement with the strain obtained using a contact method (extensometer).

An elasto-viscoplastic model has been identified using relaxation tests, under small strain condition. It presents a good agreement on the stress level, but still shows some mismatch near the elasto-viscoplastic transition.

The FG microstructure seems to show enhanced mechanical properties for forming processes compared to CG microstructure, with higher elongation at a given temperature and strain rate than the CG microstructure.

**Acknowledgment**

The authors would like to begin by thanking the IMT Mines Albi engineering school and the Occitanie Region for financial aid towards this study ("Program: 00090414 / Folder: 21012548"). We would also like to offer special thanks to all the researchers and technicians in the ICA laboratories for their technical help and scientific advice.

**References**

[1] Y. Combres, « Propriétés du titane et de ses alliages », *Métaux Alliages Non Ferr.*, août 2015. https://doi.org/10.51257/a-v1-m4780

[2] G. Gilles *et al.*, « Experimental characterization and elasto-plastic modeling of the quasi-static mechanical response of TA-6V at room temperature », *Int. J. Solids Struct.*, 48 (2011) 1277-1289. https://doi.org/10.1016/j.ijsolstr.2011.01.011

[3] N. Kotkunde, H. N. Krishnamurthy, P. Puranik, A. K. Gupta, et S. K. Singh, « Microstructure study and constitutive modeling of Ti-6Al-4V alloy at elevated temperatures », *Mater. Des.* 54 (2014) 96-103. https://doi.org/10.1016/j.matdes.2013.08.006

[4] L. Despax, « Étude des mécanismes de plasticité lors de la mise en forme de l'alliage de titane Ti-6Al-4V : influence de la microstructure initiale et des conditions de sollicitations thermomécaniques », *Thesis*, Ecole des Mines d'Albi-Carmaux, 2021. https://tel.archives-ouvertes.fr/tel-03579314
<figure><img src="image_2101.png" /><figcaption>Figure 7: Comparison between the experimental loading curves and the model at each temperature for the FG microstructure</figcaption></figure>

## Conclusion

This study investigated the mechanical behavior of two different microstructures of Ti-6Al-4V alloy at high temperature. The DIC method allowed to acquire information until very large strain (at least up to a true strain of 2). The results are in good agreement with the strain obtained using a contact method (extensometer).

An elasto-viscoplastic model has been identified using relaxation tests, under small strain condition. It presents a good agreement on the stress level, but still shows some mismatch near the elasto-viscoplastic transition.

The FG microstructure seems to show enhanced mechanical properties for forming processes compared to CG microstructure, with higher elongation at a given temperature and strain rate than the CG microstructure.

## Acknowledgment

The authors would like to begin by thanking the IMT Mines Albi engineering school and the Occitanie Region for financial aid towards this study ("Program: 00090414 / Folder: 21012548"). We would also like to offer special thanks to all the researchers and technicians in the ICA laboratories for their technical help and scientific advice.

## References

[1] Y. Combres, « Propriétés du titane et de ses alliages », *Métaux Alliages Non Ferr.*, août 2015. https://doi.org/10.51257/a-v1-m4780

[2] G. Gilles *et al.*, « Experimental characterization and elasto-plastic modeling of the quasi-static mechanical response of TA-6V at room temperature », *Int. J. Solids Struct.*, 48 (2011) 1277-1289. https://doi.org/10.1016/j.ijsolstr.2011.01.011

[3] N. Kotkunde, H. N. Krishnamurthy, P. Puranik, A. K. Gupta, et S. K. Singh, « Microstructure study and constitutive modeling of Ti-6Al-4V alloy at elevated temperatures », *Mater. Des.* 54 (2014) 96-103. https://doi.org/10.1016/j.matdes.2013.08.006

[4] L. Despax, « Étude des mécanismes de plasticité lors de la mise en forme de l'alliage de titane Ti-6Al-4V : influence de la microstructure initiale et des conditions de sollicitations thermomécaniques », *Thesis*, Ecole des Mines d'Albi-Carmaux, 2021. https://tel.archives-ouvertes.fr/tel-03579314
[5] R. Julien, V. Velay, V. Vidal, Y. Dahan, R. Forestier, et F. Rézaï-Aria, « Characterization and modeling of forged Ti-6Al-4V Titanium alloy with microstructural considerations during quenching process », *Int. J. Mech. Sci.*, 142-143 (2018) 456-467.
https://doi.org/10.1016/j.ijmecsci.2018.05.023

[6] Q. Sirvin, « Etude du comportement mécanique de tôles en alliage de titane et des paramètres procédé dans les opérations d'emboutissage à hautes températures », Thesis, Ecole des Mines d'Albi-Carmaux, 2018.

[7] H. J. Bong, D. H. Yoo, D. Kim, Y.-N. Kwon, et J. Lee, « Correlative Study on Plastic Response and Formability of Ti-6Al-4V Sheets under Hot Forming Conditions », *J. Manuf. Process.*, 58 (2020) 775-786. https://doi.org/10.1016/j.jmapro.2020.08.053

[8] G. R. Johnson et W. H. Cook, « A Constitutive Model And Data For Metals », *Proc. 7th Int. Symp. Ballist.*, p. 541-547, 1983.

[9] N. Kotkunde, A. D. Deole, A. K. Gupta, et S. K. Singh, « Comparative study of constitutive modeling for Ti-6Al-4V alloy at low strain rates and elevated temperatures », *Mater. Des.*, 55 (2014) 999-1005, doi: 10.1016/j.matdes.2013.10.089

[10] S. Gao, Y. Sang, Q. Li, Y. Sun, Y. Wu, et H. Wang, « Constitutive modeling and microstructure research on the deformation mechanism of Ti-6Al-4V alloy under hot forming condition », *J. Alloys Compounds*, 892 (2022) 162128.
https://doi.org/10.1016/j.jallcom.2021.162128

[11] P. Gao, M. Fu, M. Zhan, Z. Lei, et Y. Li, « Deformation behavior and microstructure evolution of titanium alloys with lamellar microstructure in hot working process: A review », *J. Mater. Sci. Technol.*, 39 (2020) 56-73. https://doi.org/10.1016/j.jmst.2019.07.052

[12] L. Despax, V. Vidal, D. Delagnes, M. Dehmas, H. Matsumoto, V. Velay, « Influence of strain rate and temperature on the deformation mechanisms of a fine-grained Ti-6Al-4V alloy », *Mater. Sci. Eng. A*, 790 (2020) 139718. https://doi.org/10.1016/j.msea.2020.139718

[13] P. Luong, R. Bonnaire, J. Périé, Q. Sirvin, et L. Penazzi, « Speckle pattern creation methods for two-dimensional digital image correlation strain measurements applied to mechanical tensile tests up to 700°C », *Strain*, 57 (2021) e12388. https://doi.org/10.1111/str.12388

[14] L. P. Luong, « Development of the two-dimension Digital Image Correlation method (2D-DIC) at high temperatures applicable for titanium alloys forming process », Thesis, Ecole des Mines d'Albi-Carmaux.

[15] E.-L. Odenberger, R. Pederson, et M. Oldenburg, « Finite element modeling and validation of springback and stress relaxation in the thermo-mechanical forming of thin Ti-6Al-4V sheets », *Int. J. Adv. Manuf. Technol.*, 104 (2019) 3439-3455. https://doi.org/10.1007/s00170-019-04071-9

[16] M. Rollin, « Study and Optimization of Cooling Conditions of Titanium Alloy Parts Manufactured by Superplastic Forming », Thèse, Ecole des Mines d'Albi-Carmaux, 2021.
[5] R. Julien, V. Velay, V. Vidal, Y. Dahan, R. Forestier, et F. Rézaï-Aria, « Characterization and modeling of forged Ti-6Al-4V Titanium alloy with microstructural considerations during quenching process », *Int. J. Mech. Sci.*, 142-143 (2018) 456-467.
https://doi.org/10.1016/j.ijmecsci.2018.05.023

[6] Q. Sirvin, « Etude du comportement mécanique de tôles en alliage de titane et des paramètres procédé dans les opérations d'emboutissage à hautes températures », Thesis, Ecole des Mines d'Albi-Carmaux, 2018.

[7] H. J. Bong, D. H. Yoo, D. Kim, Y.-N. Kwon, et J. Lee, « Correlative Study on Plastic Response and Formability of Ti-6Al-4V Sheets under Hot Forming Conditions », *J. Manuf. Process.*, 58 (2020) 775-786. https://doi.org/10.1016/j.jmapro.2020.08.053

[8] G. R. Johnson et W. H. Cook, « A Constitutive Model And Data For Metals », *Proc. 7th Int. Symp. Ballist.*, p. 541-547, 1983.

[9] N. Kotkunde, A. D. Deole, A. K. Gupta, et S. K. Singh, « Comparative study of constitutive modeling for Ti-6Al-4V alloy at low strain rates and elevated temperatures », *Mater. Des.*, 55 (2014) 999-1005, doi: 10.1016/j.matdes.2013.10.089

[10] S. Gao, Y. Sang, Q. Li, Y. Sun, Y. Wu, et H. Wang, « Constitutive modeling and microstructure research on the deformation mechanism of Ti-6Al-4V alloy under hot forming condition », *J. Alloys Compounds*, 892 (2022) 162128.
https://doi.org/10.1016/j.jallcom.2021.162128

[11] P. Gao, M. Fu, M. Zhan, Z. Lei, et Y. Li, « Deformation behavior and microstructure evolution of titanium alloys with lamellar microstructure in hot working process: A review », *J. Mater. Sci. Technol.*, 39 (2020) 56-73. https://doi.org/10.1016/j.jmst.2019.07.052

[12] L. Despax, V. Vidal, D. Delagnes, M. Dehmas, H. Matsumoto, V. Velay, « Influence of strain rate and temperature on the deformation mechanisms of a fine-grained Ti-6Al-4V alloy », *Mater. Sci. Eng. A*, 790 (2020) 139718. https://doi.org/10.1016/j.msea.2020.139718

[13] P. Luong, R. Bonnaire, J. Périé, Q. Sirvin, et L. Penazzi, « Speckle pattern creation methods for two-dimensional digital image correlation strain measurements applied to mechanical tensile tests up to 700°C », *Strain*, 57 (2021) e12388. https://doi.org/10.1111/str.12388

[14] L. P. Luong, « Development of the two-dimension Digital Image Correlation method (2D-DIC) at high temperatures applicable for titanium alloys forming process », Thesis, Ecole des Mines d'Albi-Carmaux.

[15] E.-L. Odenberger, R. Pederson, et M. Oldenburg, « Finite element modeling and validation of springback and stress relaxation in the thermo-mechanical forming of thin Ti-6Al-4V sheets », *Int. J. Adv. Manuf. Technol.*, 104 (2019) 3439-3455. https://doi.org/10.1007/s00170-019-04071-9

[16] M. Rollin, « Study and Optimization of Cooling Conditions of Titanium Alloy Parts Manufactured by Superplastic Forming », Thèse, Ecole des Mines d'Albi-Carmaux, 2021.
# Comparison between different methods to determine material constants of the ZK60 Mg alloy from hot bulge tests data

CUSANNO Angela<fnref n="1" />,ª*, PICCININI Antonio<fnref n="1" />,ᵇ, GUGLIELMI Pasquale<fnref n="1" />,ᶜ, SORGENTE Donato<fnref n="2" />,ᵈ, QIAO Jun<fnref n="3" />,ᵉ and PALUMBO Gianfranco<fnref n="1" />,*f

**Keywords:** Aluminium-Free Mg Alloy, ZK60, Superplastic Behavior, Bulge Tests, Inverse Analysis

**Abstract.** ZK60 is an Aluminium-free Magnesium alloy able to achieve high ductility in superplastic conditions. Those properties make the alloy suitable to be investigated for temporary prostheses production by means of the Superplastic Forming (SPF) process. To properly design the manufacturing process by numerical simulations, an accurate material constitutive model is needed. In the present work, bulge tests at 400°C were conducted on a ZK60 Magnesium sheet using different loading conditions, namely two different levels of constant pressure (CP) and jump pressure (JP) between the two pressure levels. The dome height evolution was acquired during each test and used to calibrate the material constants of the Backofen constitutive equation ($\sigma = C \cdot \dot{\varepsilon}^m$) using three different methodologies: (i) an analytical approach (ii) an inverse methodology based on a single CP test and (iii) an inverse methodology based on the JP test. The obtained constants were validated by subsequent numerical simulations. Comparing the numerical/experimental dome height curves, it was found that the sets of constants determined using the inverse methodology based on the JP test are able to describe the material superplastic behaviour over a wider span of loading conditions.

**Introduction**

Magnesium (Mg) based implants are gaining interest for temporary prostheses applications, since Mg mechanical properties are more similar to that of natural bones than the other metallic materials commonly used for implants, such as titanium alloys and stainless steels [1]. The ZK60 magnesium alloy can represent a viable alternative to the well known AZ31 for biomedical applications, since it is much safer for the human body due to the almost absence of Aluminium (Al < 0.05%), which is well-known for its toxicology related to neurodegenerative diseases (e.g., the Alzheimer disease) [2,3]. Furthermore, ZK60 showed good superplastic ductility in warm and hot conditions [4,5]. This ability to achieve high deformation make the alloy suitable to be processed by means of the Superplastic Forming (SPF) process [6], that uses gas pressure to form the sheet metal blank according to the die cavity shape. The resulting flexibility makes this technology suitable for the manufacturing of highly customised products, such as prostheses [7,8].

The adoption of the numerical approach allows to determine the gas pressure profile able to optimally exploit the material superplastic behaviour [8]. This requires a reliable FE model, based on an accurate constitutive material behavior. Tensile tests are commonly used to characterize materials superplastic behavior [9]. However, bulge tests can effectively reproduce strain conditions close to the ones occurring during the SPF process [10]. There are several types of bulge

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

<footnote n="1">¹DMMM - Politecnico di Bari, Via E. Orabona, 4, Bari, 70125, Italy</footnote>

<footnote n="2">² School of Engineering, Università degli Studi della Basilicata, Potenza, 85100, Italy</footnote>

<footnote n="3">³ School of Materials and Metallurgy, University of Science and Technology Liaoning, Qianshan Middlen Road 185, Anshan, 114051, P. R. China</footnote>

<footnote>ªangela.cusanno@poliba.it, ᵇantonio.piccininni@poliba.it, ᶜpasquale.guglielmi@poliba.it, ᵈdonato.sorgente@unibas.it, ᵉqiaojun@ustl.edu.cn, ᵇgianfranco.palumbo@poliba.it</footnote>
# Comparison between different methods to determine material constants of the ZK60 Mg alloy from hot bulge tests data

CUSANNO Angela<fnref n="1" />,ª*, PICCININI Antonio<fnref n="1" />,ᵇ, GUGLIELMI Pasquale<fnref n="1" />,ᶜ, SORGENTE Donato<fnref n="2" />,ᵈ, QIAO Jun<fnref n="3" />,ᵉ and PALUMBO Gianfranco<fnref n="1" />,ᵉ

**Keywords:** Aluminium-Free Mg Alloy, ZK60, Superplastic Behavior, Bulge Tests, Inverse Analysis

**Abstract.** ZK60 is an Aluminium-free Magnesium alloy able to achieve high ductility in superplastic conditions. Those properties make the alloy suitable to be investigated for temporary prostheses production by means of the Superplastic Forming (SPF) process. To properly design the manufacturing process by numerical simulations, an accurate material constitutive model is needed. In the present work, bulge tests at 400°C were conducted on a ZK60 Magnesium sheet using different loading conditions, namely two different levels of constant pressure (CP) and jump pressure (JP) between the two pressure levels. The dome height evolution was acquired during each test and used to calibrate the material constants of the Backofen constitutive equation ($\sigma = C \cdot \dot{\varepsilon}^m$) using three different methodologies: (i) an analytical approach (ii) an inverse methodology based on a single CP test and (iii) an inverse methodology based on the JP test. The obtained constants were validated by subsequent numerical simulations. Comparing the numerical/experimental dome height curves, it was found that the sets of constants determined using the inverse methodology based on the JP test are able to describe the material superplastic behaviour over a wider span of loading conditions.

**Introduction**

Magnesium (Mg) based implants are gaining interest for temporary prostheses applications, since Mg mechanical properties are more similar to that of natural bones than the other metallic materials commonly used for implants, such as titanium alloys and stainless steels [1]. The ZK60 magnesium alloy can represent a viable alternative to the well known AZ31 for biomedical applications, since it is much safer for the human body due to the almost absence of Aluminium (Al < 0.05%), which is well-known for its toxicology related to neurodegenerative diseases (e.g., the Alzheimer disease) [2,3]. Furthermore, ZK60 showed good superplastic ductility in warm and hot conditions [4,5]. This ability to achieve high deformation make the alloy suitable to be processed by means of the Superplastic Forming (SPF) process [6], that uses gas pressure to form the sheet metal blank according to the die cavity shape. The resulting flexibility makes this technology suitable for the manufacturing of highly customised products, such as prostheses [7,8].

The adoption of the numerical approach allows to determine the gas pressure profile able to optimally exploit the material superplastic behaviour [8]. This requires a reliable FE model, based on an accurate constitutive material behavior. Tensile tests are commonly used to characterize materials superplastic behavior [9]. However, bulge tests can effectively reproduce strain conditions close to the ones occurring during the SPF process [10]. There are several types of bulge

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

<footnote n="1">¹DMMM - Politecnico di Bari, Via E. Orabona, 4, Bari, 70125, Italy</footnote>

<footnote n="2">² School of Engineering, Università degli Studi della Basilicata, Potenza, 85100, Italy</footnote>

<footnote n="3">³ School of Materials and Metallurgy, University of Science and Technology Liaoning, Qianshan Middlen Road 185, Anshan, 114051, P. R. China</footnote>

<footnote>ªangela.cusanno@poliba.it, ᵇantonio.piccininni@poliba.it, ᶜpasquale.guglielmi@poliba.it, ᵈdonato.sorgente@unibas.it, ᵉqiaojun@ustl.edu.cn, ᵇgianfranco.palumbo@poliba.it</footnote>
tests available in literature, e.g., the Constant Pressure (CP) test [11], in which a constant gas pressure is applied and the Jump Pressure (JP) test [12], in which a stepped pressure profile changing between two values is used. The experimental data in terms of dome-height evolutions according to time and the thickness profiles at certain instants are then used for determining the material constants according to different constitutive equations, by using analytical methods or by inverse analysis [13].

In the present work the ZK60 magnesium alloy (Zr 0.85%, Zn 5.56%) was investigated using both CP and JP tests. The dome height evolutions according to time were acquired and used for calibrating the material constitutive equation. In particular, the behaviour of the alloy in superplastic conditions was modelled using the Backofen power law (Eq. 1):

$$ \sigma = C \cdot \dot{\varepsilon}^m \quad (1) $$

where $\sigma$ is the equivalent stress, $\dot{\varepsilon}$ is the equivalent strain rate, $C$ is the strength coefficient and $m$ is the strain rate sensitivity index. Three different methodologies were used to obtain the two constants of the model ($C$ and $m$): (i) by means of the analytical approach described in [11]; (ii) by means of an inverse methodology based on single CP tests and (iii) by means of an inverse methodology based on a single JP test [13]. The effectiveness of the approaches was finally evaluated by means of further numerical simulations using different loading conditions. In this way, the best set of material constants and the most effective methodology to determine them could be evaluated.

## Material and methodology

Material. Bulge tests were performed on circular Mg alloy ZK60 specimens with a 75 mm diameter and 1 mm thickness. The chemical composition of the investigated alloy is reported in Table 1.

<table><caption>Table 1 Chemical composition limits [Weight %] of the investigated Mg alloy ZK60</caption><thead><tr><th>Element</th><th>Al</th><th>Zr</th><th>Cu</th><th>Zn</th><th>Mn</th><th>Si</th><th>Ni</th><th>Fe</th><th>Mg</th></tr></thead><tbody><tr><td>Content</td><td>0.0487</td><td>0.85</td><td>0.041</td><td>5.558</td><td>0.01</td><td>0.046</td><td>0.0045</td><td>0.00399</td><td>Bal.</td></tr></tbody></table>

Experimental bulge tests. Figure 1 shows the experimental setup used for conducting hot bulge tests. The circular Mg specimen was clamped between the blankholder and the die. Then, Argon gas was used to deform the ZK60 Mg alloy blank at 400°C. Argon pressure was controlled by a proportional valve (0-10V). To reach the target temperature, a solid-state induction heating system managed with a Proportional-Integral-Derivative (PID) controller was used; temperature was constantly monitored by means of K-type wire thermocouples welded close to the test area. A BlankHolder Force (BHF) value of 12 kN was set to avoid gas leakages. A laser sensor was used to acquire the dome height evolution according to time.

Bulge tests were conducted under different loading conditions: (i) a constant pressure (CP) test, in which a constant gas pressure equal to 0.25 MPa was applied (indicated in the following as CP0.25), (ii) a CP test applying a constant gas pressure equal to 0.5 MPa (indicated in the following as CP0.5) and (iii) a Jump Pressure (JP) test in which a stepped pressure profile changing between the two values (0.5 MPa and 0.25 MPa) was used. Tests were carried out up to failure.
tests available in literature, e.g., the Constant Pressure (CP) test [11], in which a constant gas pressure is applied and the Jump Pressure (JP) test [12], in which a stepped pressure profile changing between two values is used. The experimental data in terms of dome-height evolutions according to time and the thickness profiles at certain instants are then used for determining the material constants according to different constitutive equations, by using analytical methods or by inverse analysis [13].

In the present work the ZK60 magnesium alloy (Zr 0.85%, Zn 5.56%) was investigated using both CP and JP tests. The dome height evolutions according to time were acquired and used for calibrating the material constitutive equation. In particular, the behaviour of the alloy in superplastic conditions was modelled using the Backofen power law (Eq. 1):

$$ \sigma = C \cdot \dot{\epsilon}^m \qquad (1) $$

where $\sigma$ is the equivalent stress, $\dot{\epsilon}$ is the equivalent strain rate, $C$ is the strength coefficient and $m$ is the strain rate sensitivity index. Three different methodologies were used to obtain the two constants of the model ($C$ and $m$): (i) by means of the analytical approach described in [11]; (ii) by means of an inverse methodology based on single CP tests and (iii) by means of an inverse methodology based on a single JP test [13]. The effectiveness of the approaches was finally evaluated by means of further numerical simulations using different loading conditions. In this way, the best set of material constants and the most effective methodology to determine them could be evaluated.

## Material and methodology

Material. Bulge tests were performed on circular Mg alloy ZK60 specimens with a 75 mm diameter and 1 mm thickness. The chemical composition of the investigated alloy is reported in Table 1.

<table><caption>Table 1 Chemical composition limits [Weight %] of the investigated Mg alloy ZK60</caption><thead><tr><th>Element</th><th>Al</th><th>Zr</th><th>Cu</th><th>Zn</th><th>Mn</th><th>Si</th><th>Ni</th><th>Fe</th><th>Mg</th></tr></thead><tbody><tr><td>Content</td><td>0.0487</td><td>0.85</td><td>0.041</td><td>5.558</td><td>0.01</td><td>0.046</td><td>0.0045</td><td>0.00399</td><td>Bal.</td></tr></tbody></table>

Experimental bulge tests. Figure 1 shows the experimental setup used for conducting hot bulge tests. The circular Mg specimen was clamped between the blankholder and the die. Then, Argon gas was used to deform the ZK60 Mg alloy blank at 400°C. Argon pressure was controlled by a proportional valve (0-10V). To reach the target temperature, a solid-state induction heating system managed with a Proportional-Integral-Derivative (PID) controller was used; temperature was constantly monitored by means of K-type wire thermocouples welded close to the test area. A BlankHolder Force (BHF) value of 12 kN was set to avoid gas leakages. A laser sensor was used to acquire the dome height evolution according to time.

Bulge tests were conducted under different loading conditions: (i) a constant pressure (CP) test, in which a constant gas pressure equal to 0.25 MPa was applied (indicated in the following as CP0.25), (ii) a CP test applying a constant gas pressure equal to 0.5 MPa (indicated in the following as CP0.5) and (iii) a Jump Pressure (JP) test in which a stepped pressure profile changing between the two values (0.5 MPa and 0.25 MPa) was used. Tests were carried out up to failure.
<figure><img src="image_2112.png" /><figcaption>Figure 1: Bulge tests experimental setup</figcaption></figure>

Analytical approach. The procedure proposed by Enikeev and Kruglov [11] was adopted for determining the material constants with an analytical approach. Data from dome-height evolutions of the experimental CP tests were used to determine the strain rate sensitivity index m as follows:

$$m = \frac{\ln \frac{p_1}{p_2}}{\ln \frac{t_2}{t_1}} \qquad (2)$$

Where $t_1$ and $t_2$ are the forming times corresponding to the $p_1$ and $p_2$ pressure values (0.25 and 0.5 MPa, respectively). Forming times refer to a dome height target value (the same for both tests), falling into the constant slope portion of the dome-height curve. The strength coefficient was calculated as follows:

$$C = \left( \frac{p_f R_0}{2S_0} \right) \cdot \left[ t_f / 2I_m(\alpha) \right]^m \qquad (3)$$

Where $t_f$ is the forming time (referring to the target dome height value) at a constant pressure $p_f$, $R_0$ and $S_0$ are the die radius (22.5 mm) and initial thickness (1.0 mm) of the undeformed blank, respectively. The term $I_m(\alpha)$ is an integral that can be calculated as follows:

$$I_m(\alpha) = \int_0^\alpha [(\sin^3 x/x^2)]^{1/m} \cdot (1/x - \cot x)dx \qquad (4)$$

Where $\alpha$ is half the angle subtended by a dome surface at its center of curvature.

FE model for bulge tests. The commercial software Abaqus/Standard was used to numerically reproduce the testing conditions; to reduce the computational cost, an axisymmetric FE model was used; it is shown in Figure 2. The die geometry was modelled by means of an analytical rigid surface, whereas the ZK60 blank was meshed with 280 SAX1 elements having a size equal to 0.1 mm. The blankholder was not considered and its clamping action was simulated by blocking the translational degree of freedom along the radial direction. Friction between the blank and the die was modelled according to the Penalty formulation for contact and Coulomb friction and its coefficient was set to 0.1.
<figure><img src="image_2103.png" /><figcaption>Figure 1: Bulge tests experimental setup</figcaption></figure>

Analytical approach. The procedure proposed by Enikeev and Kruglov [11] was adopted for determining the material constants with an analytical approach. Data from dome-height evolutions of the experimental CP tests were used to determine the strain rate sensitivity index m as follows:

$$m = \frac{\ln \frac{p_1}{p_2}}{\ln \frac{t_2}{t_1}} \qquad (2)$$

Where $t_1$ and $t_2$ are the forming times corresponding to the $p_1$ and $p_2$ pressure values (0.25 and 0.5 MPa, respectively). Forming times refer to a dome height target value (the same for both tests), falling into the constant slope portion of the dome-height curve. The strength coefficient was calculated as follows:

$$C = \left( \frac{p_f R_0}{2S_0} \right) \cdot \left[ t_f / 2I_m(\alpha) \right]^m \qquad (3)$$

Where $t_f$ is the forming time (referring to the target dome height value) at a constant pressure $p_f$, $R_0$ and $S_0$ are the die radius (22.5 mm) and initial thickness (1.0 mm) of the undeformed blank, respectively. The term $I_m(\alpha)$ is an integral that can be calculated as follows:

$$I_m(\alpha) = \int_0^\alpha [(\sin^3 x/x^2)]^{1/m} \cdot (1/x - \cot x) dx \qquad (4)$$

Where $\alpha$ is half the angle subtended by a dome surface at its center of curvature.

FE model for bulge tests. The commercial software Abaqus/Standard was used to numerically reproduce the testing conditions; to reduce the computational cost, an axisymmetric FE model was used; it is shown in Figure 2. The die geometry was modelled by means of an analytical rigid surface, whereas the ZK60 blank was meshed with 280 SAX1 elements having a size equal to 0.1 mm. The blankholder was not considered and its clamping action was simulated by blocking the translational degree of freedom along the radial direction. Friction between the blank and the die was modelled according to the Penalty formulation for contact and Coulomb friction and its coefficient was set to 0.1.
<figure><img src="image_2113.png" /><figcaption>Figure 2: FE model used for the numerical simulations of bulge tests</figcaption></figure>

The material behaviour was modelled using the Bailey-Norton power law ($\dot{\epsilon} = A \cdot \sigma^n$), in which the constants ($A, n$) are related to the one of the Backofen power law ($\sigma = C \cdot \dot{\epsilon}^m$) by the following equations:

$$n = 1/m \qquad (5)$$

$$A = (1/C)^{1/m} \qquad (6)$$

Pressure was applied to the upper surface of the blank. The dome-height evolution of the central point belonging to the bottom surface of the blank was recorded as output. The described numerical model was used both in the inverse analysis procedure and for validation purposes.

Inverse analysis approach. In the present work, material constants ($C, m$) in the Backofen constitutive equation (Eq. 1) were the unknown inputs to be determined by inverse analysis (IA). The numerical model described in the previous section was used to recreate the experimental working conditions; the experimental dome height curves were fitted by properly changing the value of the material constants ($C, m$) in the numerical simulations, minimizing the error between the numerical and the experimental data. The inverse analysis was driven by the MOGA-II genetic algorithm [14], that starts from an initial population created using the Sobol algorithm [15]. Additional details of the inverse analysis approach can be found in [13]. Each design represents a single numerical simulation characterized by specific values of the two constants ($C, m$); once each run was completed, the numerical evolution of the dome height was extracted by means of a python script file and written into an ASCII file. Subsequently, the numerical curve was compared with the experimental one and the sum of square residuals (SSR) calculated according to Eq. 7:

$$SSR = \sum_{i=0}^{N} (h_{i,num} - h_{i,exp})^2 \qquad (7)$$

Where $N$ is the total number of measured time instants, $h_{i,num}$ and $h_{i,exp}$ are the numerical and the experimental value of the dome height at the $i$-th instant of time, respectively. Results were analyzed in terms of the history charts, able to describe the evolution of the input variables throughout the created generations.

## Results and discussion

Experimental results. Figure 3 shows the experimental dome height evolutions of the ZK60 under the different pressure values, as well as the forming time and final height values.
<figure><img src="image_2104.png" /><figcaption>Figure 2: FE model used for the numerical simulations of bulge tests</figcaption></figure>

The material behaviour was modelled using the Bailey-Norton power law ($\dot{\epsilon} = A \cdot \sigma^n$), in which the constants ($A, n$) are related to the one of the Backofen power law ($\sigma = C \cdot \dot{\epsilon}^m$) by the following equations:

$$n = 1/m \qquad (5)$$

$$A = (1/C)^{1/m} \qquad (6)$$

Pressure was applied to the upper surface of the blank. The dome-height evolution of the central point belonging to the bottom surface of the blank was recorded as output. The described numerical model was used both in the inverse analysis procedure and for validation purposes.

Inverse analysis approach. In the present work, material constants ($C, m$) in the Backofen constitutive equation (Eq. 1) were the unknown inputs to be determined by inverse analysis (IA). The numerical model described in the previous section was used to recreate the experimental working conditions; the experimental dome height curves were fitted by properly changing the value of the material constants ($C, m$) in the numerical simulations, minimizing the error between the numerical and the experimental data. The inverse analysis was driven by the MOGA-II genetic algorithm [14], that starts from an initial population created using the Sobol algorithm [15]. Additional details of the inverse analysis approach can be found in [13]. Each design represents a single numerical simulation characterized by specific values of the two constants ($C, m$); once each run was completed, the numerical evolution of the dome height was extracted by means of a python script file and written into an ASCII file. Subsequently, the numerical curve was compared with the experimental one and the sum of square residuals (SSR) calculated according to Eq. 7:

$$SSR = \sum_{i=0}^{N} (h_{i,num} - h_{i,exp})^2 \qquad (7)$$

Where $N$ is the total number of measured time instants, $h_{i,num}$ and $h_{i,exp}$ are the numerical and the experimental value of the dome height at the $i$-th instant of time, respectively. Results were analyzed in terms of the history charts, able to describe the evolution of the input variables throughout the created generations.

## Results and discussion

Experimental results. Figure 3 shows the experimental dome height evolutions of the ZK60 under the different pressure values, as well as the forming time and final height values.
<figure><img src="image_2114.png" /><figcaption>Figure 3: (a) Experimental dome-height evolution of ZK60 during bulge tests at 400°C under different loading conditions (b) Loading conditions, forming time and final height of each test ID</figcaption></figure>

Material constants from the investigated approaches. Three different methodologies were used to obtain the two constants (C and m) of the Backofen material model, namely by means of (i) the analytical approach, (ii) inverse methodology based on the CP0.25 and CP0.5 tests and (iii) inverse methodology based on the JP test. Values of C and m resulting from the analytical approach were 2224 MPa * s^m and 0.63, respectively. For their determination a target dome height equal to 17 mm was considered, falling into the constant slope portion of both the CP 0.25 and the CP 0.5 dome-height curves as evidenced by Figure 3a.

Based on the results coming from the analytical approach, in the inverse analysis the strength coefficient C could initially vary between 100 and 4000 MPa * s^m, while the strain rate sensitivity index m between 0.3 and 0.7. The optimization loop was carried out considering an initial population of 50 individuals and 40 successive generations. Figure 4 shows the history charts resulting from the IA on the CP 0.25 test. Both C and m tended toward a stable and constant condition characterized by low SSR values after the 24th generation (i.e., after the ID 1200).

<figure><img src="image_2115.png" /><figcaption>Figure 4: History charts for the IA on the CP0.25 test: a) strength coefficient C and b) strain rate sensitivity index m.</figcaption></figure>

Figure 5 shows the history charts resulting from the IA on the CP 0.5 test. Unlike the previous case, C and m value tended toward a stable value almost at the end of the optimization loop, after the 32nd generation (i.e., after the ID 1600).
<figure><img src="image_2105.png" /><figcaption>Figure 3: (a) Experimental dome-height evolution of ZK60 during bulge tests at 400°C under different loading conditions (b) Loading conditions, forming time and final height of each test ID</figcaption></figure>

Material constants from the investigated approaches. Three different methodologies were used to obtain the two constants (C and m) of the Backofen material model, namely by means of (i) the analytical approach, (ii) inverse methodology based on the CP0.25 and CP0.5 tests and (iii) inverse methodology based on the JP test. Values of C and m resulting from the analytical approach were 2224 MPa * s^m and 0.63, respectively. For their determination a target dome height equal to 17 mm was considered, falling into the constant slope portion of both the CP 0.25 and the CP 0.5 dome-height curves as evidenced by Figure 3a.

Based on the results coming from the analytical approach, in the inverse analysis the strength coefficient C could initially vary between 100 and 4000 MPa * s^m, while the strain rate sensitivity index m between 0.3 and 0.7. The optimization loop was carried out considering an initial population of 50 individuals and 40 successive generations. Figure 4 shows the history charts resulting from the IA on the CP 0.25 test. Both C and m tended toward a stable and constant condition characterized by low SSR values after the 24th generation (i.e., after the ID 1200).

<figure><img src="image_2106.png" /><figcaption>Figure 4: History charts for the IA on the CP0.25 test: a) strength coefficient C and b) strain rate sensitivity index m.</figcaption></figure>

Figure 5 shows the history charts resulting from the IA on the CP 0.5 test. Unlike the previous case, C and m value tended toward a stable value almost at the end of the optimization loop, after the 32nd generation (i.e., after the ID 1600).
<figure><img src="image_2116.png" /><figcaption>Figure 5 History charts for the IA on the CP0.5 test: a) strength coefficient C and b) strain rate sensitivity index m.</figcaption></figure>

Figure 6 shows the history charts resulting from the IA on the JP test. After the 16th generation both the values of C and m tended already to a stable and constant condition, with minimum SSR values.

<figure><img src="image_2117.png" /><figcaption>Figure 6 History charts for the IA on the JP test: a) strength coefficient C and b) strain rate sensitivity index m.</figcaption></figure>

Figure 7a shows the comparison between the dome-height evolution of best design resulting from the IA and the experimental one (i.e., the target curve) for all the tested conditions. The best fitting of the target curve was obtained for the JP test, for which the input variables stably converged quicker. To compare results coming from different loading conditions (i.e., different forming times and, in turn, different numbers of measured time instants), the overall fitting was quantified by means of the Root Mean Square Error (RMSE), calculated as follows:

$$RMSE = \sqrt{\frac{\sum_{i=1}^{N} (h_{i,num} - h_{i,exp})^2}{N}} = \sqrt{\frac{SSR}{N}} \quad (8)$$

Where $h_{exp}$ and $h_{num}$ are the experimental and numerical dome height values and N is the total number of measured time instants. They are reported in Figure 7b. According to the results listed in Figure 7b, it appears that a lower applied pressure determined an increase of the calculated m value.
<figure><img src="image_2107.png" /><figcaption>Figure 5 History charts for the IA on the CP0.5 test: a) strength coefficient C and b) strain rate sensitivity index m.</figcaption></figure>

Figure 6 shows the history charts resulting from the IA on the JP test. After the 16th generation both the values of C and m tended already to a stable and constant condition, with minimum SSR values.

<figure><img src="image_2108.png" /><figcaption>Figure 6 History charts for the IA on the JP test: a) strength coefficient C and b) strain rate sensitivity index m.</figcaption></figure>

Figure 7a shows the comparison between the dome-height evolution of best design resulting from the IA and the experimental one (i.e., the target curve) for all the tested conditions. The best fitting of the target curve was obtained for the JP test, for which the input variables stably converged quicker. To compare results coming from different loading conditions (i.e., different forming times and, in turn, different numbers of measured time instants), the overall fitting was quantified by means of the Root Mean Square Error (RMSE), calculated as follows:

$$RMSE = \sqrt{\frac{\sum_{i=1}^{N} (h_{i,num} - h_{i,exp})^2}{N}} = \sqrt{\frac{SSR}{N}} \quad (8)$$

Where $h_{exp}$ and $h_{num}$ are the experimental and numerical dome height values and N is the total number of measured time instants. They are reported in Figure 7b. According to the results listed in Figure 7b, it appears that a lower applied pressure determined an increase of the calculated m value.
<figure><img src="image_2118.png" /><figcaption>Figure 7: (a) Comparison between the dome-height evolution of best design resulting from the IA and the experimental one for all the tested conditions; (b) Material constants from the optimization runs, sum of square residuals (SSR) and Root Mean Square Error (RMSE).</figcaption></figure>

<table><thead><tr><th>Test ID</th><th>C [MPa*s<sup>m</sup>]</th><th>m</th><th>SSR [mm<sup>2</sup>]</th><th>RMSE [mm]</th></tr></thead><tbody><tr><td>CP 0.25</td><td>794</td><td>0.524</td><td>354</td><td>0.6</td></tr><tr><td>CP 0.5</td><td>496</td><td>0.461</td><td>18.5</td><td>0.2</td></tr><tr><td>JP</td><td>888</td><td>0.532</td><td>3.7</td><td>0.0</td></tr></tbody></table>

Validation of the material constants. The previously determined constants able to describe the superplastic behavior of the ZK60 alloy were validated by means of further numerical simulations. In particular, the loading condition of the test CP0.25 was simulated implementing the material constants obtained (i) through the analytical (AN) approach (ii) by means of the inverse analysis carried out using data from the CP0.5 test and (iii) the set obtained through the inverse analysis on the JP test.

Figure 8 shows the numerical results in terms of dome height evolutions compared with the correspondent experimental ones (Figure 8a) as well as the RMSE between numerical and experimental dome-height curves (Figure 8b).

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

<figure><img src="image_2120.png" /><figcaption>Figure 8: Validation runs on CP0.25 (a) Numerical/experimental comparison of dome height evolution (b) Root Mean Square Error (RMSE) between numerical and experimental dome-height curves</figcaption></figure>

It can be noticed that implementing the material constants coming from the inverse analysis using data from the CP0.5 test (dotted curve in Figure 8a), the fitting was poor; in fact, as reported in Figure 8b, the RMSE value is close to 8 mm. This result suggests that when the constants are implemented to simulate a loading condition far from the one at which they had been determined using the inverse analysis, the fitting with the experimental dome height curve is not satisfactory. On the contrary, when implementing the material constants coming from the JP test (dashed
<figure><img src="image_2109.png" /><figcaption>Figure 7: (a) Comparison between the dome-height evolution of best design resulting from the IA and the experimental one for all the tested conditions; (b) Material constants from the optimization runs, sum of square residuals (SSR) and Root Mean Square Error (RMSE).</figcaption></figure>

<table><thead><tr><th>Test ID</th><th>C [MPa*s<sup>m</sup>]</th><th>m</th><th>SSR [mm<sup>2</sup>]</th><th>RMSE [mm]</th></tr></thead><tbody><tr><td>CP 0.25</td><td>794</td><td>0.524</td><td>354</td><td>0.6</td></tr><tr><td>CP 0.5</td><td>496</td><td>0.461</td><td>18.5</td><td>0.2</td></tr><tr><td>JP</td><td>888</td><td>0.532</td><td>3.7</td><td>0.0</td></tr></tbody></table>

Validation of the material constants. The previously determined constants able to describe the superplastic behavior of the ZK60 alloy were validated by means of further numerical simulations. In particular, the loading condition of the test CP0.25 was simulated implementing the material constants obtained (i) through the analytical (AN) approach (ii) by means of the inverse analysis carried out using data from the CP0.5 test and (iii) the set obtained through the inverse analysis on the JP test.

Figure 8 shows the numerical results in terms of dome height evolutions compared with the correspondent experimental ones (Figure 8a) as well as the RMSE between numerical and experimental dome-height curves (Figure 8b).

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

<figure><img src="image_2111.png" /><figcaption>Figure 8: Validation runs on CP0.25 (a) Numerical/experimental comparison of dome height evolution (b) Root Mean Square Error (RMSE) between numerical and experimental dome-height curves</figcaption></figure>

It can be noticed that implementing the material constants coming from the inverse analysis using data from the CP0.5 test (dotted curve in Figure 8a), the fitting was poor; in fact, as reported in Figure 8b, the RMSE value is close to 8 mm. This result suggests that when the constants are implemented to simulate a loading condition far from the one at which they had been determined using the inverse analysis, the fitting with the experimental dome height curve is not satisfactory. On the contrary, when implementing the material constants coming from the JP test (dashed
curve), the fitting to the experimental data was better, with an RMSE value equal to 1.3. However, in this case a much better fitting was achieved when implementing the material constants determined through the analytical approach (solid black curve), that are characterized by a RMSE value equal to 0.7, slightly higher than the minimum value of RMSE (0.6), obtained when implementing the material constants coming from the inverse analysis using data from the CP0.25 test.

The loading condition of the CP0.5 test was similarly investigated implementing the material constants coming from (i) the analytical (AN) approach (ii) the inverse analysis carried out using data from the CP0.25 test and (iii) the set obtained through the inverse analysis on the JP test. Unlike the previous case, for the CP0.5 validation runs the constants calculated through the analytical approach provided the worst fitting (RMSE=1.3). On the contrary, the JP test also provided in this case a good fitting. In fact, the resulting RMSE (0.3) is slightly higher than the minimum value of RMSE (0.2), obtained when implementing the material constants coming from the inverse analysis using data from the CP0.5 test.

<figure><img src="image_2121.png" /><figcaption>Figure 9: Validation runs on CP0.5 (a) Numerical/experimental comparison of dome height evolution (b) Root Mean Square Error (RMSE) between numerical and experimental dome-height curves</figcaption></figure>

Considering both the validation runs cases (applying a bulging pressure of 0.25 and 0.5 MPa), it can be stated that the set of constants obtained through the inverse analysis using data from the JP test (IA-JP) showed a good fitting for both the pressure values, meaning that the constants were able to effectively describe the material behavior over a wider span of loading conditions. On the contrary, the set of constants obtained through the analytical approach (AN) provided a good fitting only at lower pressure values. Eventually, constants obtained by means of the inverse analysis using data from the CP tests (IA-CP0.25 and IA-CP0.5) can be representative only if used to reproduce loading conditions close to the ones at which they had been determined using the inverse analysis, otherwise they won't provide good results.

## Conclusions

The present work investigates different methodologies for determining the material constants of the Backofen constitutive model able to describe the superplastic behaviour of the ZK60 alloy, starting from hot bulge tests data. Different sets of constants were obtained both through an analytical approach and by means of an inverse analysis approach. In this case, different sets of constants were obtained minimizing the difference between the numerical and the experimental dome height curves coming from both the constant pressure tests and the jump pressure test. The effectiveness of the approaches was evaluated by means of further simulations. Results showed
curve), the fitting to the experimental data was better, with an RMSE value equal to 1.3. However, in this case a much better fitting was achieved when implementing the material constants determined through the analytical approach (solid black curve), that are characterized by a RMSE value equal to 0.7, slightly higher than the minimum value of RMSE (0.6), obtained when implementing the material constants coming from the inverse analysis using data from the CP0.25 test.

The loading condition of the CP0.5 test was similarly investigated implementing the material constants coming from (i) the analytical (AN) approach (ii) the inverse analysis carried out using data from the CP0.25 test and (iii) the set obtained through the inverse analysis on the JP test. Unlike the previous case, for the CP0.5 validation runs the constants calculated through the analytical approach provided the worst fitting (RMSE=1.3). On the contrary, the JP test also provided in this case a good fitting. In fact, the resulting RMSE (0.3) is slightly higher than the minimum value of RMSE (0.2), obtained when implementing the material constants coming from the inverse analysis using data from the CP0.5 test.

<figure><img src="image_2112.png" /><figcaption>Figure 9: Validation runs on CP0.5 (a) Numerical/experimental comparison of dome height evolution (b) Root Mean Square Error (RMSE) between numerical and experimental dome-height curves</figcaption></figure>

Considering both the validation runs cases (applying a bulging pressure of 0.25 and 0.5 MPa), it can be stated that the set of constants obtained through the inverse analysis using data from the JP test (IA-JP) showed a good fitting for both the pressure values, meaning that the constants were able to effectively describe the material behavior over a wider span of loading conditions. On the contrary, the set of constants obtained through the analytical approach (AN) provided a good fitting only at lower pressure values. Eventually, constants obtained by means of the inverse analysis using data from the CP tests (IA-CP0.25 and IA-CP0.5) can be representative only if used to reproduce loading conditions close to the ones at which they had been determined using the inverse analysis, otherwise they won't provide good results.

## Conclusions

The present work investigates different methodologies for determining the material constants of the Backofen constitutive model able to describe the superplastic behaviour of the ZK60 alloy, starting from hot bulge tests data. Different sets of constants were obtained both through an analytical approach and by means of an inverse analysis approach. In this case, different sets of constants were obtained minimizing the difference between the numerical and the experimental dome height curves coming from both the constant pressure tests and the jump pressure test. The effectiveness of the approaches was evaluated by means of further simulations. Results showed
that the set of constants determined by means of the inverse analysis approach using data from the jump pressure test (C = 888 MPa * s^m and m = 0.532) satisfactorily described in a wider pressure ranges the superplastic behavior of the ZK60 alloy in terms of dome height evolution. On the contrary, material constants determined using data from the constant pressure tests provided excellent results for loading conditions close to the one from which they had been determined, but they lost their capability to describe the superplastic behavior for different loading conditions. The analytical approach provided a set of constants able to well describe only one loading condition. In addition, unlike the inverse analysis approaches, it requires two different tests with different loading conditions for determining only one set of material constants. Thus, it can be stated that the most effective methodology to adopt for determining reliable ZK60 material constants of the Backofen constitutive model is the one based on the inverse analysis approach using data from the jump pressure test.

Future works will be aimed at expanding the investigation to different temperatures and loading conditions, as well as evaluating the reliability of the determined material constants considering the thickness distribution of the bulged specimens.

## References

[1] B. Venkateswarlu, B.R. Sunil, R.S. Kumar, Magnesium based alloys and composites: Revolutionized biodegradable temporary implants and strategies to enhance their performance, *Materialia* 27 (2023) 101680. https://doi.org/10.1016/j.mtla.2023.101680

[2] J.R. Walton, Aluminum disruption of calcium homeostasis and signal transduction resembles change that occurs in aging and alzheimer's disease, *Journal of Alzheimer's Disease* 29 (2012) 255-273. https://doi.org/10.3233/JAD-2011-111712

[3] R.B. Heimann, Magnesium alloys for biomedical application: Advanced corrosion control through surface coating, *Surf Coat Technol* 405 (2021) 126521. https://doi.org/10.1016/j.surfcoat.2020.126521

[4] R. Lapovok, R. Cottam, P.F. Thomson, Y. Estrin, Extraordinary superplastic ductility of magnesium alloy ZK60, *J Mater Res* 20 (2005) 1375-1378. https://doi.org/10.1557/JMR.2005.0180

[5] W.J. Kim, M.J. Kim, J.Y. Wang, Superplastic behavior of a fine-grained ZK60 magnesium alloy processed by high-ratio differential speed rolling, *Mat. Sci. Eng: A* 527 (2009) 322-327. https://doi.org/10.1016/j.msea.2009.08.064

[6] Z. Savaedi, R. Motallebi, H. Mirzadeh, R. Mehdinavaz Aghdam, R. Mahmudi, Superplasticity of fine-grained magnesium alloys for biomedical applications: A comprehensive review, *Curr Opin Solid State Mater Sci* 27 (2023). https://doi.org/10.1016/j.cossms.2023.101058

[7] P. Guglielmi, A. Cusanno, I. Bagudanch, G. Centeno, I. Ferrer, M.L. Garcia-Romeu, G. Palumbo, Experimental and numerical analysis of innovative processes for producing a resorbable cheekbone prosthesis, *J Manuf Process* 70 (2021) 1-14. https://doi.org/10.1016/j.jmapro.2021.07.060

[8] D. Sorgente, G. Palumbo, A. Piccininni, P. Guglielmi, S.A. Aksenov, Investigation on the thickness distribution of highly customized titanium biomedical implants manufactured by superplastic forming, *CIRP J Manuf Sci Technol* 20 (2018) 29-35. https://doi.org/10.1016/j.cirpj.2017.09.004

[9] F. Nazeer, J. Long, Z. Yang, C. Li, Superplastic deformation behavior of Mg alloys: A-review, *Journal of Magnesium and Alloys* 10 (2022) 97-109. https://doi.org/10.1016/j.jma.2021.07.012
that the set of constants determined by means of the inverse analysis approach using data from the jump pressure test (C = 888 MPa * s^m and m = 0.532) satisfactorily described in a wider pressure ranges the superplastic behavior of the ZK60 alloy in terms of dome height evolution. On the contrary, material constants determined using data from the constant pressure tests provided excellent results for loading conditions close to the one from which they had been determined, but they lost their capability to describe the superplastic behavior for different loading conditions. The analytical approach provided a set of constants able to well describe only one loading condition. In addition, unlike the inverse analysis approaches, it requires two different tests with different loading conditions for determining only one set of material constants. Thus, it can be stated that the most effective methodology to adopt for determining reliable ZK60 material constants of the Backofen constitutive model is the one based on the inverse analysis approach using data from the jump pressure test.

Future works will be aimed at expanding the investigation to different temperatures and loading conditions, as well as evaluating the reliability of the determined material constants considering the thickness distribution of the bulged specimens.

## References

[1] B. Venkateswarlu, B.R. Sunil, R.S. Kumar, Magnesium based alloys and composites: Revolutionized biodegradable temporary implants and strategies to enhance their performance, *Materialia* 27 (2023) 101680. https://doi.org/10.1016/j.mtla.2023.101680

[2] J.R. Walton, Aluminum disruption of calcium homeostasis and signal transduction resembles change that occurs in aging and alzheimer's disease, *Journal of Alzheimer's Disease* 29 (2012) 255-273. https://doi.org/10.3233/JAD-2011-111712

[3] R.B. Heimann, Magnesium alloys for biomedical application: Advanced corrosion control through surface coating, *Surf Coat Technol* 405 (2021) 126521. https://doi.org/10.1016/j.surfcoat.2020.126521

[4] R. Lapovok, R. Cottam, P.F. Thomson, Y. Estrin, Extraordinary superplastic ductility of magnesium alloy ZK60, *J Mater Res* 20 (2005) 1375-1378. https://doi.org/10.1557/JMR.2005.0180

[5] W.J. Kim, M.J. Kim, J.Y. Wang, Superplastic behavior of a fine-grained ZK60 magnesium alloy processed by high-ratio differential speed rolling, *Mat. Sci. Eng: A* 527 (2009) 322-327. https://doi.org/10.1016/j.msea.2009.08.064

[6] Z. Savaedi, R. Motallebi, H. Mirzadeh, R. Mehdinavaz Aghdam, R. Mahmudi, Superplasticity of fine-grained magnesium alloys for biomedical applications: A comprehensive review, *Curr Opin Solid State Mater Sci* 27 (2023). https://doi.org/10.1016/j.cossms.2023.101058

[7] P. Guglielmi, A. Cusanno, I. Bagudanch, G. Centeno, I. Ferrer, M.L. Garcia-Romeu, G. Palumbo, Experimental and numerical analysis of innovative processes for producing a resorbable cheekbone prosthesis, *J Manuf Process* 70 (2021) 1-14. https://doi.org/10.1016/j.jmapro.2021.07.060

[8] D. Sorgente, G. Palumbo, A. Piccininni, P. Guglielmi, S.A. Aksenov, Investigation on the thickness distribution of highly customized titanium biomedical implants manufactured by superplastic forming, *CIRP J Manuf Sci Technol* 20 (2018) 29-35. https://doi.org/10.1016/j.cirpj.2017.09.004

[9] F. Nazeer, J. Long, Z. Yang, C. Li, Superplastic deformation behavior of Mg alloys: A-review, *Journal of Magnesium and Alloys* 10 (2022) 97-109. https://doi.org/10.1016/j.jma.2021.07.012
[10] S. Aksenov, D. Sorgente, Determination of biaxial stress–strain curves for superplastic materials by means of bulge forming tests at constant stress, CIRP J Manuf Sci Technol 31 (2020) 618–627. https://doi.org/10.1016/j.cirpj.2020.09.002

[11] F.U. Enikeev, A.A. Kruglov, An analysis of the superplastic forming of a thin circular diaphragm, Int. J. Mech. Sci 37 (1995) 473–483. https://doi.org/10.1016/0020-7403(94)00081-T

[12] S. Aksenov, D. Sorgente, Characterization of Stress-Strain Behavior of Superplastic Titanium Alloy by Free Bulging Tests with Pressure Jumps, Defect & Diffusion Forum 385 (2018) 443–448. https://doi.org/10.4028/www.scientific.net/DDF.385.443

[13] A. Piccininni, D. Sorgente, G. Palumbo, Genetic algorithm based inverse analysis for the superplastic characterization of a Ti-6Al-4V biomedical grade, Finite Elements in Analysis and Design 148 (2018) 27–37. https://doi.org/10.1016/j.finel.2018.06.002

[14] S. Poles, E. Rigoni, T. Robič, MOGA-II performance on noisy optimization problems, Proc. Int. Conf. Bioinspired Optim. Methods Their Appl (2004) 51–62.

[15] P. Bratley, B.L. Fox, Implementing Sobol's Quasirandom Sequence Generator, ACM Transactions on Mathematical Software 14 (1988) 88–100. https://doi.org/10.1145/42288.214372
[10] S. Aksenov, D. Sorgente, Determination of biaxial stress–strain curves for superplastic materials by means of bulge forming tests at constant stress, CIRP J Manuf Sci Technol 31 (2020) 618–627. https://doi.org/10.1016/j.cirpj.2020.09.002

[11] F.U. Enikeev, A.A. Kruglov, An analysis of the superplastic forming of a thin circular diaphragm, Int. J. Mech. Sci 37 (1995) 473–483. https://doi.org/10.1016/0020-7403(94)00081-T

[12] S. Aksenov, D. Sorgente, Characterization of Stress-Strain Behavior of Superplastic Titanium Alloy by Free Bulging Tests with Pressure Jumps, Defect & Diffusion Forum 385 (2018) 443–448. https://doi.org/10.4028/www.scientific.net/DDF.385.443

[13] A. Piccininni, D. Sorgente, G. Palumbo, Genetic algorithm based inverse analysis for the superplastic characterization of a Ti-6Al-4V biomedical grade, Finite Elements in Analysis and Design 148 (2018) 27–37. https://doi.org/10.1016/j.finel.2018.06.002

[14] S. Poles, E. Rigoni, T. Robič, MOGA-II performance on noisy optimization problems, Proc. Int. Conf. Bioinspired Optim. Methods Their Appl (2004) 51–62.

[15] P. Bratley, B.L. Fox, Implementing Sobol's Quasirandom Sequence Generator, ACM Transactions on Mathematical Software 14 (1988) 88–100. https://doi.org/10.1145/42288.214372
# Methodology for the automated spatial mapping of heterogeneous elastoplastic properties of welded joints

HAMILL Robert¹, a*, MAREK Aleksander¹, HARTE Allan² and
PIERRON Fabrice¹, 3, b

¹Faculty of Engineering and Physical Sciences, University of Southampton, UK

²UKAEA, Culham Science Centre, Oxfordshire, UK

³MatchID NV, Ghent, Belgium

a.r.hamill@soton.ac.uk, bfabrice.pierron@matchid.eu

**Keywords:** Welding, Virtual Fields Method, Identification, Full-Field Measurements, Heterogeneity, Materials Testing 2.0

**Abstract.** A methodology to map the heterogeneous elastoplastic mechanical properties of welded joints is presented. The approach and results of the proposed methodology are demonstrated using numerical data representing an off-axis butt-weld. This work extends the Virtual Fields Method (VFM) by introducing automated spatial parameterisation of the constitutive parameters. This extension enables the novel characterisation of welds with more complex geometries, loading conditions and dissimilar materials.

## Introduction

Inverse identification methods have been previously developed to take advantage of the rich data provided through full-field measurements. The VFM is one such method which uses full-field strain measurements to determine constitutive parameters, and benefits from computational efficiency compared to other methods as no resolution of the direct problem is required [1,2]. The VFM has previously been used to determine the constitutive parameters of welded joints with simple butt-joint geometry, however existing work requires a priori spatial parameterisation of constitutive parameters [3,4]. Automated spatial parameterisation removes the necessity of a priori knowledge, hence enabling more complex tests to be performed.

## Theory and key concepts

The Virtual Fields Method. Based on the principle of virtual work (PVW), the virtual fields method (VFM) seeks to identify constitutive parameters that satisfy the weak form of the stress equilibrium equation for a given strain field. Stress fields can be reconstructed from measured strain fields using a postulated constitutive model and an initial guess of the constitutive parameters for each. The static admissibility of these stress fields can then be checked using the PVW. The constitutive parameter map is iterated until global equilibrium is satisfied, at which point the parameter map is said to be identified.

Eq. 1 shows the simplified equation for the PVW assuming quasi-static loading of a specimen with negligible body forces. This equation expresses the local stress equilibrium equation in the weak form with the force boundary conditions accounted for. It states that the internal virtual work (IVW) and the external virtual work (EVW) must balance. $\sigma_{ij}$ is the internal stress within the volume, $V$. $T_i$ is the traction vector acting on the boundary surface, $S$.

$$ \frac{-\int_V \sigma_{ij} \epsilon_{ij}^* dV}{\text{IVW}} + \frac{\int_{S_f} T_i u_i^* dS}{\text{EVW}} = 0 \quad (1) $$

<figure><img src="image_2122.png" /></figure>
# Methodology for the automated spatial mapping of heterogeneous elastoplastic properties of welded joints

HAMILL Robert¹, a*, MAREK Aleksander¹, HARTE Allan² and
PIERRON Fabrice¹, 3, b

¹Faculty of Engineering and Physical Sciences, University of Southampton, UK

²UKAEA, Culham Science Centre, Oxfordshire, UK

³MatchID NV, Ghent, Belgium

a.r.hamill@soton.ac.uk, bfabrice.pierron@matchid.eu

**Keywords:** Welding, Virtual Fields Method, Identification, Full-Field Measurements, Heterogeneity, Materials Testing 2.0

**Abstract.** A methodology to map the heterogeneous elastoplastic mechanical properties of welded joints is presented. The approach and results of the proposed methodology are demonstrated using numerical data representing an off-axis butt-weld. This work extends the Virtual Fields Method (VFM) by introducing automated spatial parameterisation of the constitutive parameters. This extension enables the novel characterisation of welds with more complex geometries, loading conditions and dissimilar materials.

**Introduction**

Inverse identification methods have been previously developed to take advantage of the rich data provided through full-field measurements. The VFM is one such method which uses full-field strain measurements to determine constitutive parameters, and benefits from computational efficiency compared to other methods as no resolution of the direct problem is required [1,2]. The VFM has previously been used to determine the constitutive parameters of welded joints with simple butt-joint geometry, however existing work requires a priori spatial parameterisation of constitutive parameters [3,4]. Automated spatial parameterisation removes the necessity of a priori knowledge, hence enabling more complex tests to be performed.

**Theory and key concepts**

The Virtual Fields Method. Based on the principle of virtual work (PVW), the virtual fields method (VFM) seeks to identify constitutive parameters that satisfy the weak form of the stress equilibrium equation for a given strain field. Stress fields can be reconstructed from measured strain fields using a postulated constitutive model and an initial guess of the constitutive parameters for each. The static admissibility of these stress fields can then be checked using the PVW. The constitutive parameter map is iterated until global equilibrium is satisfied, at which point the parameter map is said to be identified.

Eq. 1 shows the simplified equation for the PVW assuming quasi-static loading of a specimen with negligible body forces. This equation expresses the local stress equilibrium equation in the weak form with the force boundary conditions accounted for. It states that the internal virtual work (IVW) and the external virtual work (EVW) must balance. $\sigma_{ij}$ is the internal stress within the volume, $V$. $T_i$ is the traction vector acting on the boundary surface, $S$.

$$ \frac{-\int_V \sigma_{ij} \epsilon_{ij}^* dV}{\text{IVW}} + \frac{\int_{S_f} T_i u_i^* dS}{\text{EVW}} = 0 \quad (1) $$

<figure><img src="image_2113.png" /></figure>
The term $u_i^*$ is a test function, known as the virtual displacement. This term does not require physical interpretation and merely acts as a spatial weighting function over the body. The spatial derivatives of the virtual displacement are known as the virtual strains, $\epsilon_{ij}^*$. The virtual strains are unrelated to the physical strains and again only act as spatial weighting functions.

Ensuring constant virtual displacements at the boundaries of any unknown traction distributions allows the term $u_i^*$ to be factored outside the integral. The tractions in the EVW term can then be simplified to only require the resultant force – a value provided by the load cell of the test machine.

In this work, the stress is assumed constant throughout the thickness and the strain data provided using digital image correlation (DIC) will take the form of discrete data points. The full-field data should have sufficient spatial resolution to approximate the surface integral as a discrete sum using the mid-point rule. Hence, the PVW can be further simplified to produce Eq. 2, where $F_i$ is the applied force vector measured by the load cell and $s^p$ is the surface area associated with each strain data point. From this equation, a cost function can be formed as the squared difference between external and internal virtual work, and minimized for the sought parameters.

$$ -\sum_{p=1}^{nPts} (\sigma_{ij}^p \epsilon_{ij}^{*p} s^p) + u_i^* F_i = 0 \quad (2) $$

Sensitivity based virtual fields (SBVFs, [5]). They are automatically generated virtual fields which adapt their spatial weighting according to each parameter's stress sensitivity in order to enhance parameter identifiability.

The equilibrium gap indicator (EGI, [6,7]). A particular formulation of virtual fields can be used to identify local discrepancies or 'equilibrium gaps' in the stress fields. Assuming kinematically admissible strain fields and correctly defined thicknesses, these equilibrium gaps (local regions of inadmissible stress states) arise due to incorrect constitutive model or parameters. As such, the EGI can be used as a metric to assess spatial regions in which equilibrium is not satisfied. Previously, this formulation has been used to detect damage and to assess heterogeneity in specimens of smoothly varying stiffness. In this work the EGI is used both to identify the location of discrepancies in the parameter map and their relative magnitude. To do this, an inspection window is rastered across the stress field with virtual fields defined to negate any external work contribution from outside of the window boundary. The EGI only assesses stress equilibrium within each window, hence it evaluates if the spatial distribution of stress satisfies static admissibility. However, it is ignorant of the magnitude of applied force and hence stress levels (i.e. all datapoints in the stress field may be wrong by some factor but the EGI will not detect any issue providing the distribution is correct).

The force reconstruction error (FRE). Like the EGI, the FRE is another virtual field formulation which checks for the equilibrium of any cross-section with the applied load.

Eq. 3 shows how the reconstructed force for a given slice, $F_r^{slice}$, can be computed from the average longitudinal stresses of that slice, $\sigma_{ii}$, and its cross-sectional area, $A^{slice}$. Here, index i denotes the longitudinal direction of the specimen, in which the resultant load $F_a$ is measured. The FRE value is expressed as in Eq. 4.

$$ F_r^{slice} = \int_A \sigma_{ii}^p dA = h \int_L \sigma_{ii}^p dL = A^{slice} * \overline{\sigma_{ii}} \quad (3) $$

$$ FRE = F_a - F_r^{slice} \quad (4) $$

Spatial parameterisation. In order to reconstruct the stress fields from the measured strain fields, we need to assign a value for each constitutive parameter at every point in our dataset. However, this spatial map of parameters is unknown and to be identified, so a parameterisation scheme is
The term $u_i^*$ is a test function, known as the virtual displacement. This term does not require physical interpretation and merely acts as a spatial weighting function over the body. The spatial derivatives of the virtual displacement are known as the virtual strains, $\epsilon_{ij}^*$. The virtual strains are unrelated to the physical strains and again only act as spatial weighting functions.

Ensuring constant virtual displacements at the boundaries of any unknown traction distributions allows the term $u_i^*$ to be factored outside the integral. The tractions in the EVW term can then be simplified to only require the resultant force – a value provided by the load cell of the test machine.

In this work, the stress is assumed constant throughout the thickness and the strain data provided using digital image correlation (DIC) will take the form of discrete data points. The full-field data should have sufficient spatial resolution to approximate the surface integral as a discrete sum using the mid-point rule. Hence, the PVW can be further simplified to produce Eq. 2, where $F_i$ is the applied force vector measured by the load cell and $s^p$ is the surface area associated with each strain data point. From this equation, a cost function can be formed as the squared difference between external and internal virtual work, and minimized for the sought parameters.

$$ - \sum_{p=1}^{nPts} (\sigma_{ij}^p \epsilon_{ij}^{*p} s^p) + u_i^* F_i = 0 \quad (2) $$

Sensitivity based virtual fields (SBVFs, [5]). They are automatically generated virtual fields which adapt their spatial weighting according to each parameter's stress sensitivity in order to enhance parameter identifiability.

The equilibrium gap indicator (EGI, [6,7]). A particular formulation of virtual fields can be used to identify local discrepancies or 'equilibrium gaps' in the stress fields. Assuming kinematically admissible strain fields and correctly defined thicknesses, these equilibrium gaps (local regions of inadmissible stress states) arise due to incorrect constitutive model or parameters. As such, the EGI can be used as a metric to assess spatial regions in which equilibrium is not satisfied. Previously, this formulation has been used to detect damage and to assess heterogeneity in specimens of smoothly varying stiffness. In this work the EGI is used both to identify the location of discrepancies in the parameter map and their relative magnitude. To do this, an inspection window is rastered across the stress field with virtual fields defined to negate any external work contribution from outside of the window boundary. The EGI only assesses stress equilibrium within each window, hence it evaluates if the spatial distribution of stress satisfies static admissibility. However, it is ignorant of the magnitude of applied force and hence stress levels (i.e. all datapoints in the stress field may be wrong by some factor but the EGI will not detect any issue providing the distribution is correct).

The force reconstruction error (FRE). Like the EGI, the FRE is another virtual field formulation which checks for the equilibrium of any cross-section with the applied load.

Eq. 3 shows how the reconstructed force for a given slice, $F_r^{slice}$, can be computed from the average longitudinal stresses of that slice, $\sigma_{ii}$, and its cross-sectional area, $A^{slice}$. Here, index i denotes the longitudinal direction of the specimen, in which the resultant load $F_a$ is measured. The FRE value is expressed as in Eq. 4.

$$ F_r^{slice} = \int_A \sigma_{ii}^p dA = h \int_L \sigma_{ii}^p dL = A^{slice} * \overline{\sigma_{ii}} \quad (3) $$

$$ FRE = F_a - F_r^{slice} \quad (4) $$

Spatial parameterisation. In order to reconstruct the stress fields from the measured strain fields, we need to assign a value for each constitutive parameter at every point in our dataset. However, this spatial map of parameters is unknown and to be identified, so a parameterisation scheme is
required to reduce the number of unknowns to an acceptable level. The parameterisation scheme will reconstruct a spatial map of values for each constitutive parameter from a number of variables (less than the number of datapoints). To date, this problem has been tackled by manually discretising the specimen into subregions *a priori* and then assuming each subregion is homogeneous (e.g. [3,4]). However, *a priori* information about the spatial distribution of properties is not always known, and automating spatial parameterisation would enable novel characterisation of specimens with more complex geometries and loading conditions.

A mesh is perhaps the most intuitive parameterisation approach. Each datapoint is assigned a value depending on its position within the element it lies in. Degrees of freedom define the nodal values and shape functions are used to specify the distribution of values throughout the element. The mesh can be refined by changing the number of elements or by changing the order of the shape functions. However, meshed methods inherently required a rigid structure and the number of unknowns increases rapidly with refinement.

Meshless approaches can provide greater flexibility, and radial basis functions (RBFs) are used in this work due to their flexibility and relatively low degree of freedom count. A radial function is a univariate function which returns a value dependent only on the distance from the kernel origin to the point being evaluated. Multiple basis functions are combined, using the principle of superposition, to approximate the desired interpolated solution. This is the method used here.

## Method

The key concepts outlined above can be combined in such a way to extend the VFM to include the automated spatial parameterisation of the constitutive parameters. This process is outlined in Fig. 1. The dashed box encloses the core VFM which, as previously described, identifies a map of constitutive parameters that best satisfies global equilibrium for a particular spatial parameterisation scheme. The EGI and FRE metrics can then be used to assess and update the parameterisation scheme. The particulars of how this is done could be modified, however the pseudo-code below outlines the approach used in this work. The core methodology will be presented through a worked example, in which synthetic data is used to mimic that obtained experimentally.

<table><caption>ALGORITHM 1: PSEUDO-CODE FOR HETEROGENEOUS IDENTIFICATION METHODOLOGY</caption><tr><td>1</td><td>prepare test data</td></tr><tr><td>2</td><td>prepare homogeneous baseline data</td></tr><tr><td>3</td><td>perform initial identification with mesh parameterisation and SBVF</td></tr><tr><td>4</td><td>initialise RBF parameterisation using initial identified results</td></tr><tr><td>5</td><td>perform identification with RBF parameterisation and local equilibrium metric VFs</td></tr><tr><td>6</td><td>while optimisation not converged do</td></tr><tr><td>7</td><td>update RBF parameterisation and proceed with identification</td></tr><tr><td>8</td><td>end while</td></tr></table>
required to reduce the number of unknowns to an acceptable level. The parameterisation scheme will reconstruct a spatial map of values for each constitutive parameter from a number of variables (less than the number of datapoints). To date, this problem has been tackled by manually discretising the specimen into subregions *a priori* and then assuming each subregion is homogeneous (e.g. [3,4]). However, *a priori* information about the spatial distribution of properties is not always known, and automating spatial parameterisation would enable novel characterisation of specimens with more complex geometries and loading conditions.

A mesh is perhaps the most intuitive parameterisation approach. Each datapoint is assigned a value depending on its position within the element it lies in. Degrees of freedom define the nodal values and shape functions are used to specify the distribution of values throughout the element. The mesh can be refined by changing the number of elements or by changing the order of the shape functions. However, meshed methods inherently required a rigid structure and the number of unknowns increases rapidly with refinement.

Meshless approaches can provide greater flexibility, and radial basis functions (RBFs) are used in this work due to their flexibility and relatively low degree of freedom count. A radial function is a univariate function which returns a value dependent only on the distance from the kernel origin to the point being evaluated. Multiple basis functions are combined, using the principle of superposition, to approximate the desired interpolated solution. This is the method used here.

## Method

The key concepts outlined above can be combined in such a way to extend the VFM to include the automated spatial parameterisation of the constitutive parameters. This process is outlined in Fig. 1. The dashed box encloses the core VFM which, as previously described, identifies a map of constitutive parameters that best satisfies global equilibrium for a particular spatial parameterisation scheme. The EGI and FRE metrics can then be used to assess and update the parameterisation scheme. The particulars of how this is done could be modified, however the pseudo-code below outlines the approach used in this work. The core methodology will be presented through a worked example, in which synthetic data is used to mimic that obtained experimentally.

<table><caption>ALGORITHM 1: PSEUDO-CODE FOR HETEROGENEOUS IDENTIFICATION METHODOLOGY</caption><tr><td>1</td><td>prepare test data</td></tr><tr><td>2</td><td>prepare homogeneous baseline data</td></tr><tr><td>3</td><td>perform initial identification with mesh parameterisation and SBVF</td></tr><tr><td>4</td><td>initialise RBF parameterisation using initial identified results</td></tr><tr><td>5</td><td>perform identification with RBF parameterisation and local equilibrium metric VFs</td></tr><tr><td>6</td><td>while optimisation not converged do</td></tr><tr><td>7</td><td>update RBF parameterisation and proceed with identification</td></tr><tr><td>8</td><td>end while</td></tr></table>
<figure><img src="image_2123.png" /><figcaption>Figure 1 – The virtual fields method with automated spatial parameterisation.</figcaption></figure>

Worked example

Test data generation. A notched, off-axis, butt-weld geometry is selected for the worked example (see Fig. 2). The weld region and heat affected zones have a width of 5.2 mm and 1.7 mm respectively. The off-axis angle is 60 degrees from the specimen longitudinal axis, and two offset notches penetrate the specimen 2 mm on either side. A linear hardening model is used together with Von Mises yield surface. Throughout the identification the elastic properties are assumed known, with a Young's modulus of 190 GPa and Poisson's ratio of 0.28. The yield stress will be identified as a heterogeneous parameter, and has a value of 360 MPa in the base metal and 420 MPa in the weld metal, with a linear variation throughout the heat affected zone. The hardening modulus, which will be identified as homogeneous, is assigned a value of 3700 MPa. Finite Element was used to simulate the test above, and synthetic image deformation [8] was used to realistically simulate DIC measurements. In this example, grey level noise is added directly to the images, leading to a longitudinal strain noise floor of 140 microstrain.

Initial identification. In order to get an initial indication of the heterogeneous parameter maps, a coarse, zero-order mesh is used, for which all points in an element share a single value. Hence, for the yield stress, the number of unknowns is equal to the number of elements. Initial mesh size is a user-choice, and here is set to be two elements in the vertical direction and six elements in the horizontal direction, giving a total of twelve unknowns governing the yield stress parameter map. The hardening modulus is assigned a single degree of freedom, giving thirteen unknowns in total. An initial guess of 350 MPa and 3000 MPa is assigned to the yield stress and hardening modulus respectively. The initial identification is performed using SBVFs and a gradient-descent optimizer. This results in a hardening modulus of 3296 MPa and a yield stress map as shown in the top row of Fig. 3.
<figure><img src="image_2114.png" /><figcaption>Figure 1 – The virtual fields method with automated spatial parameterisation.</figcaption></figure>

Worked example

Test data generation. A notched, off-axis, butt-weld geometry is selected for the worked example (see Fig. 2). The weld region and heat affected zones have a width of 5.2 mm and 1.7 mm respectively. The off-axis angle is 60 degrees from the specimen longitudinal axis, and two offset notches penetrate the specimen 2 mm on either side. A linear hardening model is used together with Von Mises yield surface. Throughout the identification the elastic properties are assumed known, with a Young's modulus of 190 GPa and Poisson's ratio of 0.28. The yield stress will be identified as a heterogeneous parameter, and has a value of 360 MPa in the base metal and 420 MPa in the weld metal, with a linear variation throughout the heat affected zone. The hardening modulus, which will be identified as homogeneous, is assigned a value of 3700 MPa. Finite Element was used to simulate the test above, and synthetic image deformation [8] was used to realistically simulate DIC measurements. In this example, grey level noise is added directly to the images, leading to a longitudinal strain noise floor of 140 microstrain.

Initial identification. In order to get an initial indication of the heterogeneous parameter maps, a coarse, zero-order mesh is used, for which all points in an element share a single value. Hence, for the yield stress, the number of unknowns is equal to the number of elements. Initial mesh size is a user-choice, and here is set to be two elements in the vertical direction and six elements in the horizontal direction, giving a total of twelve unknowns governing the yield stress parameter map. The hardening modulus is assigned a single degree of freedom, giving thirteen unknowns in total. An initial guess of 350 MPa and 3000 MPa is assigned to the yield stress and hardening modulus respectively. The initial identification is performed using SBVFs and a gradient-descent optimizer. This results in a hardening modulus of 3296 MPa and a yield stress map as shown in the top row of Fig. 3.
<figure><img src="image_2124.png" /><figcaption>Figure 2 –Notched, off-axis butt-weld geometry. Dimensions in mm.</figcaption></figure>

Initialize meshless parameterisation. The initial identification results can now be used to initialize the RBFs by fitting the RBFs to the identified parameter map (minimising the root mean squared error). The process of initially identifying using a coarse mesh and then fitting RBFs is not necessary, however, it quickly provides an approximate parameter map from which to start the global optimisation with RBFs.

Identification using RBF, EGI and FRE. Using the fitted RBF map as a starting point, identification is performed using a global optimiser to minimise a cost function composed of the normalised EGI and FRE terms. The EGI returns a 2D spatial array of values for each timestep. In this work, the root mean square of all spatial-temporal values is normalised by the baseline EGI obtained from a homogeneous case. The FRE returns a 1D spatial array, which can be extrapolated to a 2D array for each load step. Once again, the root mean square of all spatial-temporal values is normalised by the baseline FRE value. In this work, the pattern search algorithm is used, however, other optimisers such as genetic algorithms have also been used successfully.

Updating parameterisation and checking convergence. If there is no improvement after an iteration, a second attempt is made with the same number of kernels but a different starting point. If a second attempt has already been tried, the identification is ended. If there is an improvement in the cost function and the EGI and FRE metrics are below the baseline threshold, then the identification is ended, as no further improvement is feasible. Alternatively, if there is an improvement in the cost function, and the baseline threshold has not yet been reached, the parameterisation scheme is updated. To update the parameterisation, a new kernel is added to any heterogeneous parameters. This new kernel is positioned using the EGI and FRE, which highlight regions with the largest inconsistencies in the stress field.

## Results and discussion

The results are reported in Fig. 3 below. They clearly show how the process is able to automatically converge towards the solution without any *a priori* information on the distribution of the properties. This is a great strength of the technique compared to previously published work where such a complex distribution would not have been possible to identify without any *a priori* parameterisation, as in [3]. The method seems robust to realistic camera noise, and computational times are reasonable, about 4 hours on an Intel® Core™ i7-9700 CPU at 3 GHz with 32 Gb of RAM, with an unoptimized and unparallelised Matlab code.
<figure><img src="image_2115.png" /><figcaption>Figure 2 –Notched, off-axis butt-weld geometry. Dimensions in mm.</figcaption></figure>

Initialize meshless parameterisation. The initial identification results can now be used to initialize the RBFs by fitting the RBFs to the identified parameter map (minimising the root mean squared error). The process of initially identifying using a coarse mesh and then fitting RBFs is not necessary, however, it quickly provides an approximate parameter map from which to start the global optimisation with RBFs.

Identification using RBF, EGI and FRE. Using the fitted RBF map as a starting point, identification is performed using a global optimiser to minimise a cost function composed of the normalised EGI and FRE terms. The EGI returns a 2D spatial array of values for each timestep. In this work, the root mean square of all spatial-temporal values is normalised by the baseline EGI obtained from a homogeneous case. The FRE returns a 1D spatial array, which can be extrapolated to a 2D array for each load step. Once again, the root mean square of all spatial-temporal values is normalised by the baseline FRE value. In this work, the pattern search algorithm is used, however, other optimisers such as genetic algorithms have also been used successfully.

Updating parameterisation and checking convergence. If there is no improvement after an iteration, a second attempt is made with the same number of kernels but a different starting point. If a second attempt has already been tried, the identification is ended. If there is an improvement in the cost function and the EGI and FRE metrics are below the baseline threshold, then the identification is ended, as no further improvement is feasible. Alternatively, if there is an improvement in the cost function, and the baseline threshold has not yet been reached, the parameterisation scheme is updated. To update the parameterisation, a new kernel is added to any heterogeneous parameters. This new kernel is positioned using the EGI and FRE, which highlight regions with the largest inconsistencies in the stress field.

## Results and discussion

The results are reported in Fig. 3 below. They clearly show how the process is able to automatically converge towards the solution without any *a priori* information on the distribution of the properties. This is a great strength of the technique compared to previously published work where such a complex distribution would not have been possible to identify without any *a priori* parameterisation, as in [3]. The method seems robust to realistic camera noise, and computational times are reasonable, about 4 hours on an Intel® Core™ i7-9700 CPU at 3 GHz with 32 Gb of RAM, with an unoptimized and unparallelised Matlab code.
<figure><img src="image_2125.png" /><figcaption>Figure 3 - Results.</figcaption></figure>

## Conclusion and future work

The present results are first steps towards an operational tool. Improvements can include the use of bivariate RBFs (with different ‘lengths’ in two orthogonal directions) instead of univariate ones. There is a question about the possibility to also identify a hardening modulus map but considering the low sensitivity of hardening, this is an open question. Finally, experimental validation is currently being performed on a butt weld with notches and then, on a T-joint.

## References

[1] F. Pierron and M. Grédiac, The Virtual Fields Method, Springer, 2012. https://doi.org/10.1007/978-1-4614-1824-5

[2] J. Martins, A. Andrade-Campos, and S. Thuillier, “Comparison of inverse identification strategies for constitutive mechanical models using full-field measurements,” International Journal of Mechanical Sciences, 145 (2018) 330–345, 2018, https://doi.org/10.1016/j.ijmecsci.2018.07.013

[3] M. A. Sutton, J. H. Yan, S. Avril, F. Pierron, and S. M. Adeeb, “Identification of heterogeneous constitutive parameters in a welded specimen: Uniform stress and virtual fields methods for
<figure><img src="image_2116.png" /><figcaption>Figure 3 - Results.</figcaption></figure>

## Conclusion and future work

The present results are first steps towards an operational tool. Improvements can include the use of bivariate RBFs (with different ‘lengths’ in two orthogonal directions) instead of univariate ones. There is a question about the possibility to also identify a hardening modulus map but considering the low sensitivity of hardening, this is an open question. Finally, experimental validation is currently being performed on a butt weld with notches and then, on a T-joint.

## References

[1] F. Pierron and M. Grédiac, The Virtual Fields Method, Springer, 2012. https://doi.org/10.1007/978-1-4614-1824-5

[2] J. Martins, A. Andrade-Campos, and S. Thuillier, “Comparison of inverse identification strategies for constitutive mechanical models using full-field measurements,” International Journal of Mechanical Sciences, 145 (2018) 330–345, 2018, https://doi.org/10.1016/j.ijmecsci.2018.07.013

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

[4] G. L. Louedec, F. Pierron, M. A. Sutton, and A. P. Reynolds, “Identification of the local elasto-plastic behavior of FSW welds using the virtual fields method,” Experimental Mechanics, 55 (2013) 849–859.

[5] A. Marek, F.M. Davis and F. Pierron, “Sensitivity-based virtual fields for the non-linear virtual fields method”, Computational Mechanics, 60 (2017) 409-431, http://dx.doi.org/10.1007/s00466-017-1411-6

[6] C. Devivier, F. Pierron, and M. Wisnom, “Impact damage detection in composite plates using deflectometry and the Virtual Fields Method,” Composites Part A: Applied Science and Manufacturing, 48 (2013) 201–218, http://dx.doi.org/10.1016/j.compositesa.2013.01.011

[7] J. Considine, F. Pierron, K. Turner, P. Lava, and X. Tang, “Smoothly varying inplane stiffness heterogeneity evaluated under uniaxial tensile stress,” Strain, 53 (2017) e12237, http://dx.doi.org/10.1111/str.12237

[8] M. Rossi, P. Lava, F. Pierron, D. Debruyne, M. Sasso, “Effect of DIC spatial resolution, noise and interpolation error on identification results with the VFM”, Strain, vol. 51, n0. 3, pp. 206-222, 2015. http://dx.doi.org/10.1111/str.12134
material property estimation," Experimental Mechanics, 48 (2008) 451-464,
http://dx.doi.org/10.1007/s11340-008-9132-6

[4] G. L. Louedec, F. Pierron, M. A. Sutton, and A. P. Reynolds, “Identification of the local elasto-plastic behavior of FSW welds using the virtual fields method,” Experimental Mechanics, 55 (2013) 849–859.

[5] A. Marek, F.M. Davis and F. Pierron, “Sensitivity-based virtual fields for the non-linear virtual fields method”, Computational Mechanics, 60 (2017) 409-431, http://dx.doi.org/10.1007/s00466-017-1411-6

[6] C. Devivier, F. Pierron, and M. Wisnom, “Impact damage detection in composite plates using deflectometry and the Virtual Fields Method,” Composites Part A: Applied Science and Manufacturing, 48 (2013) 201–218, http://dx.doi.org/10.1016/j.compositesa.2013.01.011

[7] J. Considine, F. Pierron, K. Turner, P. Lava, and X. Tang, “Smoothly varying inplane stiffness heterogeneity evaluated under uniaxial tensile stress,” Strain, 53 (2017) e12237, http://dx.doi.org/10.1111/str.12237

[8] M. Rossi, P. Lava, F. Pierron, D. Debruyne, M. Sasso, “Effect of DIC spatial resolution, noise and interpolation error on identification results with the VFM”, Strain, vol. 51, n0. 3, pp. 206-222, 2015. http://dx.doi.org/10.1111/str.12134
# Investigations on the temperature and strain rate dependent behavior of a reinforced thermoplastic: application in hot incremental forming

FADDOULD Josephine¹,²,ᵃ, RAHME Pierre³,ᵇ, GUINES Dominique¹,ᶜ* and LEOTOING Lionel¹,ᵈ*

¹Univ Rennes, INSA Rennes, LGCGM (Laboratoire de Génie Civil et Génie Mécanique), Rennes, France

²Faculty of Engineering, Mechanical Engineering Department, Lebanese University, Lebanon

³Department of Industrial and Mechanical Engineering, Lebanese American University, Lebanon

**Keywords:** Reinforced Thermoplastic, Thermo-Viscoplasticity, FE Modeling, Mechanical Characterization

**Abstract.** This study focuses on characterizing the thermo-viscoplastic behavior of a glass fiber reinforced polypropylene through an advanced characterization technique. Based on a preliminary work addressing a mechanical characterization from conventional uniaxial tests, the proposed rheological model is calibrated using the database of in-plane equi-biaxial tensile tests on a dedicated cruciform specimen geometry from room temperature up to 140°C. The Finite Element Model Update (FEMU) method is employed for parameter calibration. A comparative analysis with the model calibrated from prior uniaxial tests is conducted through finite element simulations of a forming process. This investigation includes a parametric study to objectively assess the impact of the behavior law on forming force and final profile. This parametric study has direct implications for the simulation of an original forming process. The ability to identify the differences in forming outcomes based on different behavior laws enhances the predictive accuracy of the simulations. This ensures that the chosen behavior law aligns closely with the real behavior of the material during forming.

**Introduction**

Fiber reinforced thermoplastics, known for their good mechanical properties and recyclability, are extensively applied in engineering applications. Finite element (FE) simulations of the forming processes are essential to avoid trials and errors but require accurate representation of the material behavior under real conditions: multiaxial loading with large deformations, varying temperatures, and strain rates [1].

In-plane biaxial tensile tests using flat cruciform specimens are promising for achieving biaxial loading states. This frictionless technique allows direct application of loads along the two perpendicular arms of the specimen, offering versatility with multiple strain paths, from uniaxial to equi-biaxial stretching. Depending on the requirements and objectives of the mechanical characterization, the shape of the cruciform specimen can be modified and adjusted. This technique is commonly used to characterize metals. For instance, Liu et al. [2] proposed an optimized specimen design to investigate the hardening behavior of metallic sheets under conditions involving significant strains. Recently, some studies have explored the applicability of this technique to fiber reinforced polymers. Smits et al. [3] conducted a comparative study to investigate the influence of different specimen geometries on the biaxial failure strain of glass fiber reinforced epoxy composites. Their work aimed to identify a suitable specimen geometry for

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

<footnote>ᵃjosephine.faddoul@insa-rennes.fr, ᵇpierre.rahme@lau.edu.lb,</footnote>

<footnote>ᶜdominique.guines@insa-rennes.fr, ᵈlionel.leotoing@insa-rennes.fr</footnote>
# Investigations on the temperature and strain rate dependent behavior of a reinforced thermoplastic: application in hot incremental forming

FADDOULD Josephine¹,²,ᵃ, RAHME Pierre³,ᵇ, GUINES Dominique¹,ᶜ* and LEOTOING Lionel¹,ᵈ*

¹Univ Rennes, INSA Rennes, LGCGM (Laboratoire de Génie Civil et Génie Mécanique), Rennes, France

²Faculty of Engineering, Mechanical Engineering Department, Lebanese University, Lebanon

³Department of Industrial and Mechanical Engineering, Lebanese American University, Lebanon

ᵃjosephine.faddoul@insa-rennes.fr, ᵇpierre.rahme@lau.edu.lb,

ᶜdominique.guines@insa-rennes.fr, ᵈlionel.leotoing@insa-rennes.fr

**Keywords:** Reinforced Thermoplastic, Thermo-Viscoplasticity, FE Modeling, Mechanical Characterization

**Abstract.** This study focuses on characterizing the thermo-viscoplastic behavior of a glass fiber reinforced polypropylene through an advanced characterization technique. Based on a preliminary work addressing a mechanical characterization from conventional uniaxial tests, the proposed rheological model is calibrated using the database of in-plane equi-biaxial tensile tests on a dedicated cruciform specimen geometry from room temperature up to 140°C. The Finite Element Model Update (FEMU) method is employed for parameter calibration. A comparative analysis with the model calibrated from prior uniaxial tests is conducted through finite element simulations of a forming process. This investigation includes a parametric study to objectively assess the impact of the behavior law on forming force and final profile. This parametric study has direct implications for the simulation of an original forming process. The ability to identify the differences in forming outcomes based on different behavior laws enhances the predictive accuracy of the simulations. This ensures that the chosen behavior law aligns closely with the real behavior of the material during forming.

**Introduction**

Fiber reinforced thermoplastics, known for their good mechanical properties and recyclability, are extensively applied in engineering applications. Finite element (FE) simulations of the forming processes are essential to avoid trials and errors but require accurate representation of the material behavior under real conditions: multiaxial loading with large deformations, varying temperatures, and strain rates [1].

In-plane biaxial tensile tests using flat cruciform specimens are promising for achieving biaxial loading states. This frictionless technique allows direct application of loads along the two perpendicular arms of the specimen, offering versatility with multiple strain paths, from uniaxial to equi-biaxial stretching. Depending on the requirements and objectives of the mechanical characterization, the shape of the cruciform specimen can be modified and adjusted. This technique is commonly used to characterize metals. For instance, Liu et al. [2] proposed an optimized specimen design to investigate the hardening behavior of metallic sheets under conditions involving significant strains. Recently, some studies have explored the applicability of this technique to fiber reinforced polymers. Smits et al. [3] conducted a comparative study to investigate the influence of different specimen geometries on the biaxial failure strain of glass fiber reinforced epoxy composites. Their work aimed to identify a suitable specimen geometry for

<figure><img src="image_2117.png" /></figure>
conducting biaxial tests on fiber reinforced composites providing valuable insights into the testing methodology for these materials.

This study focuses on characterizing the thermo-viscoplastic behavior of a glass fiber reinforced polypropylene using in-plane biaxial tensile tests conducted on a dedicated cruciform specimen. The material parameters for the modified G'Sell and Jonas model proposed in a prior study [4] are identified relying on the biaxial experimental results. The calibrated model is then used to simulate a hot incremental sheet forming process. The forming outcomes are compared with simulations using material parameters identified from conventional uniaxial tensile tests.

## Material and experimental setup

The material used in this work is a 40% discontinuous long glass fiber reinforced polypropylene thermoplastic. Injected molded plates of 510×300×2 mm³ are supplied by SABIC (Stamax, 40YM240 [5]). The cruciform specimen proposed by Liu et al. [2] is considered in this study. It was proven that this shape can guarantee large equivalent plastic strains in the central zone under biaxial tensile loadings. The design and dimensions of the specimen are presented in Fig. 1. The samples are extracted by water-jet cutting. From a 2mm initial thickness, a thinning of the central zone is achieved by milling (down to a reduced thickness of 0.625 mm) to ensure failure in this zone.

<figure><img src="image_2127.png" /><figcaption>Figure 1: Geometry of cruciform specimen designed by Liu et al. [2]</figcaption></figure>

A specialized quasi-static/dynamic biaxial tensile device is used in this work to perform the in-plane biaxial tensile tests. Fig. 2 shows the experimental apparatus for a test performed at a high temperature.
conducting biaxial tests on fiber reinforced composites providing valuable insights into the testing methodology for these materials.

This study focuses on characterizing the thermo-viscoplastic behavior of a glass fiber reinforced polypropylene using in-plane biaxial tensile tests conducted on a dedicated cruciform specimen. The material parameters for the modified G'Sell and Jonas model proposed in a prior study [4] are identified relying on the biaxial experimental results. The calibrated model is then used to simulate a hot incremental sheet forming process. The forming outcomes are compared with simulations using material parameters identified from conventional uniaxial tensile tests.

## Material and experimental setup

The material used in this work is a 40% discontinuous long glass fiber reinforced polypropylene thermoplastic. Injected molded plates of 510×300×2 mm³ are supplied by SABIC (Stamax, 40YM240 [5]). The cruciform specimen proposed by Liu et al. [2] is considered in this study. It was proven that this shape can guarantee large equivalent plastic strains in the central zone under biaxial tensile loadings. The design and dimensions of the specimen are presented in Fig. 1. The samples are extracted by water-jet cutting. From a 2mm initial thickness, a thinning of the central zone is achieved by milling (down to a reduced thickness of 0.625 mm) to ensure failure in this zone.

<figure><img src="image_2118.png" /><figcaption>Figure 1: Geometry of cruciform specimen designed by Liu et al. [2]</figcaption></figure>

A specialized quasi-static/dynamic biaxial tensile device is used in this work to perform the in-plane biaxial tensile tests. Fig. 2 shows the experimental apparatus for a test performed at a high temperature.
<figure><img src="image_2128.png" /><figcaption>Figure 2: Experimental apparatus of the in-plane biaxial tensile test</figcaption></figure>

Two load sensors, which connect the grip system and the sliding bar, are placed to measure the forces along both perpendicular directions X and Y. For tests conducted at elevated temperatures, a heating system is added to the apparatus. It comprises an airflow generator and an insulated box. A high-speed camera (PHOTRON FASTCAM NOVA S9) is positioned along the central vertical axis of the biaxial bench to capture sequential images of the specimen during the test. The in-plane strain fields are then measured using the Digital Image Correlation DIC technique. The images are post-processed using the software “GOM Correlate”.

## Experimental results and discussion

The temporal evolution of the biaxial forces ($F_X$ and $F_Y$) as well as the in-plane principal $\varepsilon_1$, $\varepsilon_2$ and equivalent $\bar{\varepsilon}$ strains, are mainly concerned. Fig. 3 shows the force and strain curves resulting from the in-plane biaxial tensile test at room temperature (RT) and quasi-static velocity of 0.1mm/s. As shown in the figure, maximum forces of $F_{xf} = 1400$ N and $F_{yf} = 1100$ N are obtained at failure (the subscript f designates the value at failure). The difference between both curves is mainly attributed to the material anisotropy. The influence is also visible on the curves showing the temporal evolution of $\varepsilon_1$ and $\varepsilon_2$. The maximum level of strain reached just before failure is $\varepsilon_{1f} = 1\%$ and $\varepsilon_{2f} = 2\%$, respectively. $\bar{\varepsilon}$ is found to have a maximum value of $\bar{\varepsilon}_f = 3.2\%$ at failure.

<figure><img src="image_2129.png" /><figcaption>Figure 3: In-plane biaxial tensile test at RT and $V=0.1$mm/s</figcaption></figure>

Consistent patterns in the evolution of both forces and strains were observed across all temperature conditions, up to 140°C. It was noted that the maximum force decreases in approximately a linear trend with temperature increase, whereas the maximum equivalent strain increases exponentially. The values are summarized in Table 1. $F_f$ is the average of $F_{xf}$ and $F_{yf}$
<figure><img src="image_2119.png" /><figcaption>Figure 2: Experimental apparatus of the in-plane biaxial tensile test</figcaption></figure>

Two load sensors, which connect the grip system and the sliding bar, are placed to measure the forces along both perpendicular directions X and Y. For tests conducted at elevated temperatures, a heating system is added to the apparatus. It comprises an airflow generator and an insulated box. A high-speed camera (PHOTRON FASTCAM NOVA S9) is positioned along the central vertical axis of the biaxial bench to capture sequential images of the specimen during the test. The in-plane strain fields are then measured using the Digital Image Correlation DIC technique. The images are post-processed using the software “GOM Correlate”.

## Experimental results and discussion

The temporal evolution of the biaxial forces ($F_X$ and $F_Y$) as well as the in-plane principal $\varepsilon_1$, $\varepsilon_2$ and equivalent $\bar{\varepsilon}$ strains, are mainly concerned. Fig. 3 shows the force and strain curves resulting from the in-plane biaxial tensile test at room temperature (RT) and quasi-static velocity of 0.1mm/s. As shown in the figure, maximum forces of $F_{Xf} = 1400$ N and $F_{Yf} = 1100$ N are obtained at failure (the subscript f designates the value at failure). The difference between both curves is mainly attributed to the material anisotropy. The influence is also visible on the curves showing the temporal evolution of $\varepsilon_1$ and $\varepsilon_2$. The maximum level of strain reached just before failure is $\varepsilon_{1f} = 1\%$ and $\varepsilon_{2f} = 2\%$, respectively. $\bar{\varepsilon}$ is found to have a maximum value of $\bar{\varepsilon}_f = 3.2\%$ at failure.

<figure><img src="image_2120.png" /><figcaption>Figure 3: In-plane biaxial tensile test at RT and $V=0.1$mm/s</figcaption></figure>

Consistent patterns in the evolution of both forces and strains were observed across all temperature conditions, up to 140°C. It was noted that the maximum force decreases in approximately a linear trend with temperature increase, whereas the maximum equivalent strain increases exponentially. The values are summarized in Table 1. $F_f$ is the average of $F_{Xf}$ and $F_{Yf}$
measured along both directions just before failure. The temperature effect on the maximum equivalent strain is clearly pronounced at the highest temperatures. A 5°C gap causes a significant increase in the maximum equivalent strain level, i.e. $\bar{\varepsilon}_f$ increases by 50% when the temperature increases from 135 to 140°C.

<table><caption>Table 1: Temperature sensitivity of maximum biaxial force and equivalent strain</caption><thead><tr><th>T (°C)</th><th>RT</th><th>70</th><th>100</th><th>120</th><th>130</th><th>135</th><th>140</th></tr></thead><tbody><tr><th>ε̄<sub>f</sub>(%)</th><td>3.1</td><td>4.6</td><td>5.6</td><td>15.8</td><td>17</td><td>22</td><td>33</td></tr><tr><th>F<sub>f</sub>(N)</th><td>1140</td><td>833</td><td>654</td><td>364</td><td>314</td><td>200</td><td>185</td></tr></tbody></table>

### Thermo-viscoplastic rheological behavior model

The phenomenological constitutive law of G'Sell and Jonas (1979, 1983) [6,7] is selected to model the behavior of the studied material. This model is considered as an important contribution to understanding the plastic behavior of polymers. It describes the uniaxial stress-strain relation in dependence of strain and strain rate. The model is also capable of representing the hardening of the material after yielding. Several studies relied on this law to predict the stress-strain response of thermoplastics as well as fiber reinforced thermoplastics. For instance, in the study of Schossig et al. [8], the strain-rate-dependent behavior of glass fiber reinforced polypropylene and polybutene 1 was empirically examined and modeled using G'Sell and Jonas formulation.

In this work, the model is presented through one-dimensional elasto-viscoplastic formulation (Eq. 1). The modification of the constitutive equation was proposed in a previous work [4] for a better description of the observed experimental behavior. Additionally, the effect of the temperature is incorporated in the model by performing a temperature sensitivity analysis for the material parameters and expressing them as a function of temperature.

$$ \begin{cases} \bar{\sigma}(\bar{\varepsilon}, \dot{\bar{\varepsilon}}, T) = E(\dot{\bar{\varepsilon}}, T) \cdot \bar{\varepsilon} \\ \bar{\sigma}(\bar{\varepsilon}, \dot{\bar{\varepsilon}}, T) = \sigma_y(\dot{\bar{\varepsilon}}, T) + K(T) \cdot (1 - \exp(-W(T)\bar{\varepsilon}_p)) \cdot \dot{\bar{\varepsilon}}^m \end{cases} \quad \begin{array}{l} \text{for } \bar{\sigma} < \sigma_y \\ \text{for } \bar{\sigma} > \sigma_y \end{array} \quad (1) $$

$\bar{\sigma}$ is the flow stress, $\bar{\varepsilon}$ and $\bar{\varepsilon}_p$ are the total and plastic equivalent strains, respectively. $\dot{\bar{\varepsilon}}$ and $T$ refer to the equivalent strain rate and the temperature, respectively. $K$ (scaling factor), and $W$ are material parameters. $m$ represents the strain rate sensitivity. $E$ and $\sigma_y$ are the Young's modulus and the yield stress, respectively.

Based on the stress-strain curves determined from the uniaxial tensile tests [4], the formulations of $m$ and $E$ are found as stated in Eq. 2 and Eq. 3, respectively.

$$ m = \left. \frac{\partial \ln(\bar{\sigma}(\bar{\varepsilon}, \dot{\bar{\varepsilon}}, T) - \sigma_y)}{\partial \ln \dot{\bar{\varepsilon}}} \right|_{\bar{\varepsilon}_p, T} \quad (2) $$

$$ E(\dot{\bar{\varepsilon}}, T) = E_0(T) \cdot \dot{\bar{\varepsilon}}^m = E_1 \cdot \exp \left( A_E \left( \frac{1}{T} - \frac{1}{T_g} \right) \right) \cdot \dot{\bar{\varepsilon}}^m \quad (3) $$

Where $E_I$ and $A_E$ are material constants and $T_g$ is the glass transition temperature (-10°C).
measured along both directions just before failure. The temperature effect on the maximum equivalent strain is clearly pronounced at the highest temperatures. A 5°C gap causes a significant increase in the maximum equivalent strain level, i.e. $\bar{\varepsilon}_f$ increases by 50% when the temperature increases from 135 to 140°C.

<table><caption>Table 1: Temperature sensitivity of maximum biaxial force and equivalent strain</caption><thead><tr><th>T (°C)</th><th>RT</th><th>70</th><th>100</th><th>120</th><th>130</th><th>135</th><th>140</th></tr></thead><tbody><tr><th>ε̄<sub>f</sub>(%)</th><td>3.1</td><td>4.6</td><td>5.6</td><td>15.8</td><td>17</td><td>22</td><td>33</td></tr><tr><th>F<sub>f</sub>(N)</th><td>1140</td><td>833</td><td>654</td><td>364</td><td>314</td><td>200</td><td>185</td></tr></tbody></table>

### Thermo-viscoplastic rheological behavior model

The phenomenological constitutive law of G'Sell and Jonas (1979, 1983) [6,7] is selected to model the behavior of the studied material. This model is considered as an important contribution to understanding the plastic behavior of polymers. It describes the uniaxial stress-strain relation in dependence of strain and strain rate. The model is also capable of representing the hardening of the material after yielding. Several studies relied on this law to predict the stress-strain response of thermoplastics as well as fiber reinforced thermoplastics. For instance, in the study of Schossig et al. [8], the strain-rate-dependent behavior of glass fiber reinforced polypropylene and polybutene 1 was empirically examined and modeled using G'Sell and Jonas formulation.

In this work, the model is presented through one-dimensional elasto-viscoplastic formulation (Eq. 1). The modification of the constitutive equation was proposed in a previous work [4] for a better description of the observed experimental behavior. Additionally, the effect of the temperature is incorporated in the model by performing a temperature sensitivity analysis for the material parameters and expressing them as a function of temperature.

$$ egin{cases} \bar{\sigma}(\bar{\varepsilon}, \dot{\bar{\varepsilon}}, T) = E(\dot{\bar{\varepsilon}}, T) \cdot \bar{\varepsilon} \\ \bar{\sigma}(\bar{\varepsilon}, \dot{\bar{\varepsilon}}, T) = \sigma_y(\dot{\bar{\varepsilon}}, T) + K(T) \cdot (1 - \exp(-W(T)\bar{\varepsilon}_p)) \cdot \dot{\bar{\varepsilon}}^m \end{cases} \quad \begin{array}{l} \text{for } \bar{\sigma} < \sigma_y \\ \text{for } \bar{\sigma} > \sigma_y \end{array} \quad (1) $$

$\bar{\sigma}$ is the flow stress, $\bar{\varepsilon}$ and $\bar{\varepsilon}_p$ are the total and plastic equivalent strains, respectively. $\dot{\bar{\varepsilon}}$ and $T$ refer to the equivalent strain rate and the temperature, respectively. $K$ (scaling factor), and $W$ are material parameters. $m$ represents the strain rate sensitivity. $E$ and $\sigma_y$ are the Young's modulus and the yield stress, respectively.

Based on the stress-strain curves determined from the uniaxial tensile tests [4], the formulations of $m$ and $E$ are found as stated in Eq. 2 and Eq. 3, respectively.

$$ m = \left. \frac{\partial \ln(\bar{\sigma}(\bar{\varepsilon}, \dot{\bar{\varepsilon}}, T) - \sigma_y)}{\partial \ln \dot{\bar{\varepsilon}}} \right|_{\bar{\varepsilon}_p, T} \quad (2) $$

$$ E(\dot{\bar{\varepsilon}}, T) = E_0(T) \cdot \dot{\bar{\varepsilon}}^m = E_1 \cdot \exp \left( A_E \left( \frac{1}{T} - \frac{1}{T_g} \right) \right) \cdot \dot{\bar{\varepsilon}}^m \quad (3) $$

Where $E_I$ and $A_E$ are material constants and $T_g$ is the glass transition temperature (-10°C).
The proposed model was coupled with Hill 48 plasticity criterion to account for anisotropy [4]. The anisotropic parameters were calibrated based on uniaxial tensile tests for different material orientations.

## FE model of the in-plane biaxial test

The FE model, implemented in Abaqus software (implicit solver), is optimized for computational efficiency by modeling only one-quarter of the specimen, taking advantage of symmetry considerations (see Fig. 4).

<figure><img src="image_2130.png" /><figcaption>Figure 4: FE model of the cruciform specimen</figcaption></figure>

Symmetry conditions are imposed on sides 1 and 2, and the temporal evolution of forces ($F_x/2$ and $F_y/2$) obtained experimentally is applied along both perpendicular directions. Applying half of the experimental force reproduces the experimental conditions. The 3D mesh of the cruciform specimen is built with 8-nodes brick elements (C3D8R), with a refined element size of 0.1 mm and three elements arranged along the thickness direction in the central zone, focusing on the in-plane biaxial tensile state.

## Calibration of material parameters

For mechanical characterization through uniaxial tensile tests, stresses can be directly computed. In contrast, for in-plane biaxial tests with cruciform specimens, direct calculation of stresses and strains is not available. Consequently, based on full-field measurements, an inverse analysis procedure (Finite Element Model Update) becomes necessary for parameter identification. This method involves minimizing the gap between finite element simulation outputs and experimental observations by optimizing material parameters. The SIMPLEX algorithm is employed for the optimization. The temporal evolution of the experimental equivalent strain, from Digital Image Correlation DIC method, and averaged on the central zone of the specimen, is then compared with the simulated equivalent strain. The cost function $Q$ can be calculated (Eq. 4).

$$Q = \sqrt{\frac{\frac{1}{n}\sum_{i}^{n}(\bar{\varepsilon}_{exp}(t_i) - \bar{\varepsilon}_{sim}(t_i))^2}{\sum_{i}^{n}(\bar{\varepsilon}_{exp}(t_i))^2}} \quad (4)$$

Where the simulated and experimental equivalent strains $\bar{\varepsilon}_{sim}$ and $\bar{\varepsilon}_{exp}$ are calculated in the same manner using Eq. 5.

$$\bar{\varepsilon} = \frac{1}{3}(\sqrt{2[(\varepsilon_1 - \varepsilon_2)^2 + (\varepsilon_2 - \varepsilon_3)^2 + (\varepsilon_3 - \varepsilon_1)^2]}) \quad (5)$$

In this equation $\varepsilon_3$ is the out of plane principal strain calculated based on the volume conservation.

In this calibration strategy, the parameters $K$, $\sigma_y$, and $W$ are iteratively adjusted for a better fit with the experimental strain results. The optimized variables are listed in Table 2. The values of
The proposed model was coupled with Hill 48 plasticity criterion to account for anisotropy [4]. The anisotropic parameters were calibrated based on uniaxial tensile tests for different material orientations.

## FE model of the in-plane biaxial test

The FE model, implemented in Abaqus software (implicit solver), is optimized for computational efficiency by modeling only one-quarter of the specimen, taking advantage of symmetry considerations (see Fig. 4).

<figure><img src="image_2121.png" /><figcaption>Figure 4: FE model of the cruciform specimen</figcaption></figure>

Symmetry conditions are imposed on sides 1 and 2, and the temporal evolution of forces ($F_x/2$ and $F_y/2$) obtained experimentally is applied along both perpendicular directions. Applying half of the experimental force reproduces the experimental conditions. The 3D mesh of the cruciform specimen is built with 8-nodes brick elements (C3D8R), with a refined element size of 0.1 mm and three elements arranged along the thickness direction in the central zone, focusing on the in-plane biaxial tensile state.

## Calibration of material parameters

For mechanical characterization through uniaxial tensile tests, stresses can be directly computed. In contrast, for in-plane biaxial tests with cruciform specimens, direct calculation of stresses and strains is not available. Consequently, based on full-field measurements, an inverse analysis procedure (Finite Element Model Update) becomes necessary for parameter identification. This method involves minimizing the gap between finite element simulation outputs and experimental observations by optimizing material parameters. The SIMPLEX algorithm is employed for the optimization. The temporal evolution of the experimental equivalent strain, from Digital Image Correlation DIC method, and averaged on the central zone of the specimen, is then compared with the simulated equivalent strain. The cost function $Q$ can be calculated (Eq. 4).

$$Q = \sqrt{\frac{\frac{1}{n} \sum_{i}^{n} (\bar{\varepsilon}_{exp}(t_i) - \bar{\varepsilon}_{sim}(t_i))^2}{\sum_{i}^{n} (\bar{\varepsilon}_{exp}(t_i))^2}} \quad (4)$$

Where the simulated and experimental equivalent strains $\bar{\varepsilon}_{sim}$ and $\bar{\varepsilon}_{exp}$ are calculated in the same manner using Eq. 5.

$$\bar{\varepsilon} = \frac{1}{3}(\sqrt{2[(\varepsilon_1 - \varepsilon_2)^2 + (\varepsilon_2 - \varepsilon_3)^2 + (\varepsilon_3 - \varepsilon_1)^2]}) \quad (5)$$

In this equation $\varepsilon_3$ is the out of plane principal strain calculated based on the volume conservation.

In this calibration strategy, the parameters $K$, $\sigma_y$, and $W$ are iteratively adjusted for a better fit with the experimental strain results. The optimized variables are listed in Table 2. The values of
the cost function *Q* for each parameter set (determined at a specific temperature) show that at the end of the optimization procedure, the model's predictions closely match the experimental observations. For example, at the end of the optimization procedure at T=135°C, the prediction deviates from the experimental data by approximately 5% which is reasonably good.

A temperature sensitivity analysis shows a linear decrease of these parameters as temperature increases.

*E* and *m* are excluded from the optimization, based on their prior identification from uniaxial tensile tests [4]. The sensitivity of *E* to temperature and strain rate is expressed by the physical formulation presented in Eq. 3 with *E*<sub>l</sub>=12255MPa and *A*<sub>E</sub>=915. *m* is assumed to be constant (*m*=0.07) because of the negligible temperature effect on the strain rate sensitivity observed at high temperatures.

<table><caption>Table 2: Identified parameters based on the experimental equi-biaxial tensile tests</caption><thead><tr><td>T(°C)</td><td>K (MPa)</td><td>σ<sub>y</sub>(MPa)</td><td>W (l)</td><td>Q (%)</td></tr></thead><tbody><tr><td>20</td><td>78.9</td><td>25</td><td>134.4</td><td>5.2</td></tr><tr><td>70</td><td>36.2</td><td>29.3</td><td>117.9</td><td>7.6</td></tr><tr><td>100</td><td>30.3</td><td>23.9</td><td>72.8</td><td>6.6</td></tr><tr><td>120</td><td>24.3</td><td>10.4</td><td>44</td><td>2.6</td></tr><tr><td>130</td><td>22.5</td><td>8.9</td><td>28.75</td><td>6.7</td></tr><tr><td>135</td><td>19.9</td><td>5.6</td><td>11.2</td><td>5.2</td></tr><tr><td>140</td><td>17.5</td><td>6</td><td>7</td><td>2.5</td></tr></tbody></table>

FE simulations of the in-plane biaxial tensile tests are performed using the generalized formulation of the model. A well description of the temporal principal and equivalent strains evolution is observed for the tested temperature conditions. Fig. 5 presents the results obtained from the simulation at T=135°C and a quasi-static velocity. The correlation of the simulated in-plane principal strains ε<sub>sim1</sub> and ε<sub>sim2</sub> is relatively acceptable, even though these parameters are not directly targeted in the cost function.

<figure><img src="image_2131.png" /><figcaption>Figure 5: Comparison of experimental and numerical strain curves at T=135°C</figcaption></figure>

Finite element modeling of Incremental Sheet Forming

Based on the FE modeling of the Incremental Sheet Forming (ISF) process, a parametric study is proposed to investigate the impact of parameters like temperature and behavior law on forming forces and final shape of the part. The FE model is created using ABAQUS software (implicit solver), based on common assumptions from the literature [9]. These assumptions include the use
the cost function *Q* for each parameter set (determined at a specific temperature) show that at the end of the optimization procedure, the model's predictions closely match the experimental observations. For example, at the end of the optimization procedure at T=135°C, the prediction deviates from the experimental data by approximately 5% which is reasonably good.

A temperature sensitivity analysis shows a linear decrease of these parameters as temperature increases.

*E* and *m* are excluded from the optimization, based on their prior identification from uniaxial tensile tests [4]. The sensitivity of *E* to temperature and strain rate is expressed by the physical formulation presented in Eq. 3 with *E*<sub>l</sub>=12255MPa and *A*<sub>E</sub>=915. *m* is assumed to be constant (*m*=0.07) because of the negligible temperature effect on the strain rate sensitivity observed at high temperatures.

<table><caption>Table 2: Identified parameters based on the experimental equi-biaxial tensile tests</caption><thead><tr><td>T(°C)</td><td>K (MPa)</td><td>σ<sub>y</sub>(MPa)</td><td>W (l)</td><td>Q (%)</td></tr></thead><tbody><tr><td>20</td><td>78.9</td><td>25</td><td>134.4</td><td>5.2</td></tr><tr><td>70</td><td>36.2</td><td>29.3</td><td>117.9</td><td>7.6</td></tr><tr><td>100</td><td>30.3</td><td>23.9</td><td>72.8</td><td>6.6</td></tr><tr><td>120</td><td>24.3</td><td>10.4</td><td>44</td><td>2.6</td></tr><tr><td>130</td><td>22.5</td><td>8.9</td><td>28.75</td><td>6.7</td></tr><tr><td>135</td><td>19.9</td><td>5.6</td><td>11.2</td><td>5.2</td></tr><tr><td>140</td><td>17.5</td><td>6</td><td>7</td><td>2.5</td></tr></tbody></table>

FE simulations of the in-plane biaxial tensile tests are performed using the generalized formulation of the model. A well description of the temporal principal and equivalent strains evolution is observed for the tested temperature conditions. Fig. 5 presents the results obtained from the simulation at T=135°C and a quasi-static velocity. The correlation of the simulated in-plane principal strains ε<sub>sim1</sub> and ε<sub>sim2</sub> is relatively acceptable, even though these parameters are not directly targeted in the cost function.

<figure><img src="image_2122.png" /><figcaption>Figure 5: Comparison of experimental and numerical strain curves at T=135°C</figcaption></figure>

Finite element modeling of Incremental Sheet Forming

Based on the FE modeling of the Incremental Sheet Forming (ISF) process, a parametric study is proposed to investigate the impact of parameters like temperature and behavior law on forming forces and final shape of the part. The FE model is created using ABAQUS software (implicit solver), based on common assumptions from the literature [9]. These assumptions include the use
of shell elements and clamped sheet edges as boundary conditions. To simplify and reduce computation time, the study is conducted on a quarter of the model (see Fig. 6). The “stepped” trajectory shown in Fig. 7 is imposed on the tool. For each cycle, the tool maintains a maximum axial or radial depth increment of 1mm.

<figure><img src="image_2132.png" /><figcaption>Figure 6: Representation of the FE model</figcaption></figure>

<figure><img src="image_2133.png" /><figcaption>Figure 7: Tool trajectory with successive circular contours</figcaption></figure>

Fig. 8 illustrates the effect of the material behavior law on the tool axial force. At $T=120°C$, the difference between the simulated maximum axial force using the law identified from uniaxial tests and the one from in-plane biaxial tests is close to 11%. This difference increases to 21% at $T=130°C$, and remains stable for $T=140°C$. The results highlight the clear sensitivity of the forming force on the method to identify the parameters of the material law, which becomes more pronounced at higher temperatures. In contrast, the behavior law has a negligible impact on the simulated final shape for all the tested temperatures, as depicted in Fig. 9. For a better visualization of temperature's effect, the simulation results obtained using the behavior law calibrated from the biaxial database at different temperatures are shown in Fig. 10. Obviously, as temperature increases, the maximum axial force level at the end of forming decreases down to the lowest value of approximately 200N. On the other hand, the final formed shape of the plate is not affected by the temperature.

<figure><img src="image_2134.png" /><figcaption>Figure 8: Effect of behavior law on the tool reaction forces Fz for different test temperatures</figcaption></figure>
of shell elements and clamped sheet edges as boundary conditions. To simplify and reduce computation time, the study is conducted on a quarter of the model (see Fig. 6). The “stepped” trajectory shown in Fig. 7 is imposed on the tool. For each cycle, the tool maintains a maximum axial or radial depth increment of 1mm.

<figure><img src="image_2123.png" /><figcaption>Figure 6: Representation of the FE model</figcaption></figure>

<figure><img src="image_2124.png" /><figcaption>Figure 7: Tool trajectory with successive circular contours</figcaption></figure>

Fig. 8 illustrates the effect of the material behavior law on the tool axial force. At $T=120°C$, the difference between the simulated maximum axial force using the law identified from uniaxial tests and the one from in-plane biaxial tests is close to 11%. This difference increases to 21% at $T=130°C$, and remains stable for $T=140°C$. The results highlight the clear sensitivity of the forming force on the method to identify the parameters of the material law, which becomes more pronounced at higher temperatures. In contrast, the behavior law has a negligible impact on the simulated final shape for all the tested temperatures, as depicted in Fig. 9. For a better visualization of temperature's effect, the simulation results obtained using the behavior law calibrated from the biaxial database at different temperatures are shown in Fig. 10. Obviously, as temperature increases, the maximum axial force level at the end of forming decreases down to the lowest value of approximately 200N. On the other hand, the final formed shape of the plate is not affected by the temperature.

<figure><img src="image_2125.png" /><figcaption>Figure 8: Effect of behavior law on the tool reaction forces Fz for different test temperatures</figcaption></figure>
<figure><img src="image_2135.png" /><figcaption>Figure 9: Effect of behavior law on the profile for different test temperatures</figcaption></figure>

<figure><img src="image_2136.png" /><figcaption>Figure 10: Effect of temperature on the tool reaction force and final profile</figcaption></figure>

## Conclusion

In this work, the modified thermo-viscoplastic G'Sell and Jonas model is calibrated according to in-plane biaxial tensile tests on a dedicated cruciform specimen. The purpose of this new identification strategy is to propose a material characterization procedure based on a biaxial strain state close to the one encountered in forming processes. The effect of the identification strategy is evaluated thanks to the FE simulation of the hot incremental forming process of a truncated cone. A parametric study indicates that the simulated forming forces are significantly influenced by temperature and the choice of the rheological law i.e. calibration strategy for model parameters. However, these factors do not impact the simulated profile of the truncated cone. For future work, a heat assisted incremental forming process will be conducted experimentally, providing a basis for validation of the numerical predictive model. This will allow for a careful assessment of the degree of agreement between experimental and simulated results issued from the two different material characterization techniques.

## References

[1] T. Meinders, Emin Semih Perdahcıoğlu, M. van Riel, and H.H. Wisselink, “Numerical modeling of advanced materials,” *International Journal of Machine Tools & Manufacture*, 48 (2008) 485–498 https://doi.org/10.1016/j.ijmachtools.2007.08.005

[2] W. Liu, D. Guines, L. Leotoing, and E. Ragneau, “Identification of sheet metal hardening for large strains with an in-plane biaxial tensile test and a dedicated cross specimen,” *International Journal of Mechanical Sciences*, 101–102 (2015) 387–398. https://doi.org/10.1016/j.ijmecsci.2015.08.022

[3] A. Smits, D. Van Hemelrijck, T. P. Philippidis, and A. Cardon, “Design of a cruciform specimen for biaxial testing of fibre reinforced composite laminates”, *Composites science and technology*, 66 (2006) 964–975. https://doi.org/10.1016/j.compscitech.2005.08.011
<figure><img src="image_2126.png" /><figcaption>Figure 9: Effect of behavior law on the profile for different test temperatures</figcaption></figure>

<figure><img src="image_2127.png" /><figcaption>Figure 10: Effect of temperature on the tool reaction force and final profile</figcaption></figure>

## Conclusion

In this work, the modified thermo-viscoplastic G'Sell and Jonas model is calibrated according to in-plane biaxial tensile tests on a dedicated cruciform specimen. The purpose of this new identification strategy is to propose a material characterization procedure based on a biaxial strain state close to the one encountered in forming processes. The effect of the identification strategy is evaluated thanks to the FE simulation of the hot incremental forming process of a truncated cone. A parametric study indicates that the simulated forming forces are significantly influenced by temperature and the choice of the rheological law i.e. calibration strategy for model parameters. However, these factors do not impact the simulated profile of the truncated cone. For future work, a heat assisted incremental forming process will be conducted experimentally, providing a basis for validation of the numerical predictive model. This will allow for a careful assessment of the degree of agreement between experimental and simulated results issued from the two different material characterization techniques.

## References

[1] T. Meinders, Emin Semih Perdahcıoğlu, M. van Riel, and H.H. Wisselink, “Numerical modeling of advanced materials,” *International Journal of Machine Tools & Manufacture*, 48 (2008) 485–498 https://doi.org/10.1016/j.ijmachtools.2007.08.005

[2] W. Liu, D. Guines, L. Leotoing, and E. Ragneau, “Identification of sheet metal hardening for large strains with an in-plane biaxial tensile test and a dedicated cross specimen,” *International Journal of Mechanical Sciences*, 101–102 (2015) 387–398. https://doi.org/10.1016/j.ijmecsci.2015.08.022

[3] A. Smits, D. Van Hemelrijck, T. P. Philippidis, and A. Cardon, “Design of a cruciform specimen for biaxial testing of fibre reinforced composite laminates”, *Composites science and technology*, 66 (2006) 964–975. https://doi.org/10.1016/j.compscitech.2005.08.011
[4] J. Faddoul, P. Rahme, D. Guines, and L. Leotoing, "Thermo-visco mechanical behavior of glass fiber reinforced thermoplastic composite," *Journal of Composite Materials*, 57 (2023) 3741-3754. https://doi.org/10.1177/00219983231192147

[5] “StamaxTM,” SABIC, https://www.sabic.com/en/products/polymers/polypropylene-pp/stamax (accessed Feb. 4, 2024).

[6] C. G'Sell and J. J. Jonas, "Determination of the plastic behaviour of solid polymers at constant true strain rate," *Journal of Materials Science*, 14 (1979) 583-591. https://doi.org/10.1007/bf00772717

[7] C. G'sell, N. A. Aly-Helal, and J. J. Jonas, 'Effect of stress triaxiality on neck propagation during the tensile stretching of solid polymers', Journal of Materials Science, 18 (1983) 1731-1742. https://doi.org/10.1007/BF00542069

[8] M. Schoßig, C. Bierögel, W. Grellmann, and T. Mecklenburg, 'Mechanical behavior of glass-fiber reinforced thermoplastic materials under high strain rates', *Polymer testing*, 27 (2008) 893-900. https://doi.org/10.1016/j.polymertesting.2008.07.006.

[9] J. Belchior, D. Guines, L. Leotoing, and E. Ragneau, 'Force prediction for correction of robot tool path in single point incremental forming', *Key Engineering Materials*, 554 (2013) 1282-1289. https://doi.org/10.4028/www.scientific.net/kem.554-557.1282
[4] J. Faddoul, P. Rahme, D. Guines, and L. Leotoing, “Thermo-visco mechanical behavior of glass fiber reinforced thermoplastic composite,” *Journal of Composite Materials*, 57 (2023) 3741–3754. https://doi.org/10.1177/00219983231192147

[5] “StamaxTM,” SABIC, https://www.sabic.com/en/products/polymers/polypropylene-pp/stamax (accessed Feb. 4, 2024).

[6] C. G’Sell and J. J. Jonas, “Determination of the plastic behaviour of solid polymers at constant true strain rate,” *Journal of Materials Science*, 14 (1979) 583–591. https://doi.org/10.1007/bf00772717

[7] C. G’sell, N. A. Aly-Helal, and J. J. Jonas, ‘Effect of stress triaxiality on neck propagation during the tensile stretching of solid polymers’, Journal of Materials Science, 18 (1983) 1731–1742. https://doi.org/10.1007/BF00542069

[8] M. Schoßig, C. Bierögel, W. Grellmann, and T. Mecklenburg, ‘Mechanical behavior of glass-fiber reinforced thermoplastic materials under high strain rates’, *Polymer testing*, 27 (2008) 893–900. https://doi.org/10.1016/j.polymertesting.2008.07.006.

[9] J. Belchior, D. Guines, L. Leotoing, and E. Ragneau, ‘Force prediction for correction of robot tool path in single point incremental forming’, *Key Engineering Materials*, 554 (2013) 1282–1289. https://doi.org/10.4028/www.scientific.net/kem.554-557.1282
New and advanced numerical strategies
for material forming

<figure><img src="image_2137.png" /></figure>
New and advanced numerical strategies
for material forming

<figure><img src="image_2128.png" /></figure>
# A strategy for minimizing the computational time of simulations involving near-surface embossing of sheet metal materials

HEINZELMANN Pascal¹,ᵃ*, GÖRZ Marcel¹,ᵇ, RIEDMÜLLER Kim Rouven¹,ᶜ and LIEWALD Mathias¹, ᵈ

¹Institute for Metal Forming Technology, University of Stuttgart,

Holzgartenstraße 17, 70174 Stuttgart, Germany

ᵃpascal.heinzelmann; ᵇmarcel.goerz; ᶜkim.riedmueller; ᵈmathias.liewald}@ifu.uni-stuttgart.de

**Keywords:** Embossing, Dual-Phase Steel, Simulation, Mapping

**Abstract.** By near-surface embossing, work hardening can be introduced into sheet metal blanks, thus increasing the material's yield strength. Hence, this property-modifying process can be used to improve the lightweight potential of a component on the one hand and its crash performance on the other. For the early design of such embossing or forming processes and to predict the modification of the material properties, FE modelling is currently used. Here, the embossing of the near-surface structures is usually simulated first, followed by the forming process. However, existing simulation methods are time-consuming and computationally intensive, due to the simulation of the whole embossing process. For this purpose, this paper presents two approaches for improving the simulation time of embossing and forming processes. Compared to the conventional holistic embossing and forming simulation, these approaches are based on a new mapping strategy and experimental data. In the mapping variant, solid elements' stress and strain distribution is mapped to shell elements, thus including the embossing history for the forming simulation. For the experimental-based variant, the yield curves of embossed tensile tests are determined and used to approximate the property changes in the simulation. Both concepts are compared with each other and the conventional simulation method in terms of accuracy and calculation time. This paper finally shows that simulation can be performed faster with the new approaches than with the current sequential modelled workflow.

**Introduction**

Human-induced climate change is caused primarily by the emission of greenhouse gases such as CO₂, which is largely produced by motorised traffic. Therefore, both the US, the EU, and other countries have made commitments to significantly reduce these emissions, particularly from road transport [1]. In addition to filtering and storing such gases as they are emitted, another promising measure being driven by the automotive industry is the general reduction of fossil fuel consumption. Here, lightweight design approaches such as the use of lightweight aluminium alloys or tailoring mechanical material properties in critical areas of sheet metal components provide a promising contributions [2].

The selective modification of material properties in highly loaded areas of sheet metal components allows comparatively thin sheet metal materials to be used in manufacturing and thus to reduce weight of such components. In particular, the selective application of near-surface embossing can induce strain hardening in relevant areas of the component, which can lead to an improvement in mechanical properties [3]. Initial tests on this topic were carried out by Namoco et al. in 2006 on work-hardened aluminium alloys A6061-T4 and A5052-H34 [4]. Further, Abel et al. demonstrated the improved formability of embossed deep-drawn components made of aluminium sheets [5]. Studies with DP steels conducted by Walzer et al. revealed an increase in the tensile strength of the dual-phase (DP) steel DP600 due to near-surface embossing [6].

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

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 license. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Published under license by Materials Research Forum LLC.
# A strategy for minimizing the computational time of simulations involving near-surface embossing of sheet metal materials

HEINZELMANN Pascal¹,ᵃ*, GÖRZ Marcel¹,ᵇ, RIEDMÜLLER Kim Rouven¹,ᶜ and LIEWALD Mathias¹, ᵈ

¹Institute for Metal Forming Technology, University of Stuttgart,
Holzgartenstraße 17, 70174 Stuttgart, Germany

ᵃpascal.heinzelmann; ᵇmarcel.goerz; ᶜkim.riedmueller; ᵈmathias.liewald}@ifu.uni-stuttgart.de

**Keywords:** Embossing, Dual-Phase Steel, Simulation, Mapping

**Abstract.** By near-surface embossing, work hardening can be introduced into sheet metal blanks, thus increasing the material's yield strength. Hence, this property-modifying process can be used to improve the lightweight potential of a component on the one hand and its crash performance on the other. For the early design of such embossing or forming processes and to predict the modification of the material properties, FE modelling is currently used. Here, the embossing of the near-surface structures is usually simulated first, followed by the forming process. However, existing simulation methods are time-consuming and computationally intensive, due to the simulation of the whole embossing process. For this purpose, this paper presents two approaches for improving the simulation time of embossing and forming processes. Compared to the conventional holistic embossing and forming simulation, these approaches are based on a new mapping strategy and experimental data. In the mapping variant, solid elements' stress and strain distribution is mapped to shell elements, thus including the embossing history for the forming simulation. For the experimental-based variant, the yield curves of embossed tensile tests are determined and used to approximate the property changes in the simulation. Both concepts are compared with each other and the conventional simulation method in terms of accuracy and calculation time. This paper finally shows that simulation can be performed faster with the new approaches than with the current sequential modelled workflow.

**Introduction**

Human-induced climate change is caused primarily by the emission of greenhouse gases such as CO₂, which is largely produced by motorised traffic. Therefore, both the US, the EU, and other countries have made commitments to significantly reduce these emissions, particularly from road transport [1]. In addition to filtering and storing such gases as they are emitted, another promising measure being driven by the automotive industry is the general reduction of fossil fuel consumption. Here, lightweight design approaches such as the use of lightweight aluminium alloys or tailoring mechanical material properties in critical areas of sheet metal components provide a promising contributions [2].

The selective modification of material properties in highly loaded areas of sheet metal components allows comparatively thin sheet metal materials to be used in manufacturing and thus to reduce weight of such components. In particular, the selective application of near-surface embossing can induce strain hardening in relevant areas of the component, which can lead to an improvement in mechanical properties [3]. Initial tests on this topic were carried out by Namoco et al. in 2006 on work-hardened aluminium alloys A6061-T4 and A5052-H34 [4]. Further, Abel et al. demonstrated the improved formability of embossed deep-drawn components made of aluminium sheets [5]. Studies with DP steels conducted by Walzer et al. revealed an increase in the tensile strength of the dual-phase (DP) steel DP600 due to near-surface embossing [6].

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

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 license. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Published under license by Materials Research Forum LLC.
Additionally, three-point bending tests with embossed DP500 of Briesenick et al. showed an increase in the maximum bending force [7]. Complementary to the work of Walzer and Briesenick, tensile tests and three-point bending tests with higher-strength DP800 were conducted by Heinzelmann et al., yielding comparable results [8, 9].

However, the aforementioned studies also showed that introducing embossing as part of the development process of sheet metal components poses a challenge. This is mainly due to the fact that there is still no standardized workflow for developing embossed sheet metal components. When designing embossed components, additional embossing simulations must be conducted before the conventional metal forming simulation. To perform such simulations, solid elements are recommended to accurately represent the influences of embossing on the sheet metal structure and properties. Thus, use of solid element formulation allows the complex stress states introduced by near-surface embossing to be accurately modelled. This is of crucial importance as these stresses play a major role in the development process and associated simulations. [8]. To use the results of the embossing simulations for subsequent forming simulations, these also have to be performed with solid elements. This leads to two disadvantages: Firstly, the calculation of solid element simulations requires more computational time than equivalent shell element simulations. Secondly, solvers for the simulation of sheet metal components are optimised to calculate with shell elements only [10].

One approach to overcome this problem is to transfer the results of the simulation with solid element mesh to a shell element mesh [11]. When mapping from a solid mesh to a shell mesh, the stresses and strains of the solid elements are transcribed to the integration points of the corresponding shell elements. This allows the results of embossing simulations to be transferred to the shell elements of a forming simulation. However, this approach still involves a solid element simulation for the embossing. To further reduce the calculation time, the embossing simulation therefore needs to be further optimised. The maximum reduction would occur if only the forming simulation is necessary without a prior embossing simulation. Therefore, this paper presents a new approach that avoids the need for prior embossing simulations. This approach simulates the embossed areas with different material properties than the rest of the component. This means that the embossing simulation can be completely avoided and the computation time for the development of embossed components is reduced. This approach assumes that

* the introduced changes in mechanical properties are homogeneously distributed over the embossing area and that the material properties depend on the embossing pattern,

* the geometry of the final component is not affected by the embossing, and

* the transition between embossed and unembossed component areas can be discretised, i.e. neighbouring elements can have different material properties.

## Method

In this paper, three simulation workflows are compared regarding their computational time for forming process with embossed sheet metal materials. The first approach represents the current state-of-the-art for designing near-surface embossed parts. Here, a forming simulation is preceded by an embossing simulation which adds the results (stresses, strains, and sheet-metal thickness) of the embossing to the initial sheet material of the forming simulation before it starts. This requires both the embossing simulation and the process simulation to be performed with solid elements. The approach will therefore be referred to as the “solid approach”. The second approach also includes a preceded embossing simulation, but its results are mapped to a shell mesh. This mapped mesh is used for the forming simulation without further modification. Therefore, the forming simulation does not have to be carried out with solid elements but can retain shell element formulation. Further on, this approach will be mentioned as “mapping approach”. To further improve computational efficiency, the paper describes a final approach that eliminates the need
Additionally, three-point bending tests with embossed DP500 of Briesenick et al. showed an increase in the maximum bending force [7]. Complementary to the work of Walzer and Briesenick, tensile tests and three-point bending tests with higher-strength DP800 were conducted by Heinzelmann et al., yielding comparable results [8, 9].

However, the aforementioned studies also showed that introducing embossing as part of the development process of sheet metal components poses a challenge. This is mainly due to the fact that there is still no standardized workflow for developing embossed sheet metal components. When designing embossed components, additional embossing simulations must be conducted before the conventional metal forming simulation. To perform such simulations, solid elements are recommended to accurately represent the influences of embossing on the sheet metal structure and properties. Thus, use of solid element formulation allows the complex stress states introduced by near-surface embossing to be accurately modelled. This is of crucial importance as these stresses play a major role in the development process and associated simulations. [8]. To use the results of the embossing simulations for subsequent forming simulations, these also have to be performed with solid elements. This leads to two disadvantages: Firstly, the calculation of solid element simulations requires more computational time than equivalent shell element simulations. Secondly, solvers for the simulation of sheet metal components are optimised to calculate with shell elements only [10].

One approach to overcome this problem is to transfer the results of the simulation with solid element mesh to a shell element mesh [11]. When mapping from a solid mesh to a shell mesh, the stresses and strains of the solid elements are transcribed to the integration points of the corresponding shell elements. This allows the results of embossing simulations to be transferred to the shell elements of a forming simulation. However, this approach still involves a solid element simulation for the embossing. To further reduce the calculation time, the embossing simulation therefore needs to be further optimised. The maximum reduction would occur if only the forming simulation is necessary without a prior embossing simulation. Therefore, this paper presents a new approach that avoids the need for prior embossing simulations. This approach simulates the embossed areas with different material properties than the rest of the component. This means that the embossing simulation can be completely avoided and the computation time for the development of embossed components is reduced. This approach assumes that

* the introduced changes in mechanical properties are homogeneously distributed over the embossing area and that the material properties depend on the embossing pattern,

* the geometry of the final component is not affected by the embossing, and

* the transition between embossed and unembossed component areas can be discretised, i.e. neighbouring elements can have different material properties.

## Method

In this paper, three simulation workflows are compared regarding their computational time for forming process with embossed sheet metal materials. The first approach represents the current state-of-the-art for designing near-surface embossed parts. Here, a forming simulation is preceded by an embossing simulation which adds the results (stresses, strains, and sheet-metal thickness) of the embossing to the initial sheet material of the forming simulation before it starts. This requires both the embossing simulation and the process simulation to be performed with solid elements. The approach will therefore be referred to as the “solid approach”. The second approach also includes a preceded embossing simulation, but its results are mapped to a shell mesh. This mapped mesh is used for the forming simulation without further modification. Therefore, the forming simulation does not have to be carried out with solid elements but can retain shell element formulation. Further on, this approach will be mentioned as “mapping approach”. To further improve computational efficiency, the paper describes a final approach that eliminates the need
for prior embossing simulation. Here, the embossed areas are characterised solely by specific material properties in the modelling of the sheet metal material for the forming simulation. The material properties of the embossed areas differ from the rest of the part. This approach will be referred to as the “material approach”.

The comparison of the aforementioned approaches was carried out using two common practical examples, the tensile test and the three-point bending test. The tensile test thereby represented the “ideal” uniaxial stress state. The three-point bending test was used to investigate a more complex stress state. Both tests were performed simultaneously and experimentally. To quantify the influence of the embossing, the tests were performed with specimens having different embossing patterns. All embossings were applied in a regular hexagonal pattern, which is defined by the embossing distance *a* (distance from one embossing to the respective neighbouring embossings) and the depth of the embossings. For the study carried out, distances *a* of 1.0 mm, 1.5 mm and 2.0 mm were investigated. Each pattern was applied to the sheet material with two embossing depths *t*<sub>emb</sub> (depth of penetration of the embossing punch measured from the surface to be embossed) of 100 µm and 200 µm (see Figure 1). This resulted in six different embossing patterns. For the tests, the dual-

<figure><img src="image_2139.png" /><figcaption>Figure 1: General overview of the embossed specimens with representation of the dimensions and the variables.</figcaption></figure>
for prior embossing simulation. Here, the embossed areas are characterised solely by specific material properties in the modelling of the sheet metal material for the forming simulation. The material properties of the embossed areas differ from the rest of the part. This approach will be referred to as the “material approach”.

The comparison of the aforementioned approaches was carried out using two common practical examples, the tensile test and the three-point bending test. The tensile test thereby represented the “ideal” uniaxial stress state. The three-point bending test was used to investigate a more complex stress state. Both tests were performed simultaneously and experimentally. To quantify the influence of the embossing, the tests were performed with specimens having different embossing patterns. All embossings were applied in a regular hexagonal pattern, which is defined by the embossing distance *a* (distance from one embossing to the respective neighbouring embossings) and the depth of the embossings. For the study carried out, distances *a* of 1.0 mm, 1.5 mm and 2.0 mm were investigated. Each pattern was applied to the sheet material with two embossing depths *t*<sub>emb</sub> (depth of penetration of the embossing punch measured from the surface to be embossed) of 100 µm and 200 µm (see Figure 1). This resulted in six different embossing patterns. For the tests, the dual-

<figure><img src="image_2130.png" /><figcaption>Figure 1: General overview of the embossed specimens with representation of the dimensions and the variables.</figcaption></figure>
phase (DP) steel DP800 with a sheet thickness of 1.5 mm, which is commonly used in the automotive industry, was considered.

Test specimens for the real experiments were produced with a Trumpf TRUPunch 5000R CNC punching machine. This machine is capable of up to 1.400 strokes per minute. A round punch with a diameter *d* of 4 mm and the front surface milled into a hemisphere was used to apply the embossing patterns. For the tensile specimens, the embossing patterns were first applied onto a sheet metal blank, from which the final geometry was milled afterwards in accordance with DIN 50125-H20×80 [12]. The specimens for the three-point bending tests were first embossed and then punched to 60 x 60 mm. The tensile and three-point bending tests were performed on a Roell + Korthaus RKM100 universal testing machine. The tensile tests were carried out according to DIN EN ISO 6892-1 [13] and the three-point bending tests according to DIN EN ISO 7438 [14]. During all tests, the force and displacement curves were recorded (see Figure 2).

<figure><img src="image_2140.png" /><figcaption>Figure 2: Extrapolated flow curves of the six embossing patterns and the unembossed reference.</figcaption></figure>

For simulating the tests according to the real counterpart, the LS-DYNA solver was used. All simulations were performed using the material model "*MAT_024", as it represents the material sufficiently accurately. The flow curves used to describe the hardening behaviour are shown in Figure 2. The flow curves were included in the simulation as tabular data. The simulations were performed on an Intel® Core™ i7-7700 processor (8 cores @ 3.6 GHz). The simulation models of all approaches consisted of three parts: embossing punch (rigid), die (rigid) and specimen. The rigid parts were modelled using shell elements. As a representation of the punch, a spherical segment and for the die a plane surface was chosen. The edges' nodes of the specimen were fixed translationally and rotationally in all spatial dimensions. The movement of the punch and die was synchronised for embossing in the x and y direction. Further, in the z direction, these parts moved opposingly. Depending on the embossing pattern, the embossing punch penetrated the specimen 100 µm respectively 200 µm. The contact formulation "*CONTACT ONE WAY SURFACE_TO SURFACE FORMING" was used to model the contact between the specimen and the tool surfaces. For the solid approach, the specimen geometry was meshed using solid elements with an element edge length of 0.5 mm. The results of the embossing simulations were written as a "dynain" file and subsequently used for the forming simulations. For the shell approach, the specimen's geometry was meshed with shell elements. The results of the solid approach's embossing simulation were mapped to the shell mesh of the forming simulation using the mapping software Envyo®. For the material approach, an additional embossing simulation was not required. The material model was changed globally for the experiments regarding this paper. However, for each embossing pattern, a new material card had
phase (DP) steel DP800 with a sheet thickness of 1.5 mm, which is commonly used in the automotive industry, was considered.

Test specimens for the real experiments were produced with a Trumpf TRUPunch 5000R CNC punching machine. This machine is capable of up to 1.400 strokes per minute. A round punch with a diameter *d* of 4 mm and the front surface milled into a hemisphere was used to apply the embossing patterns. For the tensile specimens, the embossing patterns were first applied onto a sheet metal blank, from which the final geometry was milled afterwards in accordance with DIN 50125-H20×80 [12]. The specimens for the three-point bending tests were first embossed and then punched to 60 x 60 mm. The tensile and three-point bending tests were performed on a Roell + Korthaus RKM100 universal testing machine. The tensile tests were carried out according to DIN EN ISO 6892-1 [13] and the three-point bending tests according to DIN EN ISO 7438 [14]. During all tests, the force and displacement curves were recorded (see Figure 2).

<figure><img src="image_2131.png" /><figcaption>Figure 2: Extrapolated flow curves of the six embossing patterns and the unembossed reference.</figcaption></figure>

For simulating the tests according to the real counterpart, the LS-DYNA solver was used. All simulations were performed using the material model "*MAT_024", as it represents the material sufficiently accurately. The flow curves used to describe the hardening behaviour are shown in Figure 2. The flow curves were included in the simulation as tabular data. The simulations were performed on an Intel® Core™ i7-7700 processor (8 cores @ 3.6 GHz). The simulation models of all approaches consisted of three parts: embossing punch (rigid), die (rigid) and specimen. The rigid parts were modelled using shell elements. As a representation of the punch, a spherical segment and for the die a plane surface was chosen. The edges' nodes of the specimen were fixed translationally and rotationally in all spatial dimensions. The movement of the punch and die was synchronised for embossing in the x and y direction. Further, in the z direction, these parts moved opposingly. Depending on the embossing pattern, the embossing punch penetrated the specimen 100 µm respectively 200 µm. The contact formulation "*CONTACT ONE WAY SURFACE_TO SURFACE FORMING" was used to model the contact between the specimen and the tool surfaces. For the solid approach, the specimen geometry was meshed using solid elements with an element edge length of 0.5 mm. The results of the embossing simulations were written as a "dynain" file and subsequently used for the forming simulations. For the shell approach, the specimen's geometry was meshed with shell elements. The results of the solid approach's embossing simulation were mapped to the shell mesh of the forming simulation using the mapping software Envyo®. For the material approach, an additional embossing simulation was not required. The material model was changed globally for the experiments regarding this paper. However, for each embossing pattern, a new material card had
to be created for the forming simulation in order to match the specific material properties. Since the material model "*MAT_024" can be described with only a flow curve, the experimental data of the embossed tensile tests were used to generate the custom material cards. The flow curves were created using the stress-strain values of the tensile tests. As these values are only valid up to the yield point, the curves were extrapolated for higher degrees of deformation [15]. The extrapolation was performed according to the Hockett-Sherby method [16]. The flow curves respectively the material cards were generated for the unembossed material and the six embossing patterns (see Figure 2). The simulations were set up according to DIN EN ISO 6892-1 for the tensile tests and DIN EN ISO 7438 for the three-point bending tests. Like the experiments, the force and displacement curves were calculated in the simulation as well. The force in the tensile test was determined using a cross-section plane. It was placed in the centre of the specimen to calculate the force of the middle elements. For the three-point bending test, the force was given by the reaction force between the bending punch and the specimen. The displacements in both tests were determined by the node displacement.

## Results

To compare the approaches, the computation time required for the simulations was used as a benchmark. Figure 3 shows the average calculation time for the three approaches. All tensile test simulations, including the various embossing patterns and the reference test, were utilized. This figure shows the average calculation times of the different embossing patterns.

The mapping approach is 15.8%, and the material approach is 90.5% faster than the solid approach. When comparing the mapping approach with the material approach, the material approach is 87.5% faster than the mapping approach. This clearly shows how large a part of the total calculation time is due to the embossing simulation and how great the advantages of the material approach are compared to the approaches with a prior embossing simulation.

<figure><img src="image_2141.png" /><figcaption>Figure 3: Evaluation of the computation time of the three approaches for the different embossing geometries and the tensile tests.</figcaption></figure>

The assessment of the different approaches was based on force-displacement curves derived from three-point bending tests. Figure 4 illustrates these curves, which include both experimental data and results from the solid, mapping and material approaches for a specific parameter set
to be created for the forming simulation in order to match the specific material properties. Since the material model "*MAT_024" can be described with only a flow curve, the experimental data of the embossed tensile tests were used to generate the custom material cards. The flow curves were created using the stress-strain values of the tensile tests. As these values are only valid up to the yield point, the curves were extrapolated for higher degrees of deformation [15]. The extrapolation was performed according to the Hockett-Sherby method [16]. The flow curves respectively the material cards were generated for the unembossed material and the six embossing patterns (see Figure 2). The simulations were set up according to DIN EN ISO 6892-1 for the tensile tests and DIN EN ISO 7438 for the three-point bending tests. Like the experiments, the force and displacement curves were calculated in the simulation as well. The force in the tensile test was determined using a cross-section plane. It was placed in the centre of the specimen to calculate the force of the middle elements. For the three-point bending test, the force was given by the reaction force between the bending punch and the specimen. The displacements in both tests were determined by the node displacement.

## Results

To compare the approaches, the computation time required for the simulations was used as a benchmark. Figure 3 shows the average calculation time for the three approaches. All tensile test simulations, including the various embossing patterns and the reference test, were utilized. This figure shows the average calculation times of the different embossing patterns.

The mapping approach is 15.8 %, and the material approach is 90.5 % faster than the solid approach. When comparing the mapping approach with the material approach, the material approach is 87.5 % faster than the mapping approach. This clearly shows how large a part of the total calculation time is due to the embossing simulation and how great the advantages of the material approach are compared to the approaches with a prior embossing simulation.

<figure><img src="image_2132.png" /><figcaption>Figure 3: Evaluation of the computation time of the three approaches for the different embossing geometries and the tensile tests.</figcaption></figure>

The assessment of the different approaches was based on force-displacement curves derived from three-point bending tests. Figure 4 illustrates these curves, which include both experimental data and results from the solid, mapping and material approaches for a specific parameter set
(h = 100 µm, a = 2 µm). The qualitative evaluation involved comparing how well each approach agreed with the observed experimental results. Initially, simulations of the three-point bending test were set up for the unembossed case using both solid and shell elements. Apart from the use of different elements, the simulations did not differ. These simulations were then adapted for the embossed specimens and hasn't been fitted to the reality, so the computational times were comparable. As a result, there are deviations to the experiment. It's worth noting that the embossing process, which takes place before the test, induces deformation in the specimen. Consequently, the deformed specimens do not rest fully on the supports. The sub-figures in Figure 4 show the differences in the specimens, with the solid approach (upper left) taking this deformation into account, whereas the mapping and material approaches (upper right) do not. As a result, the curves of the experiment and the solid approach show good similarities in the early stages, as both the mapping and material approaches were simulated with a flat, undeformed shell mesh, while the solid approach used a deformed solid mesh. In contrast, the two shell approaches (mapping and material approach) start with full contact, skipping the initial deformation required to achieve full contact with the supports and bending punch. Additionally, challenges associated with the embossing process were identified but not considered in the simulation, including slight variations in embossing depth, changes in geometry after embossing, and interactions of the embossings. Reproducing these proved challenging and time-consuming, particularly when using a simple material model. Furthermore, the bending punch and the support is considered rigid, and the stiffness of the universal testing machine and tools used were not considered, which may not accurately represent real-world conditions, as a universal testing machine is not perfectly rigid.

<figure><img src="image_2142.png" /><figcaption>Figure 4: Example (h = 100 µm, a = 2 mm) of the force displacement curves of the experiment, solid approach, mapping approach, and material approach.</figcaption></figure>
(h = 100 µm, a = 2 µm). The qualitative evaluation involved comparing how well each approach agreed with the observed experimental results. Initially, simulations of the three-point bending test were set up for the unembossed case using both solid and shell elements. Apart from the use of different elements, the simulations did not differ. These simulations were then adapted for the embossed specimens and hasn't been fitted to the reality, so the computational times were comparable. As a result, there are deviations to the experiment. It's worth noting that the embossing process, which takes place before the test, induces deformation in the specimen. Consequently, the deformed specimens do not rest fully on the supports. The sub-figures in Figure 4 show the differences in the specimens, with the solid approach (upper left) taking this deformation into account, whereas the mapping and material approaches (upper right) do not. As a result, the curves of the experiment and the solid approach show good similarities in the early stages, as both the mapping and material approaches were simulated with a flat, undeformed shell mesh, while the solid approach used a deformed solid mesh. In contrast, the two shell approaches (mapping and material approach) start with full contact, skipping the initial deformation required to achieve full contact with the supports and bending punch. Additionally, challenges associated with the embossing process were identified but not considered in the simulation, including slight variations in embossing depth, changes in geometry after embossing, and interactions of the embossings. Reproducing these proved challenging and time-consuming, particularly when using a simple material model. Furthermore, the bending punch and the support is considered rigid, and the stiffness of the universal testing machine and tools used were not considered, which may not accurately represent real-world conditions, as a universal testing machine is not perfectly rigid.

<figure><img src="image_2133.png" /><figcaption>Figure 4: Example (h = 100 µm, a = 2 mm) of the force displacement curves of the experiment, solid approach, mapping approach, and material approach.</figcaption></figure>
## Conclusion and Outlook

In this paper, two new approaches for the forming simulation of embossed components were considered for and compared with the conventional simulation method with regard to the required computing time. For the conventional solid approach, both the embossing simulation and the process simulations were considered. For the mapping approach, the results of an embossing simulation with solid elements (solid approach) were mapped to the integration points of the forming simulations' shell mesh. Since the additional computation time of the process simulation is only due to the preceding embossing simulation, the last material approach attempted to dispense the embossing simulation by approximating embossed areas by selectively adjusting the material properties in the forming simulation. This required different material cards for the embossed and unembossed areas.

In order to compare the approaches with each other, two different material tests were carried out both simulatively and experimentally with selectively embossed specimens. Here, tensile tests according to DIN EN ISO 6892-1 and three-point bending tests according to DIN EN ISO 7438 were performed. A total of six embossing patterns were investigated, varying in embossing distance (1 mm, 1.5 mm, and 2 mm) and depth (100 µm and 200 µm). As no additional embossing simulation was required, the calculation time for the material approach revealed as the shortest. The other two approaches each include a preceding solid embossing simulation, which makes them significantly slower. As the solid approach also involves adapting the forming simulation to a solid representation of the specimen and an additional embossing simulation, this approach is the slowest. To minimise this aspect, the mapping approach was proposed, where the results of the solid simulation were mapped to the integration points of the shell elements of the forming simulation. This led to a reduction of the computation time, but was far from achieving the short calculation time of the material approach. To compare the approaches, only the computation time was considered in this paper. The preparation and adaptation of the simulations, the mapping from solid to shell elements, and the creation of the material cards were not considered in terms of time in any way. Difficulties in quantifying the individual preparation steps are the reason for this. Qualitatively, all approaches were able to reproduce the experiments well. At the start of the force-displacement curve, the solid approach was able to reproduce the initial setting of the specimen better than the other two approaches due to the additional change in geometry introduced by the embossing. The course of the experiment is also slightly better represented by the solid approach.

The approaches presented here to model embossed areas offers enormous potential for speeding up the design and development of embossed sheet metal components. The fundamental advantages have been described and quantified for the first time in this paper. Simple tests with simple stress states have been considered. Whether this approach is also suitable for more complex design methods used in industrial practice remains to be verified. The simulations considered here were also carried out with a simple material model. For more complex examples, it will also be necessary to adapt this model. This will also result in more complex material characterization tests to be carried out on embossed and unembossed materials. If for more complex cases it is found that the approximation with the material approach is not sufficiently fast or accurate, it would also be conceivable to adapt the mapping approach. For example, this approach could be used to simulate a single embossing with solid elements and the results could be copied and mapped into the corresponding embossing pattern.

## Acknowledgements

The research project "Herstellung crashoptimierter Blechbauteile mittels Prägen" of the Europäische Forschungsvereinigung EFB e. V. is carried out in the framework of the industrial collective research programme (IGF no. 21940N). It is supported by the Federal Ministry for Economic Affairs and Climate Action (BMWK) through the AiF (German Federation of Industrial Research Associations e. V) based on a decision taken by the German Bundestag. Additionally,
## Conclusion and Outlook

In this paper, two new approaches for the forming simulation of embossed components were considered for and compared with the conventional simulation method with regard to the required computing time. For the conventional solid approach, both the embossing simulation and the process simulations were considered. For the mapping approach, the results of an embossing simulation with solid elements (solid approach) were mapped to the integration points of the forming simulations' shell mesh. Since the additional computation time of the process simulation is only due to the preceding embossing simulation, the last material approach attempted to dispense the embossing simulation by approximating embossed areas by selectively adjusting the material properties in the forming simulation. This required different material cards for the embossed and unembossed areas.

In order to compare the approaches with each other, two different material tests were carried out both simulatively and experimentally with selectively embossed specimens. Here, tensile tests according to DIN EN ISO 6892-1 and three-point bending tests according to DIN EN ISO 7438 were performed. A total of six embossing patterns were investigated, varying in embossing distance (1 mm, 1.5 mm, and 2 mm) and depth (100 µm and 200 µm). As no additional embossing simulation was required, the calculation time for the material approach revealed as the shortest. The other two approaches each include a preceding solid embossing simulation, which makes them significantly slower. As the solid approach also involves adapting the forming simulation to a solid representation of the specimen and an additional embossing simulation, this approach is the slowest. To minimise this aspect, the mapping approach was proposed, where the results of the solid simulation were mapped to the integration points of the shell elements of the forming simulation. This led to a reduction of the computation time, but was far from achieving the short calculation time of the material approach. To compare the approaches, only the computation time was considered in this paper. The preparation and adaptation of the simulations, the mapping from solid to shell elements, and the creation of the material cards were not considered in terms of time in any way. Difficulties in quantifying the individual preparation steps are the reason for this. Qualitatively, all approaches were able to reproduce the experiments well. At the start of the force-displacement curve, the solid approach was able to reproduce the initial setting of the specimen better than the other two approaches due to the additional change in geometry introduced by the embossing. The course of the experiment is also slightly better represented by the solid approach.

The approaches presented here to model embossed areas offers enormous potential for speeding up the design and development of embossed sheet metal components. The fundamental advantages have been described and quantified for the first time in this paper. Simple tests with simple stress states have been considered. Whether this approach is also suitable for more complex design methods used in industrial practice remains to be verified. The simulations considered here were also carried out with a simple material model. For more complex examples, it will also be necessary to adapt this model. This will also result in more complex material characterization tests to be carried out on embossed and unembossed materials. If for more complex cases it is found that the approximation with the material approach is not sufficiently fast or accurate, it would also be conceivable to adapt the mapping approach. For example, this approach could be used to simulate a single embossing with solid elements and the results could be copied and mapped into the corresponding embossing pattern.

## Acknowledgements

The research project "Herstellung crashoptimierter Blechbauteile mittels Prägen" of the Europäische Forschungsvereinigung EFB e. V. is carried out in the framework of the industrial collective research programme (IGF no. 21940N). It is supported by the Federal Ministry for Economic Affairs and Climate Action (BMWK) through the AiF (German Federation of Industrial Research Associations e. V) based on a decision taken by the German Bundestag. Additionally,
we would like to thank our project partners DYNAmore GmbH, Eckold GmbH & Co. KG, synchropress GmbH, Engineering System International GmbH, Karosseriewerke Dresden GmbH, Zeller + Gmelin GmbH & Co. KG, Novelis Switzerland SA, Tata Steel Europe, STRACK NORMA GmbH & Co. KG, F+K Werkstoffprüfung und Labor GmbH, and KIRCHHOFF Automotive GmbH.

## References

[1] K. Shanmugam, V. Gadhamshetty, P. Yadav, D. Athanassiadis, M. Tysklind, and V. K. Upadhyayula, "Advanced High-Strength Steel and Carbon Fiber Reinforced Polymer Composite Body in White for Passenger Cars: Environmental Performance and Sustainable Return on Investment under Different Propulsion Modes," *ACS Sustainable Chem. Eng.*, vol. 7, no. 5, pp. 4951-4963, 2019. https://doi.org/10.1021/acssuschemeng.8b05588

[2] H. E. Friedrich, *Leichtbau in der Fahrzeugtechnik*. Wiesbaden: Springer Fachmedien Wiesbaden, 2017.

[3] C. S. Namoco, T. Iizuka, R. C. Sagrado, N. Takakura, and K. Yamaguchi, "Experimental and numerical investigation of restoration behavior of sheet metals subjected to bulging deformation," *Journal of Materials Processing Technology*, vol. 177, 1-3, pp. 368-372, 2006. https://doi.org/10.1016/j.jmatprotec.2006.03.208

[4] C. S. Namoco, T. Iizuka, K. Narita, N. Takakura, and K. Yamaguchi, "Effects of embossing and restoration process on the deep drawability of aluminum alloy sheets," *Journal of Materials Processing Technology*, vol. 187-188, pp. 202-206, 2007. https://doi.org/10.1016/j.jmatprotec.2006.11.182

[5] Y. Abe, K. Mori, T. Maeno, S. Ishihara, and Y. Kato, "Improvement of sheet metal formability by local work-hardening with punch indentation," *Prod. Eng. Res. Devel.*, vol. 13, no. 5, pp. 589-597, 2019. https://doi.org/10.1007/s11740-019-00910-6

[6] S. Walzer, M. Liewald, N. Simon, J. Gibmeier, H. Erdle, and T. Böhlke, "Improvement of Sheet Metal Properties by Inducing Residual Stresses into Sheet Metal Components by Embossing and Reforming," *Applied Science and Engineering Progress*, vol. 15, no. 1, 2021. https://doi.org/10.14416/j.asep.2021.09.006

[7] D. Briesenick, S. Walzer, and M. Liewald, "Study on the Effect of Embossing on the Bending Properties of High-Strength Sheet Metals," in *Forming the Future: Proceedings of the 13th International Conference on the Technology of Plasticity* (The Minerals, Metals & Materials Series), G. Daehn, J. Cao, B. Kinsey, E. Tekkaya, A. Vivek, and Y. Yoshida, Eds., 1st ed. Cham: Springer International Publishing; Imprint Springer, 2021, pp. 2585-2595

[8] P. Heinzelmann, D. Briesenick, and M. Liewald, "Characterization of the Mechanical Properties of Selectively Embossed Sheet Metal Materials Under Multi-Axial Loads," *IOP Conf. Ser.: Mater. Sci. Eng.*, vol. 1284, no. 1, p. 12043, 2023. https://doi.org/10.1088/1757-899X/1284/1/012043

[9] P. Heinzelmann, D. Briesenick, and M. Liewald, "Modifying mechanical properties of sheet metal materials by work hardening mechanisms induced by selective embossing," in *Material Forming*, Lukasz Madej, Mateusz Sitko, Konrad Perzynski, Ed., 2023, pp. 951-958. https://doi.org/10.21741/9781644902479-104

[10] E. Madenci, *The Finite Element Method and Applications in Engineering Using ANSYS®*, 2nd ed. (SpringerLink Bücher). Boston, MA: Springer, 2015.
we would like to thank our project partners DYNAmore GmbH, Eckold GmbH & Co. KG, synchropress GmbH, Engineering System International GmbH, Karosseriewerke Dresden GmbH, Zeller + Gmelin GmbH & Co. KG, Novelis Switzerland SA, Tata Steel Europe, STRACK NORMA GmbH & Co. KG, F+K Werkstoffprüfung und Labor GmbH, and KIRCHHOFF Automotive GmbH.

## References

[1] K. Shanmugam, V. Gadhamshetty, P. Yadav, D. Athanassiadis, M. Tysklind, and V. K. Upadhyayula, “Advanced High-Strength Steel and Carbon Fiber Reinforced Polymer Composite Body in White for Passenger Cars: Environmental Performance and Sustainable Return on Investment under Different Propulsion Modes,” *ACS Sustainable Chem. Eng.*, vol. 7, no. 5, pp. 4951–4963, 2019. https://doi.org/10.1021/acssuschemeng.8b05588

[2] H. E. Friedrich, *Leichtbau in der Fahrzeugtechnik*. Wiesbaden: Springer Fachmedien Wiesbaden, 2017.

[3] C. S. Namoco, T. Iizuka, R. C. Sagrado, N. Takakura, and K. Yamaguchi, “Experimental and numerical investigation of restoration behavior of sheet metals subjected to bulging deformation,” *Journal of Materials Processing Technology*, vol. 177, 1-3, pp. 368–372, 2006. https://doi.org/10.1016/j.jmatprotec.2006.03.208

[4] C. S. Namoco, T. Iizuka, K. Narita, N. Takakura, and K. Yamaguchi, “Effects of embossing and restoration process on the deep drawability of aluminum alloy sheets,” *Journal of Materials Processing Technology*, vol. 187-188, pp. 202–206, 2007. https://doi.org/10.1016/j.jmatprotec.2006.11.182

[5] Y. Abe, K. Mori, T. Maeno, S. Ishihara, and Y. Kato, “Improvement of sheet metal formability by local work-hardening with punch indentation,” *Prod. Eng. Res. Devel.*, vol. 13, no. 5, pp. 589–597, 2019. https://doi.org/10.1007/s11740-019-00910-6

[6] S. Walzer, M. Liewald, N. Simon, J. Gibmeier, H. Erdle, and T. Böhlke, “Improvement of Sheet Metal Properties by Inducing Residual Stresses into Sheet Metal Components by Embossing and Reforming,” *Applied Science and Engineering Progress*, vol. 15, no. 1, 2021. https://doi.org/10.14416/j.asep.2021.09.006

[7] D. Briesenick, S. Walzer, and M. Liewald, “Study on the Effect of Embossing on the Bending Properties of High-Strength Sheet Metals,” in *Forming the Future: Proceedings of the 13th International Conference on the Technology of Plasticity* (The Minerals, Metals & Materials Series), G. Daehn, J. Cao, B. Kinsey, E. Tekkaya, A. Vivek, and Y. Yoshida, Eds., 1st ed. Cham: Springer International Publishing; Imprint Springer, 2021, pp. 2585–2595

[8] P. Heinzelmann, D. Briesenick, and M. Liewald, “Characterization of the Mechanical Properties of Selectively Embossed Sheet Metal Materials Under Multi-Axial Loads,” *IOP Conf. Ser.: Mater. Sci. Eng.*, vol. 1284, no. 1, p. 12043, 2023. https://doi.org/10.1088/1757-899X/1284/1/012043

[9] P. Heinzelmann, D. Briesenick, and M. Liewald, “Modifying mechanical properties of sheet metal materials by work hardening mechanisms induced by selective embossing,” in *Material Forming*, Lukasz Madej, Mateusz Sitko, Konrad Perzynski, Ed., 2023, pp. 951–958. https://doi.org/10.21741/9781644902479-104

[10] E. Madenci, *The Finite Element Method and Applications in Engineering Using ANSYS®*, 2nd ed. (SpringerLink Bücher). Boston, MA: Springer, 2015.