and influence on the cutting process", Applied Thermal Engineering, vol. 50, no. 1, pp. 1311-1325, 2013. https://doi.org/10.1016/j.applthermaleng.2012.06.047

[7] S. P. F. Jaspers and J. Dautzenberg, "Material behaviour in conditions similar to metal cutting: flow stress in the primary shear zone," *J. Mater. Process. Technol.*, vol. 122, no. 2–3, pp. 322–330, Mar. 2002. https://doi.org/10.1016/S0924-0136(01)01228-6

[8] V. M. Luchesi and R. T. Coelho, "Experimental investigations of heat transfer coefficients of cutting fluids in metal cutting processes: Analysis of workpiece phenomena in a given case study," *Proc. Inst. Mech. Eng. Part B J. Eng. Manuf.*, vol. 226, no. 7, pp. 1174–1184, 2012. https://doi.org/10.1177/0954405412442459

[9] S. Kurgin, J. M. Dasch, D. L. Simon, G. C. Barber, and Q. Zou, "Evaluation of the convective heat transfer coefficient for minimum quantity lubrication (MQL)," *Ind. Lubr. Tribol.*, vol. 64, no. 6, pp. 376–386, 2012. https://doi.org/10.1108/00368791211262516

[10] J. Rech, P. J. Arrazola, C. Claudin, C. Courbon, F. Pusavec, and J. Kopac, "Characterisation of friction and heat partition coefficients at the tool-work material interface in cutting," *CIRP Ann. - Manuf. Technol.*, vol. 62, no. 1, pp. 79–82, 2013. https://doi.org/10.1016/j.cirp.2013.03.099

[11] E. Usui, T. Shirakashi, and T. Kitagawa, "Analytical Prediction of Three Dimensional Cutting Process—Part 3: Cutting Temperature and Crater Wear of Carbide Tool," *J. Eng. Ind.*, vol. 100, no. 2, p. 236, 1978. https://doi.org/10.1115/1.3439415

[12] E. Usui, T. Shirakashi, and T. Kitagawa, "Analytical prediction of cutting tool wear," *Wear*, vol. 100, no. 1–3, pp. 129–151, 1984. https://doi.org/10.1016/0043-1648(84)90010-3

[13] K. Maekawa, T. Kitagawa, T. Shirakashi, and E. Usui, "Analytical Prediction of Flank Wear of Carbide Tools in Turning Plain Carbon Steels (Part 1)- Characteristic Equation of Flank Wear," *Bull. Japan Soc. Precis. Eng.*, vol. 23, no. 2, pp. 126–133, 1988.

[14] K. Maekawa, T. Kitagawa, T. Shirakashi, and E. Usui, "Analytical prediction of flank wear of carbide tools in turning plain carbon steel. Part 2. Prediction of flank wear," *Bull. Japan Soc. Precis. Eng.*, vol. 23, no. 2, pp. 126–133, 1989.

[15] F. Klocke and P. Frank, "Simulation of tool Wear in hard turning," *Proceeding 9th CIRP Int. Work. Model. Mach. Oper.*, 2006.

[16] Y. C. Yen, J. Söhner, B. Lilly, and T. Altan, "Estimation of tool wear in orthogonal cutting using the finite element analysis," *J. Mater. Process. Technol.*, vol. 146, no. 1, pp. 82–91, 2004. https://doi.org/10.1016/S0924-0136(03)00847-1

[17] N. Senthilkumar and T. Tamizharasan, "Effect of Tool Geometry in Turning AISI 1045 Steel: Experimental Investigation and FEM Analysis," *Arab. J. Sci. Eng.*, vol. 39, no. 6, pp. 4963–4975, 2014. https://doi.org/10.1007/s13369-014-1054-2
and influence on the cutting process", *Applied Thermal Engineering*, vol. 50, no. 1, pp. 1311-1325, 2013. https://doi.org/10.1016/j.applthermaleng.2012.06.047

[7] S. P. F. Jaspers and J. Dautzenberg, "Material behaviour in conditions similar to metal cutting: flow stress in the primary shear zone," *J. Mater. Process. Technol.*, vol. 122, no. 2-3, pp. 322-330, Mar. 2002. https://doi.org/10.1016/S0924-0136(01)01228-6

[8] V. M. Luchesi and R. T. Coelho, "Experimental investigations of heat transfer coefficients of cutting fluids in metal cutting processes: Analysis of workpiece phenomena in a given case study," *Proc. Inst. Mech. Eng. Part B J. Eng. Manuf.*, vol. 226, no. 7, pp. 1174-1184, 2012. https://doi.org/10.1177/0954405412442459

[9] S. Kurgin, J. M. Dasch, D. L. Simon, G. C. Barber, and Q. Zou, "Evaluation of the convective heat transfer coefficient for minimum quantity lubrication (MQL)," *Ind. Lubr. Tribol.*, vol. 64, no. 6, pp. 376-386, 2012. https://doi.org/10.1108/00368791211262516

[10] J. Rech, P. J. Arrazola, C. Claudin, C. Courbon, F. Pusavec, and J. Kopac, "Characterisation of friction and heat partition coefficients at the tool-work material interface in cutting," *CIRP Ann. - Manuf. Technol.*, vol. 62, no. 1, pp. 79-82, 2013. https://doi.org/10.1016/j.cirp.2013.03.099

[11] E. Usui, T. Shirakashi, and T. Kitagawa, "Analytical Prediction of Three Dimensional Cutting Process—Part 3: Cutting Temperature and Crater Wear of Carbide Tool," *J. Eng. Ind.*, vol. 100, no. 2, p. 236, 1978. https://doi.org/10.1115/1.3439415

[12] E. Usui, T. Shirakashi, and T. Kitagawa, "Analytical prediction of cutting tool wear," *Wear*, vol. 100, no. 1-3, pp. 129-151, 1984. https://doi.org/10.1016/0043-1648(84)90010-3

[13] K. Maekawa, T. Kitagawa, T. Shirakashi, and E. Usui, "Analytical Prediction of Flank Wear of Carbide Tools in Turning Plain Carbon Steels (Part 1)- Characteristic Equation of Flank Wear," *Bull. Japan Soc. Precis. Eng.*, vol. 23, no. 2, pp. 126-133, 1988.

[14] K. Maekawa, T. Kitagawa, T. Shirakashi, and E. Usui, "Analytical prediction of flank wear of carbide tools in turning plain carbon steel. Part 2. Prediction of flank wear," *Bull. Japan Soc. Precis. Eng.*, vol. 23, no. 2, pp. 126-133, 1989.

[15] F. Klocke and P. Frank, "Simulation of tool Wear in hard turning," *Proceeding 9th CIRP Int. Work. Model. Mach. Oper.*, 2006.

[16] Y. C. Yen, J. Söhner, B. Lilly, and T. Altan, "Estimation of tool wear in orthogonal cutting using the finite element analysis," *J. Mater. Process. Technol.*, vol. 146, no. 1, pp. 82-91, 2004. https://doi.org/10.1016/S0924-0136(03)00847-1

[17] N. Senthilkumar and T. Tamizharasan, "Effect of Tool Geometry in Turning AISI 1045 Steel: Experimental Investigation and FEM Analysis," *Arab. J. Sci. Eng.*, vol. 39, no. 6, pp. 4963-4975, 2014. https://doi.org/10.1007/s13369-014-1054-2
Material behaviour modelling

<figure><img src="image_2001.png" /></figure>
Material behaviour modelling

<figure><img src="image_1996.png" /></figure>
# Investigations of strain rate sensitivity under different stress triaxialities for DC04

MAIER Lorenz¹,ᵃ*, MARKER Edgar¹,ᵇ, SCHULZ Fabian¹,ᶜ and VOLK Wolfram¹,ᵈ

¹Walther-Meissner-Street 4, 85748 Garching near Munich, Germany

ᵃlorenz.maier@utg.de, ᵇedgar.marker@tum.de, ᶜfabi.schulz@tum.de, ᵈwolfram.volk@utg.de

**Keywords:** Strain Rate Sensitivity, Sheet Metal Forming, Deep Drawing, Stress Triaxiality

**Abstract.** Understanding the strain rate sensitivity of materials is essential for predicting their behavior in sheet metal forming. While uniaxial tension tests are state of the art in characterizing this sensitivity, the deformation response of materials under different loading conditions can significantly deviate from uniaxial behavior. This paper presents a comprehensive study of the strain rate sensitivity of DC04 through a series of experimental investigations with different strain rates. In addition to uniaxial tension tests, the study investigates the strain rate sensitivity under shear and plane strain tests, providing a comprehensive analysis of strain rate sensitivity across different loading scenarios. The investigation aims to understand how the material responds to varying deformation rates, focusing on characterizing their deformation behavior under various loading conditions. The authors collected experimental data from the material with a DIC system. They analyzed it to derive material-specific parameters that describe their strain rate-dependent responses depending on the stress state. To explain this, the authors calibrated three models: Johnson Cook, Cowper Symonds, and Huh Kang.

**Introduction**

Efficiency and precision are crucial in modern manufacturing, especially in producing sheet metal parts and deep drawing processes [1]. Variability in batch and process conditions often results in substandard components that cause material loss and waste of resources. At a time when sustainability and resource conservation play a significant role, it is essential to find ways to reduce this waste in production. This requires a solid ability to predict the properties of the resulting parts amid such material and environmental variations. [1–3]

One promising approach lies in applying advanced constitutive models for sheet metal forming simulations. Their calibration enables accurate prediction of the process [4]. One field of those advanced constitutive models deals with the strain rate sensitivity of sheet metal materials. The calibration of these constitutive models usually requires experiments with different strain rates [5]. However, the reliability of these models and therefore simulations depends on the accuracy of the defined material parameters themselves and the comprehensive changeability of the model parameters. [1]

**State of the art**

Constitutive models are used to explain the behavior of materials during sheet metal forming. These models vary in complexity and calibration effort. In these models, the sensitivity to strain rate e.g. can also be represented, which is often achieved by calibration through uniaxial tensile experiment [3]. With this experiment Shi and Meuleman [6] also found out, that with increasing plastic deformation the strain rate sensitivity decreases. This was also observed by Jocham et al. [7]. They also found that the strain rate sensitivity is dependent on the present loading condition.

The strain rate sensitivity models, primarily tailored to emulate crash behavior, are often adjusted to exceptionally high strain rates [5, 8] and high temperature effects [9]. The focus is

<figure><img src="image_2002.png" /></figure>
# Investigations of strain rate sensitivity under different stress triaxialities for DC04

MAIER Lorenz¹,ᵃ*, MARKER Edgar¹,ᵇ, SCHULZ Fabian¹,ᶜ and VOLK Wolfram¹,ᵈ

¹Walther-Meissner-Street 4, 85748 Garching near Munich, Germany

ᵃlorenz.maier@utg.de, ᵇedgar.marker@tum.de, ᶜfabi.schulz@tum.de, ᵈwolfram.volk@utg.de

**Keywords:** Strain Rate Sensitivity, Sheet Metal Forming, Deep Drawing, Stress Triaxiality

**Abstract.** Understanding the strain rate sensitivity of materials is essential for predicting their behavior in sheet metal forming. While uniaxial tension tests are state of the art in characterizing this sensitivity, the deformation response of materials under different loading conditions can significantly deviate from uniaxial behavior. This paper presents a comprehensive study of the strain rate sensitivity of DC04 through a series of experimental investigations with different strain rates. In addition to uniaxial tension tests, the study investigates the strain rate sensitivity under shear and plane strain tests, providing a comprehensive analysis of strain rate sensitivity across different loading scenarios. The investigation aims to understand how the material responds to varying deformation rates, focusing on characterizing their deformation behavior under various loading conditions. The authors collected experimental data from the material with a DIC system. They analyzed it to derive material-specific parameters that describe their strain rate-dependent responses depending on the stress state. To explain this, the authors calibrated three models: Johnson Cook, Cowper Symonds, and Huh Kang.

**Introduction**

Efficiency and precision are crucial in modern manufacturing, especially in producing sheet metal parts and deep drawing processes [1]. Variability in batch and process conditions often results in substandard components that cause material loss and waste of resources. At a time when sustainability and resource conservation play a significant role, it is essential to find ways to reduce this waste in production. This requires a solid ability to predict the properties of the resulting parts amid such material and environmental variations. [1–3]

One promising approach lies in applying advanced constitutive models for sheet metal forming simulations. Their calibration enables accurate prediction of the process [4]. One field of those advanced constitutive models deals with the strain rate sensitivity of sheet metal materials. The calibration of these constitutive models usually requires experiments with different strain rates [5]. However, the reliability of these models and therefore simulations depends on the accuracy of the defined material parameters themselves and the comprehensive changeability of the model parameters. [1]

**State of the art**

Constitutive models are used to explain the behavior of materials during sheet metal forming. These models vary in complexity and calibration effort. In these models, the sensitivity to strain rate e.g. can also be represented, which is often achieved by calibration through uniaxial tensile experiment [3]. With this experiment Shi and Meuleman [6] also found out, that with increasing plastic deformation the strain rate sensitivity decreases. This was also observed by Jocham et al. [7]. They also found that the strain rate sensitivity is dependent on the present loading condition.

The strain rate sensitivity models, primarily tailored to emulate crash behavior, are often adjusted to exceptionally high strain rates [5, 8] and high temperature effects [9]. The focus is

<figure><img src="image_1997.png" /></figure>
usually not on the forming behavior in the manufacturing process [2]. Instead, the failure at different strain rates is being considered. [10, 11] Nevertheless, strain rate sensitivity also has an impact during the manufacturing process [1].

To investigate the influence of different stress triaxialities, Roth and Mohr [5] conducted experiments with varying loading velocities on specimens of high-strength steel sheets to examine the impact of loading velocity on the onset of fracture. Different stress triaxialities were achieved by using varying specimen geometries. The results show that the deformability increases by up to 60 % when the crosshead velocity is increased from 7.6 × 10⁻⁶ m/s to 4.2 m/s. A plastic model considering different loading rates and materials is proposed to provide an accurate description of the deformation and onset of fracture. Roth and Mohr used a split Hopkinson bar to achieve high forming speeds. [5, 8]

Vallaster et al. [1] calibrated a constitutive model for an HC340LA with uniaxial tensile tests and three strain rates from 0.0004 s⁻¹ to 0.04 s⁻¹. One aspect of their work was the validation of the simulations with strain rate sensitivity by comparing simulated strain fields with experimentally achieved results. Local changes in the strain rates were achieved by varying the blank holder force and the associated change in flow behavior. During the validation process, it was found that considering the strain rate sensitivity in the simulation could improve the representation of the force behavior. However, the simulation time was extended. [1]

One possibility to model the strain rate sensitivity is the Johnson Cook equation [9]. Here, strain rate sensitivity performs a central role in capturing the response of a material to deformation at different rates. The function's first term describes the plastic deformation with the material parameters A, B and n. This equation introduces the second term that incorporates strain rate (ε̇) in relation to the lowest strain rate (ε̇₀) to represent the change in yield stress of the material at different strain rates:

$$ \sigma(\varepsilon, \dot{\varepsilon}, T) = (A + B \cdot \varepsilon^n) \cdot \left(1 + C \ln \frac{\dot{\varepsilon}}{\varepsilon_0}\right) \cdot \left(1 - \left(\frac{T-T_0}{T_m-T_0}\right)^m\right) \quad (1) $$

This component accounts for the sensitivity of the material to the rate at which it is deformed and reflects real-world behavior where materials react differently when subjected to fast or slow deformation processes. The function's third term describes the material's temperature behavior. [9]

In addition to this model, the Cowper-Symonds model is often used, particularly for very high strain rates, with the strain rate normalized to the lowest experimental value [12]. The model proposed by Cowper and Symonds defines as follows:

$$ \sigma(\varepsilon, \dot{\varepsilon}) = k(\varepsilon_0 + \varepsilon)^n \cdot \left(1 + \left(\frac{\varepsilon^*}{C}\right)^{\frac{1}{P}}\right) \quad (2) $$

The Cowper Symonds model offers several advantages in predicting material behavior at high deformation rates. It is a useful choice in material modeling due to its simplicity and ability to describe deformation and strength behavior at high speeds. [1, 12]

Compared to the Cowper Symonds model, the term of the Huh Kang [13] model provides precise description of the material behavior at high deformation rates and shows an improved fit to experimental data [13]. The Johnson Cook model is modified by this term proposed by Huh and Kang as follows:

$$ \sigma(\varepsilon, \dot{\varepsilon}, T) = (A + B \cdot \varepsilon^n) \cdot \left(1 + C_1 \ln \frac{\dot{\varepsilon}}{\varepsilon_0} + C_2 (\ln \frac{\dot{\varepsilon}}{\varepsilon_0})^2\right) \cdot \left(1 - \left(\frac{T-T_0}{T_m-T_0}\right)^m\right) \quad (3) $$
usually not on the forming behavior in the manufacturing process [2]. Instead, the failure at different strain rates is being considered. [10, 11] Nevertheless, strain rate sensitivity also has an impact during the manufacturing process [1].

To investigate the influence of different stress triaxialities, Roth and Mohr [5] conducted experiments with varying loading velocities on specimens of high-strength steel sheets to examine the impact of loading velocity on the onset of fracture. Different stress triaxialities were achieved by using varying specimen geometries. The results show that the deformability increases by up to 60 % when the crosshead velocity is increased from 7.6 × 10⁻⁶ m/s to 4.2 m/s. A plastic model considering different loading rates and materials is proposed to provide an accurate description of the deformation and onset of fracture. Roth and Mohr used a split Hopkinson bar to achieve high forming speeds. [5, 8]

Vallaster et al. [1] calibrated a constitutive model for an HC340LA with uniaxial tensile tests and three strain rates from 0.0004 s⁻¹ to 0.04 s⁻¹. One aspect of their work was the validation of the simulations with strain rate sensitivity by comparing simulated strain fields with experimentally achieved results. Local changes in the strain rates were achieved by varying the blank holder force and the associated change in flow behavior. During the validation process, it was found that considering the strain rate sensitivity in the simulation could improve the representation of the force behavior. However, the simulation time was extended. [1]

One possibility to model the strain rate sensitivity is the Johnson Cook equation [9]. Here, strain rate sensitivity performs a central role in capturing the response of a material to deformation at different rates. The function's first term describes the plastic deformation with the material parameters A, B and n. This equation introduces the second term that incorporates strain rate (ε̇) in relation to the lowest strain rate (ε̇₀) to represent the change in yield stress of the material at different strain rates:

$$ \sigma(\varepsilon, \dot{\varepsilon}, T) = (A + B \cdot \varepsilon^n) \cdot \left(1 + C \ln \frac{\dot{\varepsilon}}{\varepsilon_0}\right) \cdot \left(1 - \left(\frac{T-T_0}{T_m-T_0}\right)^m\right) \quad (1) $$

This component accounts for the sensitivity of the material to the rate at which it is deformed and reflects real-world behavior where materials react differently when subjected to fast or slow deformation processes. The function's third term describes the material's temperature behavior. [9]

In addition to this model, the Cowper-Symonds model is often used, particularly for very high strain rates, with the strain rate normalized to the lowest experimental value [12]. The model proposed by Cowper and Symonds defines as follows:

$$ \sigma(\varepsilon, \dot{\varepsilon}) = k(\varepsilon_0 + \varepsilon)^n \cdot \left(1 + \left(\frac{\varepsilon^*}{C}\right)^{\frac{1}{P}}\right) \quad (2) $$

The Cowper Symonds model offers several advantages in predicting material behavior at high deformation rates. It is a useful choice in material modeling due to its simplicity and ability to describe deformation and strength behavior at high speeds. [1, 12]

Compared to the Cowper Symonds model, the term of the Huh Kang [13] model provides precise description of the material behavior at high deformation rates and shows an improved fit to experimental data [13]. The Johnson Cook model is modified by this term proposed by Huh and Kang as follows:

$$ \sigma(\varepsilon, \dot{\varepsilon}, T) = (A + B \cdot \varepsilon^n) \cdot \left(1 + C_1 \ln \frac{\dot{\varepsilon}}{\varepsilon_0} + C_2 (\ln \frac{\dot{\varepsilon}}{\varepsilon_0})^2\right) \cdot \left(1 - \left(\frac{T-T_0}{T_m-T_0}\right)^m\right) \quad (3) $$
## Experimental setup

The experiments were performed using an AllRound Line Z 150 universal testing machine from ZwickRoell GmbH & Co. KG, which offers a wide range of testing options and customizations. This versatile machine precisely controls various test parameters such as load speed, tensile and compressive forces, and records the elongation and load during the tests. The universal testing machine allows testing speeds from 0.0001 to 900 mm/min. The elongation is detected by the non-contact strain measurement instrument LaserXtens from ZwickRoell GmbH & Co. KG and is transmitted to the control system of the testing machine. This setup enables engineering strain rate controlled operation of the machine. The flexibility of the universal testing machine allows the testing of different materials, specimens and a range of test conditions to be simulated to gain a comprehensive understanding of material behavior under various loading conditions. The authors used different crosshead velocities, as is state of the art, to achieve different strain rates. The velocities are shown in Table 1. This range of speeds allowed the recording of force-displacement curves at different strain rates. The digital image correlation (DIC) system ARAMIS SRX from Carl Zeiss GOM Metrology GmbH was used to check the in-situ existing triaxiality of the sample geometry. The authors also used this system to calculate the effective von Mises strain rate of the sample. An iterative procedure adjusted the crosshead speed until the desired strain rate was achieved. This was necessary because the area of interest in the shear specimen is very small compared to the other specimens, and the LaserXtens offers insufficient resolution. The DIC system was used here to adjust the strain rate.

The authors conducted experiments with three different specimen geometries. This allowed the investigation of the strain rate sensitivity under variable triaxiality states. Table 2 shows the sample geometries. The authors carried out three valid repetitions for each test performed. This allowed any outliers in the force or displacement range to be identified and removed from the test series. The authors then averaged the three valid tests. The following figures show the averaged curves for each strain rate. The experimental data was processed by calculating the mean values of the curves and dividing the stresses and forces by the corresponding values of the quasi-static curve (0.00007 s⁻¹).

<table><caption>Table 1 Investigated strain rates for the three different geometries and the corresponding crosshead velocities.</caption><thead><tr><td>Strain rate</td><td>Tensile test</td><td>Shear test</td><td>Plane strain test</td></tr></thead><tbody><tr><td>0.00007 s<sup>-1</sup></td><td>Strain rate controlled</td><td>0.0002 mm/s</td><td>Strain rate controlled</td></tr><tr><td>0.002 s<sup>-1</sup></td><td>Strain rate controlled</td><td>0.0032 mm/s</td><td>Strain rate controlled</td></tr><tr><td>0.020 s<sup>-1</sup></td><td>Strain rate controlled</td><td>0.032 mm/s</td><td>Strain rate controlled</td></tr><tr><td>0.125 s<sup>-1</sup></td><td>15.0 mm/s</td><td>0.27 mm/s</td><td>15.0 mm/s</td></tr><tr><td>1.748 s<sup>-1</sup></td><td>-</td><td>15.0 mm/s</td><td>-</td></tr></tbody></table>
## Experimental setup

The experiments were performed using an AllRound Line Z 150 universal testing machine from ZwickRoell GmbH & Co. KG, which offers a wide range of testing options and customizations. This versatile machine precisely controls various test parameters such as load speed, tensile and compressive forces, and records the elongation and load during the tests. The universal testing machine allows testing speeds from 0.0001 to 900 mm/min. The elongation is detected by the non-contact strain measurement instrument LaserXtens from ZwickRoell GmbH & Co. KG and is transmitted to the control system of the testing machine. This setup enables engineering strain rate controlled operation of the machine. The flexibility of the universal testing machine allows the testing of different materials, specimens and a range of test conditions to be simulated to gain a comprehensive understanding of material behavior under various loading conditions. The authors used different crosshead velocities, as is state of the art, to achieve different strain rates. The velocities are shown in Table 1. This range of speeds allowed the recording of force-displacement curves at different strain rates. The digital image correlation (DIC) system ARAMIS SRX from Carl Zeiss GOM Metrology GmbH was used to check the in-situ existing triaxiality of the sample geometry. The authors also used this system to calculate the effective von Mises strain rate of the sample. An iterative procedure adjusted the crosshead speed until the desired strain rate was achieved. This was necessary because the area of interest in the shear specimen is very small compared to the other specimens, and the LaserXtens offers insufficient resolution. The DIC system was used here to adjust the strain rate.

The authors conducted experiments with three different specimen geometries. This allowed the investigation of the strain rate sensitivity under variable triaxiality states. Table 2 shows the sample geometries. The authors carried out three valid repetitions for each test performed. This allowed any outliers in the force or displacement range to be identified and removed from the test series. The authors then averaged the three valid tests. The following figures show the averaged curves for each strain rate. The experimental data was processed by calculating the mean values of the curves and dividing the stresses and forces by the corresponding values of the quasi-static curve (0.00007 s⁻¹).

<table><caption>Table 1 Investigated strain rates for the three different geometries and the corresponding crosshead velocities.</caption><thead><tr><td>Strain rate</td><td>Tensile test</td><td>Shear test</td><td>Plane strain test</td></tr></thead><tbody><tr><td>0.00007 s<sup>-1</sup></td><td>Strain rate controlled</td><td>0.0002 mm/s</td><td>Strain rate controlled</td></tr><tr><td>0.002 s<sup>-1</sup></td><td>Strain rate controlled</td><td>0.0032 mm/s</td><td>Strain rate controlled</td></tr><tr><td>0.020 s<sup>-1</sup></td><td>Strain rate controlled</td><td>0.032 mm/s</td><td>Strain rate controlled</td></tr><tr><td>0.125 s<sup>-1</sup></td><td>15.0 mm/s</td><td>0.27 mm/s</td><td>15.0 mm/s</td></tr><tr><td>1.748 s<sup>-1</sup></td><td>-</td><td>15.0 mm/s</td><td>-</td></tr></tbody></table>
<table><caption>Table 2 Experiment geometries conducted with the corresponding triaxiality. The strain range is marked in the respective figure of the specimen geometry.</caption><thead><tr><th></th><th>Tensile test</th><th>Shear test</th><th>Plane strain test</th></tr></thead><tbody><tr><th>Triaxiality</th><td>1/3</td><td>0.00</td><td>1/√3</td></tr><tr><th>Reference</th><td>DIN 50125 [14]</td><td>Roth et al. [8]</td><td>Khameneh et al. [2]</td></tr><tr><th>Rolling direction</th><td>↑</td><td>↗</td><td>↑</td></tr></tbody></table>

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

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

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

## Results

First, the authors carried out the uniaxial tensile test. They examined four different strain rates for this test. The resulting stress-strain curves can be seen in Fig. 1 a). A striking feature of these tests is that the LaserXtens does not allow continuous control with increasing strain rates. This occurs for the first time at a strain rate of 0.02 s⁻¹. At a strain rate of 0.125 s⁻¹, a constant crosshead speed was used, as can be seen in Table 1. At this speed, the inertia and acceleration of the LaserXtens influences the recording of the strain at very low deformation. This must be considered when further evaluating the curves.

For this DC04 it can be determined in Fig. 1 b) that the influence of strain rate sensitivity decreases with increasing plastic deformation. This is also what Jocham et al. [7] found out for their DC06. To calibrate the strain rate sensitivity models, the authors of this work used the mean value of the respective stresses of the investigated strain rates at a strain of 8-12 %. On the one hand, this can reduce the effect of the inertia of the LaserXtens. On the other hand, the noise is filtered out of the curve at a strain rate of 0.02 s⁻¹.
<table><caption>Table 2 Experiment geometries conducted with the corresponding triaxiality. The strain range is marked in the respective figure of the specimen geometry.</caption><thead><tr><th></th><th>Tensile test</th><th>Shear test</th><th>Plane strain test</th></tr></thead><tbody><tr><th>Triaxiality</th><td>1/3</td><td>0.00</td><td>1/√3</td></tr><tr><th>Reference</th><td>DIN 50125 [14]</td><td>Roth et al. [8]</td><td>Khameneh et al. [2]</td></tr><tr><th>Rolling direction</th><td>↑</td><td>↑</td><td>↑</td></tr></tbody></table>

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

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

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

## Results

First, the authors carried out the uniaxial tensile test. They examined four different strain rates for this test. The resulting stress-strain curves can be seen in Fig. 1 a). A striking feature of these tests is that the LaserXtens does not allow continuous control with increasing strain rates. This occurs for the first time at a strain rate of 0.02 s⁻¹. At a strain rate of 0.125 s⁻¹, a constant crosshead speed was used, as can be seen in Table 1. At this speed, the inertia and acceleration of the LaserXtens influences the recording of the strain at very low deformation. This must be considered when further evaluating the curves.

For this DC04 it can be determined in Fig. 1 b) that the influence of strain rate sensitivity decreases with increasing plastic deformation. This is also what Jocham et al. [7] found out for their DC06. To calibrate the strain rate sensitivity models, the authors of this work used the mean value of the respective stresses of the investigated strain rates at a strain of 8-12 %. On the one hand, this can reduce the effect of the inertia of the LaserXtens. On the other hand, the noise is filtered out of the curve at a strain rate of 0.02 s⁻¹.
<figure><img src="image_2006.png" /><figcaption>Figure 1 a) stress-strain curves of the uniaxial tensile test at four different strain rates in s⁻¹. The strain rate control was realized using the extension measurement system or a constant crosshead velocity. b) Normalized curves to the lowest strain rate.</figcaption></figure>

In addition to the uniaxial stress state, other stress states and their influence on the strain rate sensitivity were investigated. In the next step, the authors investigate the specimen geometry for the plane strain state from Table 2. The authors also conducted these investigations on the universal testing machine. The resulting force-displacement curves are shown in Fig. 2. The strain reference length for the LaserXtens centered on the specimen was 10 mm. The equivalent strain rate was verified by the DIC system. The inertia effect of the strain measurement system can again be observed in the curves. To filter out this effect, the force ratios were averaged for a displacement of eight to ten millimeters. The normalized force-displacement curves stabilize with these changes in length, as can be seen in Fig. 2 a). The curves are presented until the first curve of the experiments has reached the force maximum.

<figure><img src="image_2007.png" /><figcaption>Figure 2 a) force-displacement curves of the plane strain geometry at four different strain rates in s⁻¹. The strain rate was achieved using the extension measurement system or a constant crosshead velocity. b) Normalized force-displacement curves to the lowest strain rate.</figcaption></figure>

In addition to the comparative strain rates, the authors also utilized the DIC system to verify the desired strain state. The result of this investigation is shown in Fig. 3. By doing so, the authors
<figure><img src="image_2001.png" /><figcaption>Figure 1 a) stress-strain curves of the uniaxial tensile test at four different strain rates in s⁻¹. The strain rate control was realized using the extension measurement system or a constant crosshead velocity. b) Normalized curves to the lowest strain rate.</figcaption></figure>

In addition to the uniaxial stress state, other stress states and their influence on the strain rate sensitivity were investigated. In the next step, the authors investigate the specimen geometry for the plane strain state from Table 2. The authors also conducted these investigations on the universal testing machine. The resulting force-displacement curves are shown in Fig. 2. The strain reference length for the LaserXtens centered on the specimen was 10 mm. The equivalent strain rate was verified by the DIC system. The inertia effect of the strain measurement system can again be observed in the curves. To filter out this effect, the force ratios were averaged for a displacement of eight to ten millimeters. The normalized force-displacement curves stabilize with these changes in length, as can be seen in Fig. 2 a). The curves are presented until the first curve of the experiments has reached the force maximum.

<figure><img src="image_2002.png" /><figcaption>Figure 2 a) force-displacement curves of the plane strain geometry at four different strain rates in s⁻¹. The strain rate was achieved using the extension measurement system or a constant crosshead velocity. b) Normalized force-displacement curves to the lowest strain rate.</figcaption></figure>

In addition to the comparative strain rates, the authors also utilized the DIC system to verify the desired strain state. The result of this investigation is shown in Fig. 3. By doing so, the authors
proved that with the sample geometry used, the minor strain in the center of the sample is almost zero. In the specimen center, a plane strain condition was achieved.

<figure><img src="image_2008.png" /><figcaption>Figure 3 Evaluation of the major and minor strain in the evaluation area by the DIC system for the chosen plane strain geometry as well as the corresponding von Mises strain. The false color image shows the equivalent plastic strain according to von Mises.</figcaption></figure>

In addition to the uniaxial and the plane strain state, the authors also carried out shear tests with various strain rates. The corresponding strain rate was verified via the von Mises strain rate using the DIC system and adjusted iteratively. This allowed the authors to perform tests in the same strain rate range of the geometries already investigated. Fig. 4 a) shows the averaged force-displacement curves of the experiments conducted. It can be seen here that an oscillating behavior occurs in the force signal during the quasi-static test series. This may be due to the geometric design of the sample. Due to the two shear zones, deformations may not occur homogeneously and simultaneously in both zones. This is no longer noticeable at higher strain rates. By normalizing, as shown in Fig. 4 b), this behavior is transferred to the other experiments. To filter this behavior, the authors used the mean values of the force ratios in the displacement range from 0.5 mm to 1.0 mm between the upper and lower sections of the sample.

<figure><img src="image_2009.png" /><figcaption>Figure 4 a) force-displacement curves of the shear test at four different strain rates in s⁻¹. The strain rate was analyzed using the DIC system. b) Normalized curves to the lowest strain rate.</figcaption></figure>

The authors used the DIC system to monitor the strain distribution and to verify the desired triaxial state. Fig. 5 shows the von Mises equivalent strain in the false color diagram. The investigated area in the shear zone is also marked. Based on the major and minor strain distributions detected in this way, the authors were able to determine that the desired triaxial state dominates.
proved that with the sample geometry used, the minor strain in the center of the sample is almost zero. In the specimen center, a plane strain condition was achieved.

<figure><img src="image_2003.png" /><figcaption>Figure 3 Evaluation of the major and minor strain in the evaluation area by the DIC system for the chosen plane strain geometry as well as the corresponding von Mises strain. The false color image shows the equivalent plastic strain according to von Mises.</figcaption></figure>

In addition to the uniaxial and the plane strain state, the authors also carried out shear tests with various strain rates. The corresponding strain rate was verified via the von Mises strain rate using the DIC system and adjusted iteratively. This allowed the authors to perform tests in the same strain rate range of the geometries already investigated. Fig. 4 a) shows the averaged force-displacement curves of the experiments conducted. It can be seen here that an oscillating behavior occurs in the force signal during the quasi-static test series. This may be due to the geometric design of the sample. Due to the two shear zones, deformations may not occur homogeneously and simultaneously in both zones. This is no longer noticeable at higher strain rates. By normalizing, as shown in Fig. 4 b), this behavior is transferred to the other experiments. To filter this behavior, the authors used the mean values of the force ratios in the displacement range from 0.5 mm to 1.0 mm between the upper and lower sections of the sample.

<figure><img src="image_2004.png" /><figcaption>Figure 4 a) force-displacement curves of the shear test at four different strain rates in s⁻¹. The strain rate was analyzed using the DIC system. b) Normalized curves to the lowest strain rate.</figcaption></figure>

The authors used the DIC system to monitor the strain distribution and to verify the desired triaxial state. Fig. 5 shows the von Mises equivalent strain in the false color diagram. The investigated area in the shear zone is also marked. Based on the major and minor strain distributions detected in this way, the authors were able to determine that the desired triaxial state dominates.
<figure><img src="image_2010.png" /><figcaption>Figure 5 Evaluation of the major and minor strain in the evaluation area by the DIC system for the chosen shear test geometry as well as the corresponding von Mises strain. The false color image shows the equivalent plastic strain according to von Mises.</figcaption></figure>

Due to the narrow shear zone in this sample geometry, the authors were able to achieve a significantly higher strain rate at maximum crosshead velocity. Here, a strain rate of 1.748 s⁻¹ was achieved at 15 mm/s crosshead velocity. This strain rate could be determined by the DIC system and allowed a further strain rate to be investigated, which is almost 14 times higher compared to the 0.125 s⁻¹. This support point will serve as a validation point in the following.

The authors have tested three triaxiality states using this approach. For each sample geometry, they conducted three valid tests at four different strain rates. Using this database, the authors subsequently calibrated the three strain rate sensitivity models Johnson Cook (Eq. 1), Cowper Symonds (Eq. 2) and Huh Kang (Eq. 3). In Fig. 6, the models are calibrated to the uniaxial tensile test.

<figure><img src="image_2011.png" /><figcaption>Figure 6 Experimental data points from the uniaxial tensile tests with the three strain rate sensitivity models calibrated to them and the validation point from the shear test.</figcaption></figure>

In this calibration method, Huh Kang's model achieves the lowest sum of squares error (SSE) with 2.23 × 10⁻¹¹. The prediction of this model shows a particularly precise point of agreement during validation. Johnson Cook underestimates the stress ratio at the validation point. The SSE of this model is 2.75 × 10⁻⁴. The Cowper Symonds model overestimates the validation point. It has an SSE of 3.65 × 10⁻⁴.
<figure><img src="image_2005.png" /><figcaption>Figure 5 Evaluation of the major and minor strain in the evaluation area by the DIC system for the chosen shear test geometry as well as the corresponding von Mises strain. The false color image shows the equivalent plastic strain according to von Mises.</figcaption></figure>

Due to the narrow shear zone in this sample geometry, the authors were able to achieve a significantly higher strain rate at maximum crosshead velocity. Here, a strain rate of 1.748 s⁻¹ was achieved at 15 mm/s crosshead velocity. This strain rate could be determined by the DIC system and allowed a further strain rate to be investigated, which is almost 14 times higher compared to the 0.125 s⁻¹. This support point will serve as a validation point in the following.

The authors have tested three triaxiality states using this approach. For each sample geometry, they conducted three valid tests at four different strain rates. Using this database, the authors subsequently calibrated the three strain rate sensitivity models Johnson Cook (Eq. 1), Cowper Symonds (Eq. 2) and Huh Kang (Eq. 3). In Fig. 6, the models are calibrated to the uniaxial tensile test.

<figure><img src="image_2006.png" /><figcaption>Figure 6 Experimental data points from the uniaxial tensile tests with the three strain rate sensitivity models calibrated to them and the validation point from the shear test.</figcaption></figure>

In this calibration method, Huh Kang's model achieves the lowest sum of squares error (SSE) with 2.23 × 10⁻¹¹. The prediction of this model shows a particularly precise point of agreement during validation. Johnson Cook underestimates the stress ratio at the validation point. The SSE of this model is 2.75 × 10⁻⁴. The Cowper Symonds model overestimates the validation point. It has an SSE of 3.65 × 10⁻⁴.
Analogous to the procedure for the uniaxial tensile test, the authors used the force ratios at different strain rates in the plane strain test in the following step to calibrate the models. Fig. 7 shows the calibrated models and the force ratios.

<figure><img src="image_2012.png" /><figcaption>Figure 7 Experimental data points from the plane strain tests with the three strain rate sensitivity models calibrated to them and the validation point from the shear test.</figcaption></figure>

The implementation of the experiment encountered unforeseen complications. In this series of tests, the force ratios do not consistently remain higher, they exhibit a greater degree of variability. Despite these circumstances, however, it was possible to obtain analyzable results. Although the SSE of Huh Kang is significantly higher with $1.66 \times 10^{-4}$, it is still the lowest compared to Johnson Cook ($3.18 \times 10^{-4}$) and Cowper Symonds ($6.23 \times 10^{-4}$). Once more, Huh Kang's model exhibits the closest proximity to the validation point.

In the final step, the authors calibrated the three strain rate sensitivity models using the three strain rates of the shear test. The result of this calibration is shown graphically in Fig. 8. In the shear tests, the force ratios, similar to the plane strain experiments, deviate further from a linear arrangement in logarithmic space compared to the tensile experiments, as can be seen in Fig. 8. Consequently, the SSE of the models are also higher than those of the tensile test. Once again, the SSE of the Huh Kang model is the lowest with $2.02 \times 10^{-4}$ compared to Johnson Cook ($4.14 \times 10^{-4}$) and Cowper Symonds ($8.83 \times 10^{-4}$).

As in previous model calibrations, the Huh Kang model reaches the validation point most precisely. In this case, even the very best. It remains unclear whether the reason for this is that the validation point, similar to the calibration points, was determined by means of a shear test. However, doubts arise when considering the variance of the calibration points compared to the fit of the model.
Analogous to the procedure for the uniaxial tensile test, the authors used the force ratios at different strain rates in the plane strain test in the following step to calibrate the models. Fig. 7 shows the calibrated models and the force ratios.

<figure><img src="image_2007.png" /><figcaption>Figure 7 Experimental data points from the plane strain tests with the three strain rate sensitivity models calibrated to them and the validation point from the shear test.</figcaption></figure>

The implementation of the experiment encountered unforeseen complications. In this series of tests, the force ratios do not consistently remain higher, they exhibit a greater degree of variability. Despite these circumstances, however, it was possible to obtain analyzable results. Although the SSE of Huh Kang is significantly higher with $1.66 \times 10^{-4}$, it is still the lowest compared to Johnson Cook ($3.18 \times 10^{-4}$) and Cowper Symonds ($6.23 \times 10^{-4}$). Once more, Huh Kang's model exhibits the closest proximity to the validation point.

In the final step, the authors calibrated the three strain rate sensitivity models using the three strain rates of the shear test. The result of this calibration is shown graphically in Fig. 8. In the shear tests, the force ratios, similar to the plane strain experiments, deviate further from a linear arrangement in logarithmic space compared to the tensile experiments, as can be seen in Fig. 8. Consequently, the SSE of the models are also higher than those of the tensile test. Once again, the SSE of the Huh Kang model is the lowest with $2.02 \times 10^{-4}$ compared to Johnson Cook ($4.14 \times 10^{-4}$) and Cowper Symonds ($8.83 \times 10^{-4}$).

As in previous model calibrations, the Huh Kang model reaches the validation point most precisely. In this case, even the very best. It remains unclear whether the reason for this is that the validation point, similar to the calibration points, was determined by means of a shear test. However, doubts arise when considering the variance of the calibration points compared to the fit of the model.
<figure><img src="image_2013.png" /><figcaption>Figure 8 Experimental data points from the shear test with the three strain rate sensitivity models calibrated to them and the validation point also from the shear test.</figcaption></figure>

The model parameters from the respective strain rate sensitivity models are shown in Table 3 depending on the calibration experiment.

<table><caption>Table 3 Parameters determined from the Eq. 1, Eq. 2, and Eq. 3 for the three investigated strain rate sensitivity models of the respective experiments.</caption><thead><tr><th rowspan="2"></th><th colspan="1">Johnson Cook</th><th colspan="2">Cowper Symonds</th><th colspan="2">Huh Kang</th></tr><tr><th>C</th><th>C</th><th>P</th><th>C<sub>1</sub></th><th>C<sub>2</sub></th></tr></thead><tbody><tr><th>Tensile test</th><td>0.01743</td><td>91.06</td><td>3.407</td><td>0.00968</td><td>0.001207</td></tr><tr><th>Plane strain test</th><td>0.02032</td><td>80.16</td><td>3.599</td><td>0.01456</td><td>0.000897</td></tr><tr><th>Shear test</th><td>0.01891</td><td>99.12</td><td>3.572</td><td>0.01209</td><td>0.001061</td></tr></tbody></table>

## Conclusion

The authors carried out tests on a universal testing machine with a DC04. Three different strain conditions were achieved using different geometries: uniaxial strain, plane strain and pure shear. A DIC system was utilized to investigate the strain state. These three strain states were performed at four different strain rates between 0.00007 s⁻¹ and 0.125 s⁻¹ with three repetitions each. In addition, the shear test was carried out at a strain rate of 1.748 s⁻¹. This point served as a validation point for the calibrated strain rate sensitivity models.

As already stated by Shi and Meuleman as well as Jocham et al., a decrease in strain rate sensitivity with increasing plastic deformation was observed for the uniaxial tensile tests. The authors could not observe the plasticity dependence of strain rate sensitivity for other triaxiality states (plane strain and pure shear). This was not to be expected, as a significantly higher equivalent plastic strain was achieved in shear experiment compared to the uniaxial tensile test.

During model calibration, it was found that the Huh Kang model proves to be suitable for modeling strain rate sensitivity in the relevant strain rate range of production. The authors demonstrate this by comparison with the validation point from the shear test. The validation point was best matched here for all three strain states investigated.

In contrast to the uniaxial tensile test, no decrease in strain rate sensitivity was observed with higher plastification for other triaxialities. It therefore remains of interest to examine whether a dependence on the triaxiality can be observed at high strain rates. In this respect, the investigation did not reveal any clear findings at the strain rates considered. A model that can depict this would be particularly helpful for crash simulations.
<figure><img src="image_2008.png" /><figcaption>Figure 8 Experimental data points from the shear test with the three strain rate sensitivity models calibrated to them and the validation point also from the shear test.</figcaption></figure>

The model parameters from the respective strain rate sensitivity models are shown in Table 3 depending on the calibration experiment.

<table><caption>Table 3 Parameters determined from the Eq. 1, Eq. 2, and Eq. 3 for the three investigated strain rate sensitivity models of the respective experiments.</caption><thead><tr><th rowspan="2"></th><th colspan="1">Johnson Cook</th><th colspan="2">Cowper Symonds</th><th colspan="2">Huh Kang</th></tr><tr><th>C</th><th>C</th><th>P</th><th>C<sub>1</sub></th><th>C<sub>2</sub></th></tr></thead><tbody><tr><th>Tensile test</th><td>0.01743</td><td>91.06</td><td>3.407</td><td>0.00968</td><td>0.001207</td></tr><tr><th>Plane strain test</th><td>0.02032</td><td>80.16</td><td>3.599</td><td>0.01456</td><td>0.000897</td></tr><tr><th>Shear test</th><td>0.01891</td><td>99.12</td><td>3.572</td><td>0.01209</td><td>0.001061</td></tr></tbody></table>

## Conclusion

The authors carried out tests on a universal testing machine with a DC04. Three different strain conditions were achieved using different geometries: uniaxial strain, plane strain and pure shear. A DIC system was utilized to investigate the strain state. These three strain states were performed at four different strain rates between 0.00007 s⁻¹ and 0.125 s⁻¹ with three repetitions each. In addition, the shear test was carried out at a strain rate of 1.748 s⁻¹. This point served as a validation point for the calibrated strain rate sensitivity models.

As already stated by Shi and Meuleman as well as Jocham et al., a decrease in strain rate sensitivity with increasing plastic deformation was observed for the uniaxial tensile tests. The authors could not observe the plasticity dependence of strain rate sensitivity for other triaxiality states (plane strain and pure shear). This was not to be expected, as a significantly higher equivalent plastic strain was achieved in shear experiment compared to the uniaxial tensile test.

During model calibration, it was found that the Huh Kang model proves to be suitable for modeling strain rate sensitivity in the relevant strain rate range of production. The authors demonstrate this by comparison with the validation point from the shear test. The validation point was best matched here for all three strain states investigated.

In contrast to the uniaxial tensile test, no decrease in strain rate sensitivity was observed with higher plastification for other triaxialities. It therefore remains of interest to examine whether a dependence on the triaxiality can be observed at high strain rates. In this respect, the investigation did not reveal any clear findings at the strain rates considered. A model that can depict this would be particularly helpful for crash simulations.
## Acknowledgments

The authors would like to thank the German Research Foundation (Deutsche Forschungs- gemeinschaft, DFG) for the financial support under the grant number 455388465.

## References

[1] E. Vallaster, S. Wiesenmayer, M. Merklein, "Effect of a strain rate dependent material modeling of a steel on the prediction accuracy of a numerical deep drawing process," *Prod. Eng. Res. Devel.*, 2023. https://doi.org/10.1007/s11740-023-01222-6

[2] F. Khameneh, A. Abedini, C. Butcher, "Experimental and Numerical Fracture Characterization of DP1180 Steel in Combined Simple Shear and Uniaxial Tension," *Metals*, 7 (2023) 1305. https://doi.org/10.3390/met13071305

[3] D. Anderson, S. Winkler, A. Bardelcik, M. J. Worswick, "Influence of stress triaxiality and strain rate on the failure behavior of a dual-phase DP780 steel," *Materials & Design*, 60 (2014) 198-207. https://doi.org/10.1016/j.matdes.2014.03.073

[4] L. Maier, M. Eder, R. Norz, W. Volk, "Innovative experimental setup for the investigation of material models with regard to strain hardening behavior," *IOP Conf. Ser.: Mater. Sci. Eng.*, 1284, (2023) 12071. https://doi.org/10.1088/1757-899X/1284/1/012071

[5] C. C. Roth, D. Mohr, "Effect of strain rate on ductile fracture initiation in advanced high strength steel sheets: Experiments and modeling," *International Journal of Plasticity*, 56 (2014) 19-44. https://doi.org/10.1016/j.ijplas.2014.01.003

[6] M. F. Shi, D. J. Meuleman, "On certain aspects of strain rate sensitivity of sheet metals," *JMEP*, 4 (1995) 321-333. https://doi.org/10.1007/BF02649070

[7] D. Jocham, R. Norz, W. Volk, "Strain rate sensitivity of DC06 for high strains under biaxial stress in hydraulic bulge test and under uniaxial stress in tensile test," *Int J Mater Form*, 10 (2017) 453-461. https://doi.org/10.1007/s12289-016-1293-8

[8] C. C. Roth, D. Mohr, "Ductile fracture experiments with locally proportional loading histories," *International Journal of Plasticity*, 79 (2016) 328-354. https://doi.org/10.1016/j.ijplas.2015.08.004

[9] Y. Wang, X. Zeng, H. Chen, X. Yang, F. Wang, L. Zeng, "Modified Johnson-Cook constitutive model of metallic materials under a wide range of temperatures and strain rates," *Results in Physics*, 27 (2021) 104498. https://doi.org/10.1016/j.rinp.2021.104498

[10] P. Verleysen, J. Peirs, J. van Slycken, K. Faes, L. Duchene, "Effect of strain rate on the forming behaviour of sheet metals," *Journal of Materials Processing Technology*, 211 (2011) 1457-1464. https://doi.org/10.1016/j.jmatprotec.2011.03.018

[11] Y. Lou, S. Zhang, J. W. Yoon, "Strength modeling of sheet metals from shear to plane strain tension," *International Journal of Plasticity*, 134 (2020) 102813. https://doi.org/10.1016/j.ijplas.2020.102813

[12] A. Škrlec, J. Klemenc, "Estimating the Strain-Rate-Dependent Parameters of the Cowper-Symonds and Johnson-Cook Material Models using Taguchi Arrays," *SV-JME*, 62 (2016) 220-230. https://doi.org/10.5545/sv-jme.2015.3266

[13] H. Huh, W. J. Kang, S. S. Han, "A tension split Hopkinson bar for investigating the dynamic behavior of sheet metals," *Experimental Mechanics*, 42 (2002) 8-17. https://doi.org/10.1007/BF02411046

[14] *Testing of metallic materials: Tensile test pieces*, 50125:2022-08, German institute for standardization.
## Acknowledgments

The authors would like to thank the German Research Foundation (Deutsche Forschungs- gemeinschaft, DFG) for the financial support under the grant number 455388465.

## References

[1] E. Vallaster, S. Wiesenmayer, M. Merklein, "Effect of a strain rate dependent material modeling of a steel on the prediction accuracy of a numerical deep drawing process," *Prod. Eng. Res. Devel.*, 2023. https://doi.org/10.1007/s11740-023-01222-6

[2] F. Khameneh, A. Abedini, C. Butcher, "Experimental and Numerical Fracture Characterization of DP1180 Steel in Combined Simple Shear and Uniaxial Tension," *Metals*, 7 (2023) 1305. https://doi.org/10.3390/met13071305

[3] D. Anderson, S. Winkler, A. Bardelcik, M. J. Worswick, "Influence of stress triaxiality and strain rate on the failure behavior of a dual-phase DP780 steel," *Materials & Design*, 60 (2014) 198-207. https://doi.org/10.1016/j.matdes.2014.03.073

[4] L. Maier, M. Eder, R. Norz, W. Volk, "Innovative experimental setup for the investigation of material models with regard to strain hardening behavior," *IOP Conf. Ser.: Mater. Sci. Eng.*, 1284, (2023) 12071. https://doi.org/10.1088/1757-899X/1284/1/012071

[5] C. C. Roth, D. Mohr, "Effect of strain rate on ductile fracture initiation in advanced high strength steel sheets: Experiments and modeling," *International Journal of Plasticity*, 56 (2014) 19-44. https://doi.org/10.1016/j.ijplas.2014.01.003

[6] M. F. Shi, D. J. Meuleman, "On certain aspects of strain rate sensitivity of sheet metals," *JMEP*, 4 (1995) 321-333. https://doi.org/10.1007/BF02649070

[7] D. Jocham, R. Norz, W. Volk, "Strain rate sensitivity of DC06 for high strains under biaxial stress in hydraulic bulge test and under uniaxial stress in tensile test," *Int J Mater Form*, 10 (2017) 453-461. https://doi.org/10.1007/s12289-016-1293-8

[8] C. C. Roth, D. Mohr, "Ductile fracture experiments with locally proportional loading histories," *International Journal of Plasticity*, 79 (2016) 328-354. https://doi.org/10.1016/j.ijplas.2015.08.004

[9] Y. Wang, X. Zeng, H. Chen, X. Yang, F. Wang, L. Zeng, "Modified Johnson-Cook constitutive model of metallic materials under a wide range of temperatures and strain rates," *Results in Physics*, 27 (2021) 104498. https://doi.org/10.1016/j.rinp.2021.104498

[10] P. Verleysen, J. Peirs, J. van Slycken, K. Faes, L. Duchene, "Effect of strain rate on the forming behaviour of sheet metals," *Journal of Materials Processing Technology*, 211 (2011) 1457-1464. https://doi.org/10.1016/j.jmatprotec.2011.03.018

[11] Y. Lou, S. Zhang, J. W. Yoon, "Strength modeling of sheet metals from shear to plane strain tension," *International Journal of Plasticity*, 134 (2020) 102813. https://doi.org/10.1016/j.ijplas.2020.102813

[12] A. Škrlec, J. Klemenc, "Estimating the Strain-Rate-Dependent Parameters of the Cowper-Symonds and Johnson-Cook Material Models using Taguchi Arrays," *SV-JME*, 62 (2016) 220-230. https://doi.org/10.5545/sv-jme.2015.3266

[13] H. Huh, W. J. Kang, S. S. Han, "A tension split Hopkinson bar for investigating the dynamic behavior of sheet metals," *Experimental Mechanics*, 42 (2002) 8-17. https://doi.org/10.1007/BF02411046

[14] *Testing of metallic materials: Tensile test pieces*, 50125:2022-08, German institute for standardization.
# Fatigue resistance of deep drawn parts: A scale bridging simulative study using representative volume elements and crystal plasticity simulations

FEHLEMANN Niklas<fnref n="1" />,ᵃ*, HENRICH Manuel<fnref n="1" />, MÜLLER Martina<fnref n="2" />, KÖNEMANN Markus<fnref n="1" />, BERGS Thomas<fnref n="2" />,<fnref n="3" />, MÜNSTERMANN Sebastian<fnref n="1" />

**Keywords:** Deep Drawing, Fatigue, Representative Volume Elements, Crystal Plasticity

**Abstract.** The mechanical properties of formed components are determined by the interaction between the microstructure and the load path of the forming process. To investigate and understand these effects, micromechanical simulation concepts can be used, such as statistically Representative Volume Elements (sRVE) coupled with crystal plasticity simulations. This study presents a concept that uses sRVE simulations to quantify the influence of three different deep drawing load paths on the fatigue resistance of DP800 steel. The first step is a scale-bridging simulation approach that employs macroscopic simulations of the deep drawing process to extract the boundary conditions for the sRVE simulations with Damask. Subsequent cyclic loading is then simulated. 50 sRVE are computed for each load path to estimate fatigue resistance based on a Fatigue Indicator Parameter. The results indicate that fatigue resistance increases with increasing deformation-induced strain hardening. Additionally, a positive correlation between the martensitic ligament structures and fatigue resistance was observed.

**Introduction**

To achieve Europe's ambitious climate goals [1], it is increasingly necessary to design and construct structural components more efficiently. This requires continuous improvement of manufacturing processes, such as forming processes, in conjunction with the materials used. Several studies have already demonstrated that a well targeted and damage controlled forming process leads to improved component performance and can ultimately unlock lightweight potential [2,3]. In deep drawing, a critical sheet metal forming process, key parameters for modified process control include variables such as blank holder force, drawing velocity, friction, and the radius of the drawing ring [4,5].

To achieve a precise selection and optimization of materials for different manufacturing processes, it is essential to have a fundamental understanding of materials and their microstructures. Simulation methods are often used to achieve this understanding. In microstructure simulation, numerical crystal plasticity (CP) models are frequently used in conjunction with statistically Representative Volume Elements (sRVE), for example for the modelling of complex materials in [6,7]. Bargmann et al. [8] provide a comprehensive overview of the various methods for generating and applying sRVE. They offer a significant advantage in that they allow for the independent investigation of individual microstructure components, which can be experimentally challenging and often impractical [9]. Pütz et al. [10] quantified the influence of martensitic band structures in DP800 steel on the initiation of martensite fracture

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

<footnote n="1">¹Institute of Metal Forming, RWTH Aachen University, Intzestraße 10, 52072 Aachen, Germany</footnote>

<footnote n="2">²Manufacturing Technology Institute MTI of RWTH Aachen University, Campus-Boulevard 30, 52074 Aachen, Germany</footnote>

<footnote n="3">³Fraunhofer Institute for Production Technology IPT, Steinbachstraße 17, 52074 Aachen, Germany</footnote>

<footnote>ᵃniklas.fehlemann@ibf.rwth-aachen.de</footnote>
# Fatigue resistance of deep drawn parts: A scale bridging simulative study using representative volume elements and crystal plasticity simulations

FEHLEMANN Niklas<fnref n="1" />,ᵃ*, HENRICH Manuel<fnref n="1" />, MÜLLER Martina<fnref n="2" />, KÖNEMANN Markus<fnref n="1" />, BERGS Thomas<fnref n="2" />,<fnref n="3" />, MÜNSTERMANN Sebastian<fnref n="1" />

**Keywords:** Deep Drawing, Fatigue, Representative Volume Elements, Crystal Plasticity

**Abstract.** The mechanical properties of formed components are determined by the interaction between the microstructure and the load path of the forming process. To investigate and understand these effects, micromechanical simulation concepts can be used, such as statistically Representative Volume Elements (sRVE) coupled with crystal plasticity simulations. This study presents a concept that uses sRVE simulations to quantify the influence of three different deep drawing load paths on the fatigue resistance of DP800 steel. The first step is a scale-bridging simulation approach that employs macroscopic simulations of the deep drawing process to extract the boundary conditions for the sRVE simulations with Damask. Subsequent cyclic loading is then simulated. 50 sRVE are computed for each load path to estimate fatigue resistance based on a Fatigue Indicator Parameter. The results indicate that fatigue resistance increases with increasing deformation-induced strain hardening. Additionally, a positive correlation between the martensitic ligament structures and fatigue resistance was observed.

**Introduction**

To achieve Europe's ambitious climate goals [1], it is increasingly necessary to design and construct structural components more efficiently. This requires continuous improvement of manufacturing processes, such as forming processes, in conjunction with the materials used. Several studies have already demonstrated that a well targeted and damage controlled forming process leads to improved component performance and can ultimately unlock lightweight potential [2,3]. In deep drawing, a critical sheet metal forming process, key parameters for modified process control include variables such as blank holder force, drawing velocity, friction, and the radius of the drawing ring [4,5].

To achieve a precise selection and optimization of materials for different manufacturing processes, it is essential to have a fundamental understanding of materials and their microstructures. Simulation methods are often used to achieve this understanding. In microstructure simulation, numerical crystal plasticity (CP) models are frequently used in conjunction with statistically Representative Volume Elements (sRVE), for example for the modelling of complex materials in [6,7]. Bargmann et al. [8] provide a comprehensive overview of the various methods for generating and applying sRVE. They offer a significant advantage in that they allow for the independent investigation of individual microstructure components, which can be experimentally challenging and often impractical [9]. Pütz et al. [10] quantified the influence of martensitic band structures in DP800 steel on the initiation of martensite fracture

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

<footnote n="1">¹Institute of Metal Forming, RWTH Aachen University, Intzestraße 10, 52072 Aachen, Germany</footnote>

<footnote n="2">²Manufacturing Technology Institute MTI of RWTH Aachen University, Campus-Boulevard 30, 52074 Aachen, Germany</footnote>

<footnote n="3">³Fraunhofer Institute for Production Technology IPT, Steinbachstraße 17, 52074 Aachen, Germany</footnote>

<footnote>ᵃniklas.fehlemann@ibf.rwth-aachen.de</footnote>
under bending loads using statistically representative volume elements (sRVE). In addition, the study found that the influence is depending on the phase contrast.

When optimizing component design, it is important to consider not only the forming process but also the load boundary conditions of the component application in the next step. Many parts and components are subjected to cyclic loads, and as a result, almost 90% of all component failures can be attributed to fatigue phenomena [11]. For accurate modelling of fatigue, the representation of the underlying physical mechanisms is crucial, which therefore makes crystal plasticity models paired with sRVE a method of choice [12]. Experimental investigations have shown that fatigue cracks frequently occur at persistent slip bands (PSBs) [13,14]. Based on these underlying physical mechanisms, various so-called Fatigue Indicator Parameters (FIPs) can be derived and used for prediction of the fatigue lifetime based on simulative results [15].

In this paper we present a simulation concept for the evaluation of fatigue life. Within this concept, the effects of the microstructure of a material as well as those of the load path before fatigue can be considered. The fatigue resistance of differently deep-drawn parts, such as seat rails, is used to demonstrate this concept. The aim is to gain a better understanding of the interaction between the complex strain paths during the manufacture and component in-use with the microstructure.

## Materials and Methods

The approach presented in this study is based on three distinct steps. First, a geometric reconstruction of the material to generate sRVE. Secondly, an approach for transferring the load paths of the macroscopic forming simulation to the sRVE and third, a methodology for the material simulation and the subsequent estimation of the fatigue resistance.

<figure><img src="image_2015.png" /><figcaption>Figure 1: Exemplary sRVE with martensite band (oriented along x- and z-direction) used in this study. Left: Representation of the individual grains (the colors are only an indicator for the individual grains, they do not correspond to a specific orientation). Right: visualization of the different phases (martensite: blue, ferrite: red)</figcaption></figure>

*Material reconstruction.* The example material in this study is a commercially available dual-phase steel DP800. The material was already extensively characterized in various previous studies, especially the complex micromechanical behavior of the material [16–18]. The used sRVE from the material were created using the AI-supported sRVE-generator DRAGen [19,20], developed by the author's group. This generator can reconstruct the morphology of the microstructure with greater accuracy than conventional methods, including the distinct martensitic banded structures that frequently occur in dual-phase steel grades. For these steels in particular, an accurate
under bending loads using statistically representative volume elements (sRVE). In addition, the study found that the influence is depending on the phase contrast.

When optimizing component design, it is important to consider not only the forming process but also the load boundary conditions of the component application in the next step. Many parts and components are subjected to cyclic loads, and as a result, almost 90% of all component failures can be attributed to fatigue phenomena [11]. For accurate modelling of fatigue, the representation of the underlying physical mechanisms is crucial, which therefore makes crystal plasticity models paired with sRVE a method of choice [12]. Experimental investigations have shown that fatigue cracks frequently occur at persistent slip bands (PSBs) [13,14]. Based on these underlying physical mechanisms, various so-called Fatigue Indicator Parameters (FIPs) can be derived and used for prediction of the fatigue lifetime based on simulative results [15].

In this paper we present a simulation concept for the evaluation of fatigue life. Within this concept, the effects of the microstructure of a material as well as those of the load path before fatigue can be considered. The fatigue resistance of differently deep-drawn parts, such as seat rails, is used to demonstrate this concept. The aim is to gain a better understanding of the interaction between the complex strain paths during the manufacture and component in-use with the microstructure.

## Materials and Methods

The approach presented in this study is based on three distinct steps. First, a geometric reconstruction of the material to generate sRVE. Secondly, an approach for transferring the load paths of the macroscopic forming simulation to the sRVE and third, a methodology for the material simulation and the subsequent estimation of the fatigue resistance.

<figure><img src="image_2010.png" /><figcaption>Figure 1: Exemplary sRVE with martensite band (oriented along x- and z-direction) used in this study. Left: Representation of the individual grains (the colors are only an indicator for the individual grains, they do not correspond to a specific orientation). Right: visualization of the different phases (martensite: blue, ferrite: red)</figcaption></figure>

*Material reconstruction.* The example material in this study is a commercially available dual-phase steel DP800. The material was already extensively characterized in various previous studies, especially the complex micromechanical behavior of the material [16–18]. The used sRVE from the material were created using the AI-supported sRVE-generator DRAGen [19,20], developed by the author's group. This generator can reconstruct the morphology of the microstructure with greater accuracy than conventional methods, including the distinct martensitic banded structures that frequently occur in dual-phase steel grades. For these steels in particular, an accurate
representation of the microstructure is important for a correct modelling of the mechanical properties [21]. For more details on the underlying generation mechanisms in the generator, the reader is referred to Henrich et al. [20]. A total of 50 statistically representative volume elements (sRVE) were generated, evenly distributed between 25 specimens featuring martensite bands and 25 without. To maintain statistical representativeness while simultaneously managing computational costs, a size of 20x20x20 µm was chosen. In this context, 1 µm was discretized into 2 points, resulting in an edge length of 40. The martensite content was specified at 30%. The ultimate martensite content of the sRVE was determined as $31.2 ± 2.3$ %, the slight differences arise from the discrete sRVE-generation process. The configuration of martensite bands adhered to data from a referenced source [22]. Figure 1 provides a representative depiction of an sRVE featuring a martensite band. The left side of the figure shows the grain structure, the right side shows the distribution of the phases. For the sRVE simulation, the spectral solver of the open source framework DAMASK [23] is used.

**Macroscopic simulation approach.** To seamlessly transpose the boundary conditions of a forming process onto the sRVE in the subsequent step, a suitable scale bridging modeling approach becomes necessary. As an illustrative application, deep drawing with varied drawing die radii $r_D$

<figure><img src="image_2016.png" /><figcaption>Figure 2: Schematic structure of the deep drawing process. Figure adapted from [4].</figcaption></figure>

is considered, as detailed in [4] and depicted in Figure 2. Different drawing ring radii of 3, 6 and 9 mm were selected as exemplary process variations. To enable the transfer of boundary conditions, a specific element in a region of interest within the macroscopic simulation model must be chosen. Depending on the goals of the investigation, this selection can be based on various criteria, such as the element with the highest equivalent plastic strain or von-Mises equivalent stress. Alternatively, a specific location on the component may serve as a reference element for extracting the boundary conditions. In this study, the criterion chosen was the maximum equivalent plastic strain at the end of the load path, aimed at examining the impact of distinct strain hardening behavior on the fatigue properties. It is noteworthy that the element is not constant, it varies with different die radius. Due to the different simulation codes employed for macroscopic (Abaqus/CAE) and microscopic (Damask) simulations, the conventional Abaqus internal submodelling method for transferring boundary conditions is ruled out. Consequently, an explicit user-defined subroutine (VUMAT) is utilized to perform the elasto-plastic simulation of deep drawing, generating the deformation gradient as an output. The nine components of the deformation gradient serve as a load case for the Damask simulation, comprehensively characterize the load path. As a simplification, the entire forming process is assumed to occur at a constant strain rate of 0.001 1/s, owing to the unavailability of strain rate-dependent material parameters for the sRVE-simulations. Nonetheless, given the uniform speed at which all deep drawing tests are conducted, the impact of this simplification is anticipated to be minimal.

**Microscopic simulation approach.** The complete load path for the simulation is as follows: First, the sRVE is loaded according to the load path of the deep drawing. The sRVE is then relieved so the homogenized macroscopic stress drops to zero and only residual stresses remain in the
representation of the microstructure is important for a correct modelling of the mechanical properties [21]. For more details on the underlying generation mechanisms in the generator, the reader is referred to Henrich et al. [20]. A total of 50 statistically representative volume elements (sRVE) were generated, evenly distributed between 25 specimens featuring martensite bands and 25 without. To maintain statistical representativeness while simultaneously managing computational costs, a size of 20x20x20 µm was chosen. In this context, 1 µm was discretized into 2 points, resulting in an edge length of 40. The martensite content was specified at 30%. The ultimate martensite content of the sRVE was determined as $31.2 \pm 2.3 \%$, the slight differences arise from the discrete sRVE-generation process. The configuration of martensite bands adhered to data from a referenced source [22]. Figure 1 provides a representative depiction of an sRVE featuring a martensite band. The left side of the figure shows the grain structure, the right side shows the distribution of the phases. For the sRVE simulation, the spectral solver of the open source framework DAMASK [23] is used.

**Macroscopic simulation approach.** To seamlessly transpose the boundary conditions of a forming process onto the sRVE in the subsequent step, a suitable scale bridging modeling approach becomes necessary. As an illustrative application, deep drawing with varied drawing die radii $r_D$

<figure><img src="image_2011.png" /><figcaption>**Figure 2: Schematic structure of the deep drawing process.** Figure adapted from [4].</figcaption></figure>

is considered, as detailed in [4] and depicted in Figure 2. Different drawing ring radii of 3, 6 and 9 mm were selected as exemplary process variations. To enable the transfer of boundary conditions, a specific element in a region of interest within the macroscopic simulation model must be chosen. Depending on the goals of the investigation, this selection can be based on various criteria, such as the element with the highest equivalent plastic strain or von-Mises equivalent stress. Alternatively, a specific location on the component may serve as a reference element for extracting the boundary conditions. In this study, the criterion chosen was the maximum equivalent plastic strain at the end of the load path, aimed at examining the impact of distinct strain hardening behavior on the fatigue properties. It is noteworthy that the element is not constant, it varies with different die radius. Due to the different simulation codes employed for macroscopic (Abaqus/CAE) and microscopic (Damask) simulations, the conventional Abaqus internal submodelling method for transferring boundary conditions is ruled out. Consequently, an explicit user-defined subroutine (VUMAT) is utilized to perform the elasto-plastic simulation of deep drawing, generating the deformation gradient as an output. The nine components of the deformation gradient serve as a load case for the Damask simulation, comprehensively characterize the load path. As a simplification, the entire forming process is assumed to occur at a constant strain rate of 0.001 1/s, owing to the unavailability of strain rate-dependent material parameters for the sRVE-simulations. Nonetheless, given the uniform speed at which all deep drawing tests are conducted, the impact of this simplification is anticipated to be minimal.

**Microscopic simulation approach.** The complete load path for the simulation is as follows: First, the sRVE is loaded according to the load path of the deep drawing. The sRVE is then relieved so the homogenized macroscopic stress drops to zero and only residual stresses remain in the
microstructure. Finally, the sRVE is cyclically loaded in a strain-controlled manner with an amplitude of 0.5% for 10 cycles.

A crystal plasticity model is used to map the mechanical behavior of the ferrite. Given the unavailability of valid information about the crystallographic orientations in the martensite due to the unresolved very small substructures in the used sRVE, the martensite is assumed to exhibit isotropic plasticity. The model employed for this purpose is also accessible in Damask, and its formulation closely resembles the crystal plasticity model utilized [23]. The parameters for DP800 were taken from a prior publication and are listed in Table 1:

<table><caption>Table 1: Material model parameters used for ferrite and martensite for DP800.</caption><thead><tr><th>Phase</th><th>N<sub>Slip</sub></th><th>γ<sub>slip</sub> (-)</th><th>a (-)</th><th>n (-)</th><th>ξ<sup>0</sup><sub>Slip</sub> (MPa)</th><th>ξ<sup>inf</sup><sub>Slip</sub> (MPa)</th><th>h<sup>0</sup><sub>Slip-Slip</sub> (MPa)</th></tr></thead><tbody><tr><td>Ferrite</td><td>24</td><td>0.001</td><td>4.5</td><td>5</td><td>150</td><td>400</td><td>625</td></tr><tr><td>Martensite</td><td>--</td><td>0.001</td><td>1.25</td><td>30</td><td>250</td><td>750</td><td>2500</td></tr></tbody></table>

The fatigue modeling approach adheres to the methodologies outlined in Gillner et al. and Henrich et al. [12,24]. The critical parameter in this context is the plastic slip $γ_{sl}$ accumulated during cyclic loading, serving as an indicator of dislocation slip and the occurrence of persistent slip bands. Within the framework of crystal plasticity, $γ_{sl}$ is computed based on the plastic velocity gradient ($L_p$).

$$γ_{sl} = \int_{0}^{t} \sqrt{\frac{2}{3} L_p : L_p} dt \qquad (1)$$

The velocity gradient is directly linked to the deformation gradient:

$$\dot{F}_p = L_p F_p \qquad (2)$$

In this context, $F_p$ denotes the plastic deformation gradient, and $\dot{F}_p$ represents its derivative w.r.t. the time. In the crystal plasticity model embedded in Damask, the plastic deformation gradient is integrated into the multiplicative decomposition of the overall deformation gradient.

$$F = F_e F_i F_p \qquad (3)$$

Here, $F_e$ represents the elastic and $F_p$ denotes the plastic component of the deformation gradient. $F_i$ characterizes inelastic strains arising from factors such as damage or phase transformations. The overall deformation gradient $F$ is prescribed as a load path.

After the simulation, the grain-averaged plastic slip is extracted using a python-based postprocessing routine and the maximum is taken over each sRVE. Here, it should be noted that the plastic slip accumulation is only calculated for the cyclic loading, the accumulation during the forming process is not considered:

$$γ_{sl,max}^{rve} = \max \left( \frac{1}{N_{pts}^{gr}} \sum_{i=0}^{N_{pts}^{gr}} γ_{sl} \right) \quad | \quad gr = 1...N_{gr} \qquad (4)$$

The parameters are defined as follows: $N_{gr}$ is the number of grains per Representative Volume Element, $N_{pts}^{gr}$ is the number of points per grain, and $γ_{sl}$ represents the accumulated plastic slip per material point. The averaging process relies on the number of points, which is a common practice in other methods. Consequently, each sRVE yields on specific $γ_{sl,max}^{rve}$ value. To accommodate the
microstructure. Finally, the sRVE is cyclically loaded in a strain-controlled manner with an amplitude of 0.5% for 10 cycles.

A crystal plasticity model is used to map the mechanical behavior of the ferrite. Given the unavailability of valid information about the crystallographic orientations in the martensite due to the unresolved very small substructures in the used sRVE, the martensite is assumed to exhibit isotropic plasticity. The model employed for this purpose is also accessible in Damask, and its formulation closely resembles the crystal plasticity model utilized [23]. The parameters for DP800 were taken from a prior publication and are listed in Table 1:

<table><caption>Table 1: Material model parameters used for ferrite and martensite for DP800.</caption><thead><tr><th>Phase</th><th>N<sub>Slip</sub></th><th>γ<sub>slip</sub> (-)</th><th>a (-)</th><th>n (-)</th><th>ξ<sup>0</sup><sub>Slip</sub> (MPa)</th><th>ξ<sup>inf</sup><sub>Slip</sub> (MPa)</th><th>h<sup>0</sup><sub>Slip-Slip</sub> (MPa)</th></tr></thead><tbody><tr><td>Ferrite</td><td>24</td><td>0.001</td><td>4.5</td><td>5</td><td>150</td><td>400</td><td>625</td></tr><tr><td>Martensite</td><td>--</td><td>0.001</td><td>1.25</td><td>30</td><td>250</td><td>750</td><td>2500</td></tr></tbody></table>

The fatigue modeling approach adheres to the methodologies outlined in Gillner et al. and Henrich et al. [12,24]. The critical parameter in this context is the plastic slip $γ_{sl}$ accumulated during cyclic loading, serving as an indicator of dislocation slip and the occurrence of persistent slip bands. Within the framework of crystal plasticity, $γ_{sl}$ is computed based on the plastic velocity gradient ($L_p$).

$$γ_{sl} = \int_{0}^{t} \sqrt{\frac{2}{3} L_p : L_p} dt \qquad (1)$$

The velocity gradient is directly linked to the deformation gradient:

$$\dot{F}_p = L_p F_p \qquad (2)$$

In this context, $F_p$ denotes the plastic deformation gradient, and $\dot{F}_p$ represents its derivative w.r.t. the time. In the crystal plasticity model embedded in Damask, the plastic deformation gradient is integrated into the multiplicative decomposition of the overall deformation gradient.

$$F = F_e F_i F_p \qquad (3)$$

Here, $F_e$ represents the elastic and $F_p$ denotes the plastic component of the deformation gradient. $F_i$ characterizes inelastic strains arising from factors such as damage or phase transformations. The overall deformation gradient $F$ is prescribed as a load path.

After the simulation, the grain-averaged plastic slip is extracted using a python-based postprocessing routine and the maximum is taken over each sRVE. Here, it should be noted that the plastic slip accumulation is only calculated for the cyclic loading, the accumulation during the forming process is not considered:

$$γ_{sl,max}^{rve} = \max\left(\frac{1}{N_{pts}^{gr}} \sum_{i=0}^{N_{pts}^{gr}} γ_{sl}\right) \quad | \quad gr = 1...N_{gr} \qquad (4)$$

The parameters are defined as follows: $N_{gr}$ is the number of grains per Representative Volume Element, $N_{pts}^{gr}$ is the number of points per grain, and $γ_{sl}$ represents the accumulated plastic slip per material point. The averaging process relies on the number of points, which is a common practice in other methods. Consequently, each sRVE yields on specific $γ_{sl,max}^{rve}$ value. To accommodate the
scatter in the simulations, the values for each load case are depicted in a histogram, and a right-sided Gumbel distribution is fitted to this data. This approach originates in the weakest link theory [25]:

$$F(x; \mu, \beta) = e^{-e^{-\frac{(x-\mu)}{\beta}}} \qquad (5)$$

The parameters $\mu$ and $\beta$ characterize the location and scale of the probability density function (pdf). As per Henrich et al. [24], the product of $\mu$ and $\beta$ results in a Fatigue Indicator Parameter (FIP), offering insights into the material's resistance to fatigue. In general, a high FIP indicates a high probability of fatigue failure and vice versa.

$$FIP = \mu * \beta \qquad (6)$$

It is crucial to emphasize that this indicator is purely qualitative and does not offer insights into the actual service life of the material or component. Its suitability lies in comparing various load cases and/or material conditions.

## Results

<figure><img src="image_2017.png" /><figcaption>Figure 3: Deformation gradients for all load cases considered. The solid line shows the F11 component (x-direction). The F33 component (grey dashed line) is 1 for all 3 load paths, which means no deformation in the z-direction.</figcaption></figure>

Figure 3 displays the evolution data of the deformation gradients extracted from the macroscopic simulation. The solid lines depict the F11 component of the deformation gradient (x-direction), while the dashed grey line represents the F33 component (z), which is consistently 1.0 for all three load cases, indicating a plane strain state in the respective element. Notably, significant differences in final deformation are evident among various drawing die radii. Specifically, at a radius of 9 mm, the F11 deformation gradient at the end of the load case is around 1.22, contrasting with 1.03 and 1.06 for radii 6 and 9, respectively. The tension-compression load change is more pronounced at radius 3 than at radius 6, and it is observed that shear components of the deformation gradient are relatively low. Therefore, the shear components are not depicted here.
scatter in the simulations, the values for each load case are depicted in a histogram, and a right-sided Gumbel distribution is fitted to this data. This approach originates in the weakest link theory [25]:

$$F(x; \mu, \beta) = e^{-e^{-\frac{(x-\mu)}{\beta}}} \qquad (5)$$

The parameters $\mu$ and $\beta$ characterize the location and scale of the probability density function (pdf). As per Henrich et al. [24], the product of $\mu$ and $\beta$ results in a Fatigue Indicator Parameter (FIP), offering insights into the material's resistance to fatigue. In general, a high FIP indicates a high probability of fatigue failure and vice versa.

$$FIP = \mu * \beta \qquad (6)$$

It is crucial to emphasize that this indicator is purely qualitative and does not offer insights into the actual service life of the material or component. Its suitability lies in comparing various load cases and/or material conditions.

## Results

<figure><img src="image_2012.png" /><figcaption>Figure 3: Deformation gradients for all load cases considered. The solid line shows the F11 component (x-direction). The F33 component (grey dashed line) is 1 for all 3 load paths, which means no deformation in the z-direction.</figcaption></figure>

Figure 3 displays the evolution data of the deformation gradients extracted from the macroscopic simulation. The solid lines depict the F11 component of the deformation gradient (x-direction), while the dashed grey line represents the F33 component (z), which is consistently 1.0 for all three load cases, indicating a plane strain state in the respective element. Notably, significant differences in final deformation are evident among various drawing die radii. Specifically, at a radius of 9 mm, the F11 deformation gradient at the end of the load case is around 1.22, contrasting with 1.03 and 1.06 for radii 6 and 9, respectively. The tension-compression load change is more pronounced at radius 3 than at radius 6, and it is observed that shear components of the deformation gradient are relatively low. Therefore, the shear components are not depicted here.
<figure><img src="image_2018.png" /><figcaption>Figure 4: Histogram and fitted Gumbel distribution (orange) for drawing die radius of 3mm.</figcaption></figure>

<figure><img src="image_2019.png" /><figcaption>Figure 5: Histogram and fitted Gumbel distribution (orange) for drawing die radius of 6mm.</figcaption></figure>

<figure><img src="image_2020.png" /><figcaption>Figure 6: Histogram and fitted Gumbel distribution (orange) for drawing die radius of 9mm.</figcaption></figure>

The different load paths consequently result in distinct fatigue behavior. The figures 4 - 6 depict histograms and fitted Gumbel distributions for the 50 sRVE for all three considered load cases.
<figure><img src="image_2013.png" /><figcaption>Figure 4: Histogram and fitted Gumbel distribution (orange) for drawing die radius of 3mm.</figcaption></figure>

<figure><img src="image_2014.png" /><figcaption>Figure 5: Histogram and fitted Gumbel distribution (orange) for drawing die radius of 6mm.</figcaption></figure>

<figure><img src="image_2015.png" /><figcaption>Figure 6: Histogram and fitted Gumbel distribution (orange) for drawing die radius of 9mm.</figcaption></figure>

The different load paths consequently result in distinct fatigue behavior. The figures 4 - 6 depict histograms and fitted Gumbel distributions for the 50 sRVE for all three considered load cases.
Alongside the visual representations in the histograms, the computed fatigue indicators are presented in Table 2. A FIP was calculated for all 50 sRVE and a separate FIP for a subset with and without martensite bands from 25 sRVE each.

<table><caption>Table 2: Fatigue indicator parameter and location and scale parameter of the Gumbel distribution for all three considered load cases.</caption><thead><tr><th></th><th></th><th>FIP (-)</th><th>μ (-)</th><th>β (-)</th><th>No. sRVE</th></tr></thead><tbody><tr><td rowspan="3">Radius 9mm</td><td>All</td><td>0.0010</td><td>0.0887</td><td>0.0117</td><td>50</td></tr><tr><td>Bands</td><td>0.0009</td><td>0.0880</td><td>0.0101</td><td>25</td></tr><tr><td>No Bands</td><td>0.0012</td><td>0.0895</td><td>0.0130</td><td>25</td></tr><tr><td rowspan="3">Radius 6mm</td><td>All</td><td>0.0024</td><td>0.1316</td><td>0.0182</td><td>50</td></tr><tr><td>Bands</td><td>0.0017</td><td>0.1330</td><td>0.0128</td><td>25</td></tr><tr><td>No Bands</td><td>0.0028</td><td>0.1308</td><td>0.0210</td><td>25</td></tr><tr><td rowspan="3">Radius 3mm</td><td>All</td><td>0.0015</td><td>0.1046</td><td>0.0144</td><td>50</td></tr><tr><td>Bands</td><td>0.0012</td><td>0.1030</td><td>0.0114</td><td>25</td></tr><tr><td>No Bands</td><td>0.0018</td><td>0.1068</td><td>0.0166</td><td>25</td></tr></tbody></table>

These data indicate that the Fatigue Indicator Parameter for a drawing die radius of 9 mm is significantly lower than that for 3 mm and 6 mm, with a difference of approximately 2.4 times to radius 6 and about 1.5 times to radius 3. This discrepancy is attributed to both a higher average accumulated strain (μ) and increased scatter (β). The differences in the FIP directly correlate with the final equivalent plastic strain of the different load paths, as illustrated in Figure 7. A higher final level of strain hardening leads to a lower FIP. It is also noteworthy that, for all three load cases, the FIP for the subset without martensite bands is higher than for the subset with martensite bands. This is primarily attributed to the higher scatter in the former.

<figure><img src="image_2021.png" /><figcaption>Figure 7: Correlation between Fatigue Indicator Parameter and global equivalent plastic strain, averaged over all 50 sRVE. The FIP tends to decrease with higher strain hardening.</figcaption></figure>

Discussion

The results presented in the preceding section can be analyzed through several stages. Firstly, a distinct correlation between strain hardening and the FIP becomes evident, where higher levels of strain hardening result in a lower FIP and therefore a greater resistance to fatigue. This association is explicable by the fact that, at a higher level of strain hardening, a higher portion of the deformation at the grain level is elastic, thereby reducing the accumulation of plastic slip.
Alongside the visual representations in the histograms, the computed fatigue indicators are presented in Table 2. A FIP was calculated for all 50 sRVE and a separate FIP for a subset with and without martensite bands from 25 sRVE each.

<table><caption>Table 2: Fatigue indicator parameter and location and scale parameter of the Gumbel distribution for all three considered load cases.</caption><thead><tr><th></th><th></th><th>FIP (-)</th><th>μ (-)</th><th>β (-)</th><th>No. sRVE</th></tr></thead><tbody><tr><td rowspan="3">Radius 9mm</td><td>All</td><td>0.0010</td><td>0.0887</td><td>0.0117</td><td>50</td></tr><tr><td>Bands</td><td>0.0009</td><td>0.0880</td><td>0.0101</td><td>25</td></tr><tr><td>No Bands</td><td>0.0012</td><td>0.0895</td><td>0.0130</td><td>25</td></tr><tr><td rowspan="3">Radius 6mm</td><td>All</td><td>0.0024</td><td>0.1316</td><td>0.0182</td><td>50</td></tr><tr><td>Bands</td><td>0.0017</td><td>0.1330</td><td>0.0128</td><td>25</td></tr><tr><td>No Bands</td><td>0.0028</td><td>0.1308</td><td>0.0210</td><td>25</td></tr><tr><td rowspan="3">Radius 3mm</td><td>All</td><td>0.0015</td><td>0.1046</td><td>0.0144</td><td>50</td></tr><tr><td>Bands</td><td>0.0012</td><td>0.1030</td><td>0.0114</td><td>25</td></tr><tr><td>No Bands</td><td>0.0018</td><td>0.1068</td><td>0.0166</td><td>25</td></tr></tbody></table>

These data indicate that the Fatigue Indicator Parameter for a drawing die radius of 9 mm is significantly lower than that for 3 mm and 6 mm, with a difference of approximately 2.4 times to radius 6 and about 1.5 times to radius 3. This discrepancy is attributed to both a higher average accumulated strain (μ) and increased scatter (β). The differences in the FIP directly correlate with the final equivalent plastic strain of the different load paths, as illustrated in Figure 7. A higher final level of strain hardening leads to a lower FIP. It is also noteworthy that, for all three load cases, the FIP for the subset without martensite bands is higher than for the subset with martensite bands. This is primarily attributed to the higher scatter in the former.

<figure><img src="image_2016.png" /><figcaption>Figure 7: Correlation between Fatigue Indicator Parameter and global equivalent plastic strain, averaged over all 50 sRVE. The FIP tends to decrease with higher strain hardening.</figcaption></figure>

Discussion

The results presented in the preceding section can be analyzed through several stages. Firstly, a distinct correlation between strain hardening and the FIP becomes evident, where higher levels of strain hardening result in a lower FIP and therefore a greater resistance to fatigue. This association is explicable by the fact that, at a higher level of strain hardening, a higher portion of the deformation at the grain level is elastic, thereby reducing the accumulation of plastic slip.
However, it is observed that both the mean strain and the scattering increase with decreasing strain hardening. Experimental studies indicate that scattering tends to increase with an extended fatigue life, which deviates from the observed behavior here. Further investigations are required to elucidate the reasons behind the opposing trends exhibited by the location and scale parameters. It should also be noted that preforming should not be advantageous for parts that are already pre-fatigued. With such parts, non-critical fatigue cracks may become critical. However, as forming from an initial state takes place here, this phenomenon may not be relevant.

A final aspect to discuss is the distinction between the sRVE with and without martensite bands. It becomes evident that the absence of martensite bands primarily results in a higher $\beta$ parameter, indicating higher scatter. Since a higher proportion of the strain is absorbed by the martensite in the presence of martensitic band structures [10], this can lead to a more homogeneous distribution of strain in the ferrite and therefore less scattering. Therefore, although martensite bands favor the occurrence of ductile damage, they can potentially be beneficial in terms of the fatigue resistance.

Lastly, the analyses presented here does not account for the influence of microscopic damage. In general, micropores may develop - primarily through the fracturing of the hard martensite in the case of the current DP800 - that can also serve as initiation points for fatigue cracks. Given that the distinct load paths are expected to result in varying degrees of damage depending on the drawing die radius [4], the consideration of damage is essential for a more generalized assessment of fatigue resistance. Addressing this aspect can be achieved within the framework of representative volume element simulations, for instance, using suitable phase field models [26], even though acquiring the necessary data for the calibration of these models poses a considerable challenge.

## Conclusion

A new simulation concept was presented that uses sRVE simulations and a scalebridging approach to evaluate the effects of microstructure and load path on fatigue resistance. The application of this approach to the fatigue resistance of deep-drawn components made of DP800 with different process parameters revealed a strong dependence of the fatigue resistance on the degree of hardening. It has been demonstrated that the presence of martensitic band structures can positively affect fatigue resistance. Although variables such as deformation-induced damage have not yet been considered, there is a concept available that can be utilized in the design of new materials and manufacturing processes.

## Acknowledgements

This research was funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation; Projectnumber 278868966 – TRR 188; Damage Controlled Forming Processes, subproject B05 and A06). Simulations were performed with computing resources granted by RWTH Aachen University under project rwth0744.

## References

[1] European Commission, European Green Deal. Available at https://commission.europa.eu/strategy-and-policy/priorities-2019-2024/european-green-deal_en.

[2] G. Hirt, A.E. Tekkaya, T. Clausmeyer, and J. Lohmar, Potential and status of damage controlled forming processes, Production Engineering 14 (2020), pp. 1–4. https://doi.org/10.1007/s11740-019-00948-6

[3] A.E. Tekkaya, N. Ben Khalifa, O. Hering, R. Meya, S. Myslicki, and F. Walther, Forming-induced damage and its effects on product properties, CIRP Annals 66 (2017), pp. 281–284. https://doi.org/10.1016/j.cirp.2017.04.113

[4] M. Müller, T. Herrig, and T. Bergs, Numerical analysis of the load paths and the resulting damage evolution during deep drawing of dual-phase steel, IOP Conference Series: Materials
However, it is observed that both the mean strain and the scattering increase with decreasing strain hardening. Experimental studies indicate that scattering tends to increase with an extended fatigue life, which deviates from the observed behavior here. Further investigations are required to elucidate the reasons behind the opposing trends exhibited by the location and scale parameters. It should also be noted that preforming should not be advantageous for parts that are already pre-fatigued. With such parts, non-critical fatigue cracks may become critical. However, as forming from an initial state takes place here, this phenomenon may not be relevant.

A final aspect to discuss is the distinction between the sRVE with and without martensite bands. It becomes evident that the absence of martensite bands primarily results in a higher $\beta$ parameter, indicating higher scatter. Since a higher proportion of the strain is absorbed by the martensite in the presence of martensitic band structures [10], this can lead to a more homogeneous distribution of strain in the ferrite and therefore less scattering. Therefore, although martensite bands favor the occurrence of ductile damage, they can potentially be beneficial in terms of the fatigue resistance.

Lastly, the analyses presented here does not account for the influence of microscopic damage. In general, micropores may develop - primarily through the fracturing of the hard martensite in the case of the current DP800 - that can also serve as initiation points for fatigue cracks. Given that the distinct load paths are expected to result in varying degrees of damage depending on the drawing die radius [4], the consideration of damage is essential for a more generalized assessment of fatigue resistance. Addressing this aspect can be achieved within the framework of representative volume element simulations, for instance, using suitable phase field models [26], even though acquiring the necessary data for the calibration of these models poses a considerable challenge.

## Conclusion

A new simulation concept was presented that uses sRVE simulations and a scalebridging approach to evaluate the effects of microstructure and load path on fatigue resistance. The application of this approach to the fatigue resistance of deep-drawn components made of DP800 with different process parameters revealed a strong dependence of the fatigue resistance on the degree of hardening. It has been demonstrated that the presence of martensitic band structures can positively affect fatigue resistance. Although variables such as deformation-induced damage have not yet been considered, there is a concept available that can be utilized in the design of new materials and manufacturing processes.

## Acknowledgements

This research was funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation; Projectnumber 278868966 – TRR 188; Damage Controlled Forming Processes, subproject B05 and A06). Simulations were performed with computing resources granted by RWTH Aachen University under project rwth0744.

## References

[1] European Commission, European Green Deal. Available at https://commission.europa.eu/strategy-and-policy/priorities-2019-2024/european-green-deal_en.

[2] G. Hirt, A.E. Tekkaya, T. Clausmeyer, and J. Lohmar, Potential and status of damage controlled forming processes, Production Engineering 14 (2020), pp. 1–4. https://doi.org/10.1007/s11740-019-00948-6

[3] A.E. Tekkaya, N. Ben Khalifa, O. Hering, R. Meya, S. Myslicki, and F. Walther, Forming-induced damage and its effects on product properties, CIRP Annals 66 (2017), pp. 281–284. https://doi.org/10.1016/j.cirp.2017.04.113

[4] M. Müller, T. Herrig, and T. Bergs, Numerical analysis of the load paths and the resulting damage evolution during deep drawing of dual-phase steel, IOP Conference Series: Materials
Science and Engineering 1284 (2023), p. 12006. https://doi.org/10.1088/1757-899X/1284/1/012006

[5] P.R. Tiwari, A. Rathore, and M.G. Bodkhe, Factors affecting the deep drawing process – A review, Materials Today: Proceedings 56 (2022), pp. 2902–2908. https://doi.org/10.1016/j.matpr.2021.10.189

[6] F. Qayyum, S. Guk, M. Schmidtchen, R. Kawalla, and U. Prahl, Modeling the Local Deformation and Transformation Behavior of Cast X8CrMnNi16-6-6 TRIP Steel and 10% Mg-PSZ Composite Using a Continuum Mechanics-Based Crystal Plasticity Model, Crystals 10 (2020), p. 221. https://doi.org/10.3390/cryst10030221

[7] D. Reimann, K. Nidadavolu, H. ul Hassan, N. Vajragupta, T. Glasmachers, P. Junker, and A. Hartmaier, Modeling Macroscopic Material Behavior With Machine Learning Algorithms Trained by Micromechanical Simulations, Frontiers in Materials 6 (2019). https://doi.org/10.3389/fmats.2019.00181

[8] S. Bargmann, B. Klusemann, J. Markmann, J.E. Schnabel, K. Schneider, C. Soyarslan, and J. Wilmers, Generation of 3D representative volume elements for heterogeneous materials: A review, Progress in Materials Science 96 (2018), pp. 322–384. https://doi.org/10.1016/j.pmatsci.2018.02.003

[9] J.P.M. Hoefnagels, B.G. Vossen, and C.C. Tasan, Electron Micrographic Digital Image Correlation: Method Optimization and Microstructural banding Case Study, in: Application of Imaging Techniques to Mechanics of Materials and Structures, Volume 4, T. Proulx, ed., Conference Proceedings of the Society for Experimental Mechanics Series. Springer New York, New York, NY, 2013, pp. 71–77. https://doi.org/10.1007/978-1-4419-9796-8_10

[10] F. Pütz, N. Fehlemann, V. Göksu, M. Henrich, M. Könemann, and S. Münstermann, A data driven computational microstructure analysis on the influence of martensite banding on damage in DP-steels, Computational Materials Science 218 (2023), p. 111903. https://doi.org/10.1016/j.commatsci.2022.111903

[11] R.I. Stephens, A. Fatemi, R.R. Stephens, and H.O. Fuchs, Metal fatigue in engineering, 2nd edn. Wiley, c2001, New York, 2001.

[12] K. Gillner and S. Münstermann, Numerically predicted high cycle fatigue properties through representative volume elements of the microstructure, International Journal of Fatigue 105 (2017), pp. 219–234. https://doi.org/10.1016/j.ijfatigue.2017.09.002

[13] J. Polák and J. Man, Mechanisms of extrusion and intrusion formation in fatigued crystalline materials, Materials Science and Engineering: A 596 (2014), pp. 15–24. https://doi.org/10.1016/j.msea.2013.12.005

[14] J. Polák, J. Man, T. Vystavěl, and M. Petrenec, The shape of extrusions and intrusions and initiation of stage I fatigue cracks, Materials Science and Engineering: A 517 (2009), pp. 204–211. https://doi.org/10.1016/j.msea.2009.03.070

[15] N. Sayer, E. Natkowski, P. Sonnweber-Ribic, and S. Münstermann, A novel microscale fatigue failure indicator considering plastic irreversibility for microstructure-based lifetime simulation, International Journal of Fatigue 163 (2022), p. 107115. https://doi.org/10.1016/j.ijfatigue.2022.107115

[16] R. Meya, C. Kusche, C. Löbbe, T. Al-Samman, S. Korte-Kerzel, and A. Tekkaya, Global and High-Resolution Damage Quantification in Dual-Phase Steel Bending Samples with Varying Stress States, Metals 9 (2019), p. 319. https://doi.org/10.3390/met9030319
Science and Engineering 1284 (2023), p. 12006. https://doi.org/10.1088/1757-899X/1284/1/012006

[5] P.R. Tiwari, A. Rathore, and M.G. Bodkhe, Factors affecting the deep drawing process – A review, Materials Today: Proceedings 56 (2022), pp. 2902–2908. https://doi.org/10.1016/j.matpr.2021.10.189

[6] F. Qayyum, S. Guk, M. Schmidtchen, R. Kawalla, and U. Prahl, Modeling the Local Deformation and Transformation Behavior of Cast X8CrMnNi16-6-6 TRIP Steel and 10% Mg-PSZ Composite Using a Continuum Mechanics-Based Crystal Plasticity Model, Crystals 10 (2020), p. 221. https://doi.org/10.3390/cryst10030221

[7] D. Reimann, K. Nidadavolu, H. ul Hassan, N. Vajragupta, T. Glasmachers, P. Junker, and A. Hartmaier, Modeling Macroscopic Material Behavior With Machine Learning Algorithms Trained by Micromechanical Simulations, Frontiers in Materials 6 (2019). https://doi.org/10.3389/fmats.2019.00181

[8] S. Bargmann, B. Klusemann, J. Markmann, J.E. Schnabel, K. Schneider, C. Soyarslan, and J. Wilmers, Generation of 3D representative volume elements for heterogeneous materials: A review, Progress in Materials Science 96 (2018), pp. 322–384. https://doi.org/10.1016/j.pmatsci.2018.02.003

[9] J.P.M. Hoefnagels, B.G. Vossen, and C.C. Tasan, Electron Micrographic Digital Image Correlation: Method Optimization and Microstructural banding Case Study, in: Application of Imaging Techniques to Mechanics of Materials and Structures, Volume 4, T. Proulx, ed., Conference Proceedings of the Society for Experimental Mechanics Series. Springer New York, New York, NY, 2013, pp. 71–77. https://doi.org/10.1007/978-1-4419-9796-8_10

[10] F. Pütz, N. Fehlemann, V. Göksu, M. Henrich, M. Könemann, and S. Münstermann, A data driven computational microstructure analysis on the influence of martensite banding on damage in DP-steels, Computational Materials Science 218 (2023), p. 111903. https://doi.org/10.1016/j.commatsci.2022.111903

[11] R.I. Stephens, A. Fatemi, R.R. Stephens, and H.O. Fuchs, Metal fatigue in engineering, 2nd edn. Wiley, c2001, New York, 2001.

[12] K. Gillner and S. Münstermann, Numerically predicted high cycle fatigue properties through representative volume elements of the microstructure, International Journal of Fatigue 105 (2017), pp. 219–234. https://doi.org/10.1016/j.ijfatigue.2017.09.002

[13] J. Polák and J. Man, Mechanisms of extrusion and intrusion formation in fatigued crystalline materials, Materials Science and Engineering: A 596 (2014), pp. 15–24. https://doi.org/10.1016/j.msea.2013.12.005

[14] J. Polák, J. Man, T. Vystavěl, and M. Petrenec, The shape of extrusions and intrusions and initiation of stage I fatigue cracks, Materials Science and Engineering: A 517 (2009), pp. 204–211. https://doi.org/10.1016/j.msea.2009.03.070

[15] N. Sayer, E. Natkowski, P. Sonnweber-Ribic, and S. Münstermann, A novel microscale fatigue failure indicator considering plastic irreversibility for microstructure-based lifetime simulation, International Journal of Fatigue 163 (2022), p. 107115. https://doi.org/10.1016/j.ijfatigue.2022.107115

[16] R. Meya, C. Kusche, C. Löbbe, T. Al-Samman, S. Korte-Kerzel, and A. Tekkaya, Global and High-Resolution Damage Quantification in Dual-Phase Steel Bending Samples with Varying Stress States, Metals 9 (2019), p. 319. https://doi.org/10.3390/met9030319
[17] C. Tian, D. Ponge, L. Christiansen, and C. Kirchlechner, On the mechanical heterogeneity in dual phase steel grades: Activation of slip systems and deformation of martensite in DP800, Acta Materialia 183 (2020), pp. 274–284. https://doi.org/10.1016/j.actamat.2019.11.002

[18] M.A. Wollenweber, S. Medghalchi, L.R. Guimarães, N. Lohrey, C.F. Kusche, U. Kerzel, T. Al-Samman, and S. Korte-Kerzel, On the damage behaviour in dual-phase DP800 steel deformed in single and combined strain paths, Materials & Design 231 (2023), p. 112016. https://doi.org/10.1016/j.matdes.2023.112016

[19] N. Fehlemann, Y. Sparrer, F. Pütz, M. Konemann, and S. Münstermann, Influence of synthetically generated inclusions on the stress accumulation and concentration in X65 pipeline steel, IOP Conference Series: Materials Science and Engineering 1157 (2021), p. 12056. https://doi.org/10.1088/1757-899X/1157/1/012056

[20] M. Henrich, N. Fehlemann, F. Bexter, M. Neite, L. Kong, F. Shen, M. Könemann, M. Dölz, and S. Münstermann, DRAGen - a Deep Learning Supported sRVE Generator Framework for Complex Microstructure Models, Heliyon 9 (2023), e19003. https://doi.org/10.1016/j.heliyon.2023.e19003

[21] D. Yan, C.C. Tasan, and D. Raabe, High resolution in situ mapping of microstrain and microstructure evolution reveals damage resistance criteria in dual phase steels, Acta Materialia 96 (2015), pp. 399–409. https://doi.org/10.1016/j.actamat.2015.05.038

[22] N. Fehlemann, A.L. Suarez Aguilera, S. Sandfeld, F. Bexter, M. Neite, D. Lenz, M. Könemann, and S. Münstermann, Identification of Martensite Bands in Dual-Phase Steels: A Deep Learning Object Detection Approach Using Faster Region-Based-Convolutional Neural Network, steel research int. 94 (2023). https://doi.org/10.1002/srin.202200836

[23] F. Roters, M. Diehl, P. Shanthraj, P. Eisenlohr, C. Reuber, S.L. Wong, T. Maiti, A. Ebrahimi, T. Hochrainer, H.-O. Fabritius, S. Nikolov, M. Friák, N. Fujita, N. Grilli, K. Janssens, N. Jia, P. Kok, D. Ma, F. Meier, E. Werner, M. Stricker, D. Weygand, and D. Raabe, DAMASK – The Düsseldorf Advanced Material Simulation Kit for modeling multi-physics crystal plasticity, thermal, and damage phenomena from the single crystal up to the component scale, Computational Materials Science 158 (2019), pp. 420–478. https://doi.org/10.1016/j.commatsci.2018.04.030

[24] M. Henrich, M. Dölz, and S. Münstermann, Development of a Numerical Framework for Microstructure Sensitive Fatigue Life Investigations, in: Volume 2: Computer Technology & Bolted Joints; Design & Analysis, ASME 2023 Pressure Vessels & Piping Conference, Atlanta, Georgia, USA, Jul. 16-21, 2023. American Society of Mechanical Engineers, 2023. https://doi.org/10.1115/PVP2023-106419

[25] S. Beretta and Y. Murakami, STATISTICAL ANALYSIS OF DEFECTS FOR FATIGUE STRENGTH PREDICTION AND QUALITY CONTROL OF MATERIALS, Fatigue Fract Eng Mat Struct 21 (1998), pp. 1049–1065. https://doi.org/10.1046/j.1460-2695.1998.00104.x

[26] P. Shanthraj, L. Sharma, B. Svendsen, F. Roters, and D. Raabe, A phase field model for damage in elasto-viscoplastic materials, Computer Methods in Applied Mechanics and Engineering 312 (2016), pp. 167–185. https://doi.org/10.1016/j.cma.2016.05.006
[17] C. Tian, D. Ponge, L. Christiansen, and C. Kirchlechner, On the mechanical heterogeneity in dual phase steel grades: Activation of slip systems and deformation of martensite in DP800, Acta Materialia 183 (2020), pp. 274–284. https://doi.org/10.1016/j.actamat.2019.11.002

[18] M.A. Wollenweber, S. Medghalchi, L.R. Guimarães, N. Lohrey, C.F. Kusche, U. Kerzel, T. Al-Samman, and S. Korte-Kerzel, On the damage behaviour in dual-phase DP800 steel deformed in single and combined strain paths, Materials & Design 231 (2023), p. 112016. https://doi.org/10.1016/j.matdes.2023.112016

[19] N. Fehlemann, Y. Sparrer, F. Pütz, M. Konemann, and S. Münstermann, Influence of synthetically generated inclusions on the stress accumulation and concentration in X65 pipeline steel, IOP Conference Series: Materials Science and Engineering 1157 (2021), p. 12056. https://doi.org/10.1088/1757-899X/1157/1/012056

[20] M. Henrich, N. Fehlemann, F. Bexter, M. Neite, L. Kong, F. Shen, M. Könemann, M. Dölz, and S. Münstermann, DRAGen - a Deep Learning Supported sRVE Generator Framework for Complex Microstructure Models, Heliyon 9 (2023), e19003. https://doi.org/10.1016/j.heliyon.2023.e19003

[21] D. Yan, C.C. Tasan, and D. Raabe, High resolution in situ mapping of microstrain and microstructure evolution reveals damage resistance criteria in dual phase steels, Acta Materialia 96 (2015), pp. 399–409. https://doi.org/10.1016/j.actamat.2015.05.038

[22] N. Fehlemann, A.L. Suarez Aguilera, S. Sandfeld, F. Bexter, M. Neite, D. Lenz, M. Könemann, and S. Münstermann, Identification of Martensite Bands in Dual-Phase Steels: A Deep Learning Object Detection Approach Using Faster Region-Based-Convolutional Neural Network, steel research int. 94 (2023). https://doi.org/10.1002/srin.202200836

[23] F. Roters, M. Diehl, P. Shanthraj, P. Eisenlohr, C. Reuber, S.L. Wong, T. Maiti, A. Ebrahimi, T. Hochrainer, H.-O. Fabritius, S. Nikolov, M. Friák, N. Fujita, N. Grilli, K. Janssens, N. Jia, P. Kok, D. Ma, F. Meier, E. Werner, M. Stricker, D. Weygand, and D. Raabe, DAMASK – The Düsseldorf Advanced Material Simulation Kit for modeling multi-physics crystal plasticity, thermal, and damage phenomena from the single crystal up to the component scale, Computational Materials Science 158 (2019), pp. 420–478. https://doi.org/10.1016/j.commatsci.2018.04.030

[24] M. Henrich, M. Dölz, and S. Münstermann, Development of a Numerical Framework for Microstructure Sensitive Fatigue Life Investigations, in: Volume 2: Computer Technology & Bolted Joints; Design & Analysis, ASME 2023 Pressure Vessels & Piping Conference, Atlanta, Georgia, USA, Jul. 16-21, 2023. American Society of Mechanical Engineers, 2023. https://doi.org/10.1115/PVP2023-106419

[25] S. Beretta and Y. Murakami, STATISTICAL ANALYSIS OF DEFECTS FOR FATIGUE STRENGTH PREDICTION AND QUALITY CONTROL OF MATERIALS, Fatigue Fract Eng Mat Struct 21 (1998), pp. 1049–1065. https://doi.org/10.1046/j.1460-2695.1998.00104.x

[26] P. Shanthraj, L. Sharma, B. Svendsen, F. Roters, and D. Raabe, A phase field model for damage in elasto-viscoplastic materials, Computer Methods in Applied Mechanics and Engineering 312 (2016), pp. 167–185. https://doi.org/10.1016/j.cma.2016.05.006
# A new view on the solution of rate-independent crystal plasticity finite element models

NIJHUIS Björn¹,ᵃ*, PERDAHCIOĞLU Semih¹,ᵇ and VAN DEN BOOGAARD Ton¹,ᶜ

¹Chair of Nonlinear Solid Mechanics, University of Twente, Enschede, The Netherlands

ᵃb.nijhuis@utwente.nl, ᵇe.s.perdahcioglu@utwente.nl, ᶜa.h.vandenboogaard@utwente.nl

**Keywords:** Crystal Plasticity, Stress Update Algorithm, Active Set, Fixed-Point Iterations

**Abstract.** The crystal plasticity finite element method (CP-FEM) readily enables microstructure-based material modelling by relating macroscopic plastic deformation to dislocation slip on crystal slip systems. Rate-independent CP models provide physically accurate solutions by allowing slip only if the resolved shear stress on a slip system equals the critical resolved shear stress. However, computing the amount of slip for such models remains challenging. This work proposes a novel stable and efficient stress-update algorithm based on fixed-point iterations. These iterations trace the hypersurfaces that describe the slip state for which individual slip system's yield functions are zero, until all slip system hypersurfaces intersect. This simultaneously provides the set of active slip systems and the slip on these systems, avoiding the need for an iterative active set search algorithm without inducing spurious slip on systems on which the shear stress is below the critical resolved shear stress.

**Introduction**

The crystal plasticity finite element method (CP-FEM) is a powerful tool for incorporating microstructural effects into predictions of macroscopic mechanical behaviour of polycrystalline materials. CP-FEM models plastic deformation as the result of the combined slip activity on the material's crystallographic slip systems, requiring the amount of slip to be computed given the applied loading.

For deformations at low temperature, rate-independent CP models provide physically accurate solutions by allowing slip only if the resolved shear stress (RSS) on a slip system equals the critical resolved shear stress (CRSS) [1]. Active and inactive slip systems must then be distinguished, but the selection of the active set is nontrivial as it is not uniquely defined at start of an increment [2]. Therefore, stress update algorithms for rate-independent CP-FEM may start with a guess for the active set and compute the slip on the systems in that set, after which the active set is updated and slip computations are restarted if the resulting stress state requires slip systems to be (de)activated [3]. Such iterative search procedures clearly limit the computational efficiency of rate-independent CP models. As such, various authors introduced more sophisticated algorithms to initialize the active set [4-6].

The complexity of determining an active set led to the development of multiple CP formulations that do not require an active set at all. One approach is to relate the resolved and critical resolved shear stresses of all systems in a single yield function that defines a yield surface with rounded corners [7]. It is, however, difficult to incorporate effects of slip system interactions, e.g. latent hardening, in such a model [8].

A popular method that also does not require an active set, but does easily allow to include slip system interactions, is to employ viscoplastic regularisation and relate the slip activity of a slip system to the ratio between its RSS and CRSS [9]. This, however, renders the model rate-dependent and results in small amounts of slip on systems on which the RSS does not exceed the CRSS. An increase in plastic deformation may then cause extra hardening. Choosing model

<figure><img src="image_2022.png" /></figure>
# A new view on the solution of rate-independent crystal plasticity finite element models

NIJHUIS Björn¹,ᵃ*, PERDAHCIOĞLU Semih¹,ᵇ and VAN DEN BOOGAARD Ton¹,ᶜ

¹Chair of Nonlinear Solid Mechanics, University of Twente, Enschede, The Netherlands

ᵃb.nijhuis@utwente.nl, ᵇe.s.perdahcioglu@utwente.nl, ᶜa.h.vandenboogaard@utwente.nl

**Keywords:** Crystal Plasticity, Stress Update Algorithm, Active Set, Fixed-Point Iterations

**Abstract.** The crystal plasticity finite element method (CP-FEM) readily enables microstructure-based material modelling by relating macroscopic plastic deformation to dislocation slip on crystal slip systems. Rate-independent CP models provide physically accurate solutions by allowing slip only if the resolved shear stress on a slip system equals the critical resolved shear stress. However, computing the amount of slip for such models remains challenging. This work proposes a novel stable and efficient stress-update algorithm based on fixed-point iterations. These iterations trace the hypersurfaces that describe the slip state for which individual slip system's yield functions are zero, until all slip system hypersurfaces intersect. This simultaneously provides the set of active slip systems and the slip on these systems, avoiding the need for an iterative active set search algorithm without inducing spurious slip on systems on which the shear stress is below the critical resolved shear stress.

**Introduction**

The crystal plasticity finite element method (CP-FEM) is a powerful tool for incorporating microstructural effects into predictions of macroscopic mechanical behaviour of polycrystalline materials. CP-FEM models plastic deformation as the result of the combined slip activity on the material's crystallographic slip systems, requiring the amount of slip to be computed given the applied loading.

For deformations at low temperature, rate-independent CP models provide physically accurate solutions by allowing slip only if the resolved shear stress (RSS) on a slip system equals the critical resolved shear stress (CRSS) [1]. Active and inactive slip systems must then be distinguished, but the selection of the active set is nontrivial as it is not uniquely defined at start of an increment [2]. Therefore, stress update algorithms for rate-independent CP-FEM may start with a guess for the active set and compute the slip on the systems in that set, after which the active set is updated and slip computations are restarted if the resulting stress state requires slip systems to be (de)activated [3]. Such iterative search procedures clearly limit the computational efficiency of rate-independent CP models. As such, various authors introduced more sophisticated algorithms to initialize the active set [4-6].

The complexity of determining an active set led to the development of multiple CP formulations that do not require an active set at all. One approach is to relate the resolved and critical resolved shear stresses of all systems in a single yield function that defines a yield surface with rounded corners [7]. It is, however, difficult to incorporate effects of slip system interactions, e.g. latent hardening, in such a model [8].

A popular method that also does not require an active set, but does easily allow to include slip system interactions, is to employ viscoplastic regularisation and relate the slip activity of a slip system to the ratio between its RSS and CRSS [9]. This, however, renders the model rate-dependent and results in small amounts of slip on systems on which the RSS does not exceed the CRSS. An increase in plastic deformation may then cause extra hardening. Choosing model

<figure><img src="image_2017.png" /></figure>
parameters that result in a lower strain-rate sensitivity reduces the spurious slips, but also results in an ill-conditioned system of equations [10].

Recently proposed stress update algorithms based on the interior point method also avoid the active set search by rendering all slip systems active [11, 12]. These algorithms result in a truly rate-independent response, but alike viscoplastic models, induce slip even if the RSS on a system is below its CRSS.

This work aims at developing a novel stress-update algorithm for rate-independent CP models that does not require an active set search procedure and does not induce spurious slip on systems on which the RSS is lower than the CRSS. We first formulate the rate-independent crystal plasticity model, after which we propose a stress-update algorithm based on fixed-point iterations. Then, we demonstrate the feasibility of this algorithm and compare its stability with that of a conventional Newton-Raphson based stress update algorithm for rate-independent CP.

## Formulation of the crystal plasticity model

We take as the starting point for the formulation of the CP model the multiplicative decomposition of the deformation gradient tensor $\mathbf{F}$ into an elastic ($\mathbf{F}_e$) and a plastic part ($\mathbf{F}_p$) as

$$ \mathbf{F} = \mathbf{F}_e \cdot \mathbf{F}_p. \qquad (1) $$

Here, $\mathbf{F}_p$ defines an intermediate configuration that results from plastic slip caused by dislocation motion and in which the crystal lattice is undistorted. $\mathbf{F}_e$ accounts for the rotation and elastic stretch of the lattice and relates this intermediate configuration to the current configuration. In the following, an overset hat (') is used to denote symbols that are mapped from the intermediate to the current configuration. The total velocity gradient $\mathbf{L} = \dot{\mathbf{F}} \cdot \mathbf{F}^{-1}$ can now be additively decomposed into its elastic ($\hat{\mathbf{L}}_e$) and plastic ($\hat{\mathbf{L}}_p$) components as

$$ \mathbf{L} = \hat{\mathbf{L}}_e + \hat{\mathbf{L}}_p, \qquad (2) $$

while the plastic deformation gradient is given by a summation of the plastic slip rates on the crystal's $n_s$ slip systems as

$$ \hat{\mathbf{L}}_p = \sum_{\alpha=1}^{n_s} \dot{\mathbf{y}}^{(\alpha)} \hat{\mathbf{s}}^{(\alpha)} \otimes \hat{\mathbf{m}}^{(\alpha)}. \qquad (3) $$

Here, $\dot{\mathbf{y}}^{(\alpha)}$ is the plastic slip rate on a slip system $\alpha$ that is defined by the unit slip direction vector $\hat{\mathbf{s}}^{(\alpha)}$ and unit slip plane normal vector $\hat{\mathbf{m}}^{(\alpha)}$. In this work, we consider a face centred cubic (FCC) lattice.

In a rate-independent CP model, slip on a system $\alpha$ is allowed only if the resolved shear stress $\tau_r^{(\alpha)}$ on that system equals the critical resolved shear stress $\tau_f^{(\alpha)}$. The resolved shear stress $\tau_r^{(\alpha)}$ is calculated by projecting the stress $\sigma$ as

$$ \tau_r^{(\alpha)} = \sigma: \hat{\mathbf{P}}^{(\alpha)} = \sigma: \frac{1}{2} (\hat{\mathbf{s}}^{(\alpha)} \otimes \hat{\mathbf{m}}^{(\alpha)} + \hat{\mathbf{m}}^{(\alpha)} \otimes \hat{\mathbf{s}}^{(\alpha)}), \qquad (4) $$

where $\hat{\mathbf{P}}^{(\alpha)}$ denotes the symmetric Schmid tensor. We compute the critical resolved shear stress $\tau_f^{(\alpha)}$ using a Taylor type hardening law as

$$ \tau_f^{(\alpha)} = \tau_0 + Gb \sqrt{\sum_{\beta=1}^{n_s} Q^{(\alpha\beta)} \rho^{(\beta)}}. \qquad (5) $$
parameters that result in a lower strain-rate sensitivity reduces the spurious slips, but also results in an ill-conditioned system of equations [10].

Recently proposed stress update algorithms based on the interior point method also avoid the active set search by rendering all slip systems active [11, 12]. These algorithms result in a truly rate-independent response, but alike viscoplastic models, induce slip even if the RSS on a system is below its CRSS.

This work aims at developing a novel stress-update algorithm for rate-independent CP models that does not require an active set search procedure and does not induce spurious slip on systems on which the RSS is lower than the CRSS. We first formulate the rate-independent crystal plasticity model, after which we propose a stress-update algorithm based on fixed-point iterations. Then, we demonstrate the feasibility of this algorithm and compare its stability with that of a conventional Newton-Raphson based stress update algorithm for rate-independent CP.

## Formulation of the crystal plasticity model

We take as the starting point for the formulation of the CP model the multiplicative decomposition of the deformation gradient tensor $\mathbf{F}$ into an elastic ($\mathbf{F}_e$) and a plastic part ($\mathbf{F}_p$) as

$$ \mathbf{F} = \mathbf{F}_e \cdot \mathbf{F}_p. \qquad (1) $$

Here, $\mathbf{F}_p$ defines an intermediate configuration that results from plastic slip caused by dislocation motion and in which the crystal lattice is undistorted. $\mathbf{F}_e$ accounts for the rotation and elastic stretch of the lattice and relates this intermediate configuration to the current configuration. In the following, an overset hat (') is used to denote symbols that are mapped from the intermediate to the current configuration. The total velocity gradient $\mathbf{L} = \dot{\mathbf{F}} \cdot \mathbf{F}^{-1}$ can now be additively decomposed into its elastic ($\hat{\mathbf{L}}_e$) and plastic ($\hat{\mathbf{L}}_p$) components as

$$ \mathbf{L} = \hat{\mathbf{L}}_e + \hat{\mathbf{L}}_p, \qquad (2) $$

while the plastic deformation gradient is given by a summation of the plastic slip rates on the crystal's $n_s$ slip systems as

$$ \hat{\mathbf{L}}_p = \sum_{\alpha=1}^{n_s} \dot{\mathbf{y}}^{(\alpha)} \hat{\mathbf{s}}^{(\alpha)} \otimes \hat{\mathbf{m}}^{(\alpha)}. \qquad (3) $$

Here, $\dot{\mathbf{y}}^{(\alpha)}$ is the plastic slip rate on a slip system $\alpha$ that is defined by the unit slip direction vector $\hat{\mathbf{s}}^{(\alpha)}$ and unit slip plane normal vector $\hat{\mathbf{m}}^{(\alpha)}$. In this work, we consider a face centred cubic (FCC) lattice.

In a rate-independent CP model, slip on a system $\alpha$ is allowed only if the resolved shear stress $\tau_r^{(\alpha)}$ on that system equals the critical resolved shear stress $\tau_f^{(\alpha)}$. The resolved shear stress $\tau_r^{(\alpha)}$ is calculated by projecting the stress $\sigma$ as

$$ \tau_r^{(\alpha)} = \sigma: \hat{\mathbf{P}}^{(\alpha)} = \sigma: \frac{1}{2} (\hat{\mathbf{s}}^{(\alpha)} \otimes \hat{\mathbf{m}}^{(\alpha)} + \hat{\mathbf{m}}^{(\alpha)} \otimes \hat{\mathbf{s}}^{(\alpha)}), \qquad (4) $$

where $\hat{\mathbf{P}}^{(\alpha)}$ denotes the symmetric Schmid tensor. We compute the critical resolved shear stress $\tau_f^{(\alpha)}$ using a Taylor type hardening law as

$$ \tau_f^{(\alpha)} = \tau_0 + Gb \sqrt{\sum_{\beta=1}^{n_s} Q^{(\alpha\beta)} \rho^{(\beta)}}. \qquad (5) $$
Here, $\tau_0$ represents a constant lattice friction, $G$ is the shear modulus and $b$ is the length of the Burgers vector. The component $Q^{(\alpha\beta)}$ of the interaction matrix $Q$ determines the hardening of slip system $\alpha$ as a result of an increase in dislocation density $\rho^{(\beta)}$ on slip system $\beta$. The components of $Q$ are governed by six independent types of slip system interactions [13], as listed in Table 1, for which the coefficients are taken from results of dislocation dynamics analysis conducted by Kubin et al. [13]. The relation between the dislocation density and the accumulated plastic slip $\gamma^{(\alpha)} = \int_{t_0}^{t_{\text{end}}} |\dot{\gamma}^{(\alpha)}| dt$ on a slip system $\alpha$ is modelled with the phenomenological constitutive model given by Becker [14]:

$$ \rho^{(\alpha)} = \rho_{\infty} \left[ 1 - \left( 1 - \frac{\rho_0}{\rho_{\infty}} \right) \exp\left(-\frac{\gamma^{\alpha}}{\gamma_{\infty}}\right) \right]. \qquad (6) $$

Here, $\rho_0$, $\rho_{\infty}$ and $\gamma_{\infty}$ are the initial dislocation density, saturation dislocation density and saturation slip, respectively. The values for the constants in Eqs. 5 and 6 are provided in Table 2.

<table><caption>Table 1: Coefficients for slip system interactions in an FCC lattice [13]</caption><thead><tr><th>Self</th><th>Coplanar</th><th>Collinear</th><th>Orthogonal</th><th>Glissile</th><th>Sessile</th></tr></thead><tbody><tr><td>0.122</td><td>0.122</td><td>0.625</td><td>0.070</td><td>0.137</td><td>0.122</td></tr></tbody></table>

<table><caption>Table 2: Model parameters used in Eqs. 5 and 6 [12]</caption><thead><tr><th>E [GPa]</th><th>ν [-]</th><th>τ<sub>0</sub> [MPa]</th><th>b [mm]</th><th>ρ<sub>0</sub> [mm-2]</th><th>ρ<sub>∞</sub> [mm-2]</th><th>γ<sub>∞</sub> [-]</th></tr></thead><tbody><tr><td>72.0</td><td>0.3</td><td>18.0</td><td>2.86×10<sup>-7</sup></td><td>107</td><td>109</td><td>0.4</td></tr></tbody></table>

## Fixed-point iteration scheme for slip computations

Now, we will formulate a stress update algorithm that provides the increments of stress and plastic slip on the crystal's slip systems during a time increment from $t_n$ to $t_{n+1} = t_n + \Delta t$ as function of the total rate of deformation tensor $D$ (the symmetric part of $L$, Eq. 2). A trial stress state is given as

$$ \sigma_{n+1}^* = \sigma_n + \hat{C}_e : \Delta D. \qquad (7) $$

Here, $\hat{C}_e$ denotes the elasticity tensor of the lattice in the current configuration. If for at least one slip system the RSS resulting from the trial stress state exceeds the CRSS, the closest point projection (CPP) from the trial state onto the yield surface is computed. The nonlinear systems of equations that defines this projection can be derived using the principle of maximum dissipation. For a detailed treatment of this derivation, the reader is referred to [12]. It follows that for all active slip systems, which define the active set $A$, it has to hold that

$$ \mathbf{R}_{\sigma} = \hat{C}_{e}: \Delta \mathbf{D} - \Delta \boldsymbol{\sigma} - \hat{C}_{e}: \sum_{\alpha \in A} \hat{\mathbf{P}}^{(\alpha)} \Delta \gamma^{(\alpha)} = \mathbf{0}, \qquad (8) $$

$$ \phi^{(\alpha)} = \tau_{r}^{(\alpha)} - \tau_{f}^{(\alpha)} = 0 \quad \forall \alpha \in A, \qquad (9) $$

while for all systems that are not in the active set it has to hold that

$$ \phi^{(\alpha)} = \tau_{r}^{(\alpha)} - \tau_{f}^{(\alpha)} < 0 \quad \forall \alpha \in N. \qquad (10) $$

Here, $N = T \setminus A$ denotes the set of inactive systems, with $T$ being the set of all $n_s$ slip systems. In Eq. 8, it is used that $\Delta\gamma^{(\alpha)} = \dot{\gamma}^{(\alpha)}\Delta t$ and $\Delta D = D\Delta t$.
Here, $\tau_0$ represents a constant lattice friction, $G$ is the shear modulus and $b$ is the length of the Burgers vector. The component $Q^{(\alpha\beta)}$ of the interaction matrix $Q$ determines the hardening of slip system $\alpha$ as a result of an increase in dislocation density $\rho^{(\beta)}$ on slip system $\beta$. The components of $Q$ are governed by six independent types of slip system interactions [13], as listed in Table 1, for which the coefficients are taken from results of dislocation dynamics analysis conducted by Kubin et al. [13]. The relation between the dislocation density and the accumulated plastic slip $\gamma^{(\alpha)} = \int_{t_0}^{t_{\text{end}}} |\dot{\gamma}^{(\alpha)}| dt$ on a slip system $\alpha$ is modelled with the phenomenological constitutive model given by Becker [14]:

$$ \rho^{(\alpha)} = \rho_{\infty} \left[ 1 - \left( 1 - \frac{\rho_0}{\rho_{\infty}} \right) \exp\left(-\frac{\gamma^{\alpha}}{\gamma_{\infty}}\right) \right]. \qquad (6) $$

Here, $\rho_0$, $\rho_{\infty}$ and $\gamma_{\infty}$ are the initial dislocation density, saturation dislocation density and saturation slip, respectively. The values for the constants in Eqs. 5 and 6 are provided in Table 2.

<table><caption>Table 1: Coefficients for slip system interactions in an FCC lattice [13]</caption><thead><tr><th>Self</th><th>Coplanar</th><th>Collinear</th><th>Orthogonal</th><th>Glissile</th><th>Sessile</th></tr></thead><tbody><tr><td>0.122</td><td>0.122</td><td>0.625</td><td>0.070</td><td>0.137</td><td>0.122</td></tr></tbody></table>

<table><caption>Table 2: Model parameters used in Eqs. 5 and 6 [12]</caption><thead><tr><th>E [GPa]</th><th>ν [-]</th><th>τ<sub>0</sub> [MPa]</th><th>b [mm]</th><th>ρ<sub>0</sub> [mm-2]</th><th>ρ<sub>∞</sub> [mm-2]</th><th>γ<sub>∞</sub> [-]</th></tr></thead><tbody><tr><td>72.0</td><td>0.3</td><td>18.0</td><td>2.86×10<sup>-7</sup></td><td>107</td><td>109</td><td>0.4</td></tr></tbody></table>

## Fixed-point iteration scheme for slip computations

Now, we will formulate a stress update algorithm that provides the increments of stress and plastic slip on the crystal's slip systems during a time increment from $t_n$ to $t_{n+1} = t_n + \Delta t$ as function of the total rate of deformation tensor $D$ (the symmetric part of $L$, Eq. 2). A trial stress state is given as

$$ \sigma_{n+1}^* = \sigma_n + \hat{C}_e: \Delta D. \qquad (7) $$

Here, $\hat{C}_e$ denotes the elasticity tensor of the lattice in the current configuration. If for at least one slip system the RSS resulting from the trial stress state exceeds the CRSS, the closest point projection (CPP) from the trial state onto the yield surface is computed. The nonlinear systems of equations that defines this projection can be derived using the principle of maximum dissipation. For a detailed treatment of this derivation, the reader is referred to [12]. It follows that for all active slip systems, which define the active set $A$, it has to hold that

$$ \mathbf{R}_{\sigma} = \hat{C}_{e}: \Delta \mathbf{D} - \Delta \boldsymbol{\sigma} - \hat{C}_{e}: \sum_{\alpha \in A} \hat{\mathbf{P}}^{(\alpha)} \Delta \gamma^{(\alpha)} = \mathbf{0}, \qquad (8) $$

$$ \phi^{(\alpha)} = \tau_{r}^{(\alpha)} - \tau_{f}^{(\alpha)} = 0 \quad \forall \alpha \in A, \qquad (9) $$

while for all systems that are not in the active set it has to hold that

$$ \phi^{(\alpha)} = \tau_{r}^{(\alpha)} - \tau_{f}^{(\alpha)} < 0 \quad \forall \alpha \in N. \qquad (10) $$

Here, $N = T \setminus A$ denotes the set of inactive systems, with $T$ being the set of all $n_s$ slip systems. In Eq. 8, it is used that $\Delta\gamma^{(\alpha)} = \dot{\gamma}^{(\alpha)}\Delta t$ and $\Delta D = D\Delta t$.
A common approach to solve the system given by Eqs. 8 and 9 is to linearise these equations and iteratively update the increments of plastic slip on the active systems and the stress using the Newton-Raphson (NR) algorithm [2, 3, 10]. This, clearly, requires the active set A to be established. If the conditions in Eqs. 9 and 10 cannot be satisfied for a given active set, it has to be updated and the NR iterations have to be restarted.

To avoid having to explicitly initialise and update the active set during a time increment, we propose to combine Eqs. 9 and 10 into a single condition that has to hold for each slip system:

$$ \tilde{\phi}^{(\alpha)} = \tau_{\mathrm{r}}^{(\alpha)} - \tilde{\tau}_{\mathrm{f}}^{(\alpha)} = 0 \quad \forall \alpha \in T, \qquad (11) $$

for which a modified flow stress is defined as

$$ \tilde{\tau}_{\mathrm{f}}^{(\alpha)} = \begin{cases} \tau_{\mathrm{r}}^{(\alpha)} & \text{if } \gamma^{(\alpha)} = 0 \land |\tau_{\mathrm{R}}^{(\alpha)}(\Delta\mathbf{D}, \Delta\mathbf{\gamma})| < \tau_{\mathrm{f}}^{(\alpha)}(\mathbf{\gamma}_n, \Delta\mathbf{\gamma}) \\ \mathrm{sign}(\gamma^{(\alpha)}) \tau_{\mathrm{f}}^{(\alpha)}(\mathbf{\gamma}_n, \Delta\mathbf{\gamma}) & \text{otherwise} \end{cases}. \qquad (12) $$

The condition in Eq. 11 can be rewritten in Mandel notation [15] as

$$ \tilde{\phi}^{(\alpha)}(\Delta\mathbf{D}, \mathbf{\gamma}_n, \Delta\mathbf{\gamma}) = \hat{\mathbf{P}}^{(\alpha)}{}^{\mathrm{T}} \underline{\underline{\sigma}}_{n+1}^{*}(\Delta\mathbf{D}) - \hat{\mathbf{P}}^{(\alpha)}{}^{\mathrm{T}} \hat{\mathcal{C}}_{\mathrm{e}} \hat{\mathbf{P}} \Delta\mathbf{\gamma} - \tilde{\tau}_{\mathrm{f}}^{(\alpha)}(\mathbf{\gamma}_n, \Delta\mathbf{\gamma}) = 0 \quad \forall \alpha \in T, \qquad (13) $$

using Eqs. 4 and 8 and noting that $\sigma_{n+1} = \sigma_n + \Delta\sigma$. In Eq. 13 and the remainder of the text, all underlined symbols denote matrix and vector representations of tensors in Mandel notation. Column $\alpha$ of the matrix $\hat{\mathbf{P}}$ corresponds to the Mandel vector $\hat{\mathbf{P}}^{(\alpha)}$ of the symmetric Schmid tensor $\hat{\mathbf{P}}^{(\alpha)}$ [12]. The vector $\mathbf{\gamma}_n$ contains the accumulated slips of all slip systems.

The solution $\Delta\mathbf{\gamma}$ to Eq. 13 is located at an intersection of the $n_s$ hypersurfaces $\tilde{\phi}^{(\alpha)}(\Delta\mathbf{\gamma}) = 0$. Each of these hypersurfaces defines where $\Delta\mathbf{\gamma}$ is such that the linear function

$$ r^{(\alpha)}(\Delta\mathbf{\gamma}) = \hat{\mathbf{P}}^{(\alpha)}{}^{\mathrm{T}} \hat{\mathcal{C}}_{\mathrm{e}} \hat{\mathbf{P}} \Delta\mathbf{\gamma} \qquad (14) $$

crosses the nonlinear function

$$ g^{(\alpha)}(\Delta\mathbf{D}, \mathbf{\gamma}_n, \Delta\mathbf{\gamma}) = \hat{\mathbf{P}}^{(\alpha)}{}^{\mathrm{T}} \underline{\underline{\sigma}}_{n+1}^{*}(\Delta\mathbf{D}) - \tilde{\tau}_{\mathrm{f}}^{(\alpha)}(\mathbf{\gamma}_n, \Delta\mathbf{\gamma}). \qquad (15) $$

Due to the nature of the functions in Eqs. 14 and 15 and the definition of $\tilde{\tau}_{\mathrm{f}}^{(\alpha)}$ in Eq. 12, this crossing always occurs at $|\Delta\gamma^{(\alpha)}| > 0$ if $|\tau_{\mathrm{r}}^{(\alpha)}| = \tilde{\tau}_{\mathrm{f}}^{(\alpha)}$, and at $\gamma^{(\alpha)} = 0$ if $|\tau_{\mathrm{r}}^{(\alpha)}| < \tilde{\tau}_{\mathrm{f}}^{(\alpha)}$. Figures 1a and 1b show this schematically for an exemplary single crystal with two slip systems that is subjected to a tensile stress of 300 MPa. The two slip systems are oriented at an angle of -15° and 60° with respect to the tensile axis, respectively. The CRSS on both systems evolves as given by Eqs. 5 and 6. For this example only, $Q^{(\alpha\beta)} = 1.4 - 0.4\delta_{\alpha\beta}$, where $\delta_{\alpha\beta}$ is the Kronecker delta function. As seen in Figure 1c, an intersection of the hypersurfaces $\tilde{\phi}^{(\alpha)} = 0$ defines both the active set and the slip activity on the (active) slip systems.

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

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

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

<figure><img src="image_2026.png" /></figure>
A common approach to solve the system given by Eqs. 8 and 9 is to linearise these equations and iteratively update the increments of plastic slip on the active systems and the stress using the Newton-Raphson (NR) algorithm [2, 3, 10]. This, clearly, requires the active set A to be established. If the conditions in Eqs. 9 and 10 cannot be satisfied for a given active set, it has to be updated and the NR iterations have to be restarted.

To avoid having to explicitly initialise and update the active set during a time increment, we propose to combine Eqs. 9 and 10 into a single condition that has to hold for each slip system:

$$ \tilde{\phi}^{(\alpha)} = \tau_{\mathrm{r}}^{(\alpha)} - \tilde{\tau}_{\mathrm{f}}^{(\alpha)} = 0 \quad \forall \alpha \in T, \qquad (11) $$

for which a modified flow stress is defined as

$$ \tilde{\tau}_{\mathrm{f}}^{(\alpha)} = \begin{cases} \tau_{\mathrm{r}}^{(\alpha)} & \text{if } \gamma^{(\alpha)} = 0 \land |\tau_{\mathrm{R}}^{(\alpha)}(\Delta\mathbf{D}, \Delta\mathbf{\gamma})| < \tau_{\mathrm{f}}^{(\alpha)}(\mathbf{\gamma}_n, \Delta\mathbf{\gamma}) \\ \mathrm{sign}(\gamma^{(\alpha)}) \tau_{\mathrm{f}}^{(\alpha)}(\mathbf{\gamma}_n, \Delta\mathbf{\gamma}) & \text{otherwise} \end{cases}. \qquad (12) $$

The condition in Eq. 11 can be rewritten in Mandel notation [15] as

$$ \tilde{\phi}^{(\alpha)}(\Delta\mathbf{D}, \mathbf{\gamma}_n, \Delta\mathbf{\gamma}) = \hat{\mathbf{P}}^{(\alpha)}{}^{\mathrm{T}} \underline{\mathbf{\sigma}}_{n+1}^{*}(\Delta\mathbf{D}) - \hat{\mathbf{P}}^{(\alpha)}{}^{\mathrm{T}} \hat{\mathcal{C}}_{\mathrm{e}} \hat{\mathbf{P}} \Delta\mathbf{\gamma} - \tilde{\tau}_{\mathrm{f}}^{(\alpha)}(\mathbf{\gamma}_n, \Delta\mathbf{\gamma}) = 0 \quad \forall \alpha \in T, \qquad (13) $$

using Eqs. 4 and 8 and noting that $\sigma_{n+1} = \sigma_n + \Delta\sigma$. In Eq. 13 and the remainder of the text, all underlined symbols denote matrix and vector representations of tensors in Mandel notation. Column $\alpha$ of the matrix $\hat{\mathbf{P}}$ corresponds to the Mandel vector $\hat{\mathbf{P}}^{(\alpha)}$ of the symmetric Schmid tensor $\hat{\mathbf{P}}^{(\alpha)}$ [12]. The vector $\mathbf{\gamma}_n$ contains the accumulated slips of all slip systems.

The solution $\Delta\mathbf{\gamma}$ to Eq. 13 is located at an intersection of the $n_s$ hypersurfaces $\tilde{\phi}^{(\alpha)}(\Delta\mathbf{\gamma}) = 0$. Each of these hypersurfaces defines where $\Delta\mathbf{\gamma}$ is such that the linear function

$$ r^{(\alpha)}(\Delta\mathbf{\gamma}) = \hat{\mathbf{P}}^{(\alpha)}{}^{\mathrm{T}} \hat{\mathcal{C}}_{\mathrm{e}} \hat{\mathbf{P}} \Delta\mathbf{\gamma} \qquad (14) $$

crosses the nonlinear function

$$ g^{(\alpha)}(\Delta\mathbf{D}, \mathbf{\gamma}_n, \Delta\mathbf{\gamma}) = \hat{\mathbf{P}}^{(\alpha)}{}^{\mathrm{T}} \underline{\mathbf{\sigma}}_{n+1}^{*}(\Delta\mathbf{D}) - \tilde{\tau}_{\mathrm{f}}^{(\alpha)}(\mathbf{\gamma}_n, \Delta\mathbf{\gamma}). \qquad (15) $$

Due to the nature of the functions in Eqs. 14 and 15 and the definition of $\tilde{\tau}_{\mathrm{f}}^{(\alpha)}$ in Eq. 12, this crossing always occurs at $|\Delta\gamma^{(\alpha)}| > 0$ if $|\tau_{\mathrm{r}}^{(\alpha)}| = \tilde{\tau}_{\mathrm{f}}^{(\alpha)}$, and at $\gamma^{(\alpha)} = 0$ if $|\tau_{\mathrm{r}}^{(\alpha)}| < \tilde{\tau}_{\mathrm{f}}^{(\alpha)}$. Figures 1a and 1b show this schematically for an exemplary single crystal with two slip systems that is subjected to a tensile stress of 300 MPa. The two slip systems are oriented at an angle of -15° and 60° with respect to the tensile axis, respectively. The CRSS on both systems evolves as given by Eqs. 5 and 6. For this example only, $Q^{(\alpha\beta)} = 1.4 - 0.4\delta_{\alpha\beta}$, where $\delta_{\alpha\beta}$ is the Kronecker delta function. As seen in Figure 1c, an intersection of the hypersurfaces $\tilde{\phi}^{(\alpha)} = 0$ defines both the active set and the slip activity on the (active) slip systems.

■ : $r^{(1)}$ (Eq. 14)
- : $g^{(1)}$ (Eq. 15)
: $\tilde{\phi}^{(1)} = 0$ (Eq. 13)

■ : $r^{(2)}$ (Eq. 14)
- : $g^{(2)}$ (Eq. 15)
: $\tilde{\phi}^{(2)} = 0$ (Eq. 13)

-/- : $\tilde{\phi}^{(1)} = 0 / \tilde{\phi}^{(2)} = 0$
□ : $\Delta\gamma_k$ without line search
○ : $\Delta\gamma_k$ with line search

<figure><img src="image_2018.png" /></figure>
<figure><img src="image_2027.png" /><figcaption>Figure 1: Visualisation of (a)-(b) the hypersurfaces $\tilde{\phi}^{(\alpha)} = 0$ for the model shown top right and (c) the use of these hypersurfaces to determine the active set and slip activity.</figcaption></figure>

Based on this observation, we developed a novel fixed-point iteration (FPI) scheme that traces the hypersurfaces of the individual systems until it finds the fixed-point at which all $n_s$ hypersurfaces intersect. Table 3 shows the pseudocode for the proposed FPI scheme. Figure 1c shows the progress of the fixed-point iterations to compute the slip on the two slip systems for the exemplary single crystal introduced above.

Given an estimate $\Delta\mathbf{y}_k$ for $\Delta\bar{\mathbf{y}}$, the FPI scheme performs two main steps. In the first step, $n_s$ one-dimensional root-finding problems are solved. Each solution provides the step length $t^{(\alpha)}$ for which the point $\Delta\mathbf{y}_{k+1}^{(\alpha)} = \Delta\mathbf{y}_k + t^{(\alpha)}\mathbf{d}^{(\alpha)}$ is located on the hypersurface $\tilde{\phi}^{(\alpha)} = 0$. The search direction vector $\mathbf{d}^{(\alpha)}$ is chosen as the normalised $\alpha$-th column of the symmetric matrix $\hat{\mathbf{P}}^T\hat{\mathcal{C}}_e\hat{\mathbf{P}}$. The results $\Delta\mathbf{y}_{k+1}^{(\alpha)}$ are combined to compute the direction $\Delta\Delta\mathbf{y}_{k+1}$ along which the slip increment will be updated by setting $\Delta\Delta\mathbf{y}_{k+1}^{(\alpha)} = t^{(\alpha)}\mathbf{d}^{(\alpha\alpha)}$. Here, $\mathbf{d}^{(\alpha\alpha)}$ is the $\alpha$-th component of the search direction vector $\mathbf{d}^{(\alpha)}$.

In the second step, a line search along $\Delta\Delta\mathbf{y}_{k+1}$ is performed and $\Delta\mathbf{y}_k$ is updated as $\Delta\mathbf{y}_{k+1} = \Delta\mathbf{y}_k + q\Delta\Delta\mathbf{y}_{k+1}$. The step length $q > 0$ is chosen such as to minimise the L2-norm of the residual functions $\tilde{\phi}^{(\alpha)}(\Delta\mathbf{y}_k + q\Delta\Delta\mathbf{y}_{k+1})$ of all currently active slip systems. The positive influence of this line search step on the convergence behaviour of the FPI scheme can be clearly observed in Figure 1c. When the line search step is enabled, the number of iterations to compute the slip activity of the exemplary single crystal reduces from 17 to only 2.

The one-dimensional root-finding and the scalar minimisation problem are solved by the general and standard algorithms given by Brent [16]. These algorithms require only function evaluations and simple bracketing operations and are therefore easily implemented and numerically cheap to evaluate.

The crystal lattice rotates during deformation. This requires to update $\hat{\mathbf{P}}^{(\alpha)}$ and $\hat{\mathcal{C}}_e$ during iterations. These updates are implemented as described in [12]. Allowing the lattice to rotate during iterations causes the fixed-point to move during the iterations. The influence thereof on the convergence of the FPI scheme will be analysed in the result section below.
<figure><img src="image_2019.png" /><figcaption>Figure 1: Visualisation of (a)-(b) the hypersurfaces $\tilde{\phi}^{(\alpha)} = 0$ for the model shown top right and (c) the use of these hypersurfaces to determine the active set and slip activity.</figcaption></figure>

Based on this observation, we developed a novel fixed-point iteration (FPI) scheme that traces the hypersurfaces of the individual systems until it finds the fixed-point at which all $n_s$ hypersurfaces intersect. Table 3 shows the pseudocode for the proposed FPI scheme. Figure 1c shows the progress of the fixed-point iterations to compute the slip on the two slip systems for the exemplary single crystal introduced above.

Given an estimate $\Delta\mathbf{y}_k$ for $\Delta\bar{\mathbf{y}}$, the FPI scheme performs two main steps. In the first step, $n_s$ one-dimensional root-finding problems are solved. Each solution provides the step length $t^{(\alpha)}$ for which the point $\Delta\mathbf{y}_{k+1}^{(\alpha)} = \Delta\mathbf{y}_k + t^{(\alpha)}\mathbf{d}^{(\alpha)}$ is located on the hypersurface $\tilde{\phi}^{(\alpha)} = 0$. The search direction vector $\mathbf{d}^{(\alpha)}$ is chosen as the normalised $\alpha$-th column of the symmetric matrix $\hat{\mathbf{P}}^T\hat{\mathcal{C}}_e\hat{\mathbf{P}}$. The results $\Delta\mathbf{y}_{k+1}^{(\alpha)}$ are combined to compute the direction $\Delta\Delta\mathbf{y}_{k+1}$ along which the slip increment will be updated by setting $\Delta\Delta\mathbf{y}_{k+1}^{(\alpha)} = t^{(\alpha)}\mathbf{d}^{(\alpha\alpha)}$. Here, $\mathbf{d}^{(\alpha\alpha)}$ is the $\alpha$-th component of the search direction vector $\mathbf{d}^{(\alpha)}$.

In the second step, a line search along $\Delta\Delta\mathbf{y}_{k+1}$ is performed and $\Delta\mathbf{y}_k$ is updated as $\Delta\mathbf{y}_{k+1} = \Delta\mathbf{y}_k + q\Delta\Delta\mathbf{y}_{k+1}$. The step length $q > 0$ is chosen such that to minimise the L2-norm of the residual functions $\tilde{\phi}^{(\alpha)}(\Delta\mathbf{y}_k + q\Delta\Delta\mathbf{y}_{k+1})$ of all currently active slip systems. The positive influence of this line search step on the convergence behaviour of the FPI scheme can be clearly observed in Figure 1c. When the line search step is enabled, the number of iterations to compute the slip activity of the exemplary single crystal reduces from 17 to only 2.

The one-dimensional root-finding and the scalar minimisation problem are solved by the general and standard algorithms given by Brent [16]. These algorithms require only function evaluations and simple bracketing operations and are therefore easily implemented and numerically cheap to evaluate.

The crystal lattice rotates during deformation. This requires to update $\hat{\mathbf{P}}^{(\alpha)}$ and $\hat{\mathcal{C}}_e$ during iterations. These updates are implemented as described in [12]. Allowing the lattice to rotate during iterations causes the fixed-point to move during the iterations. The influence thereof on the convergence of the FPI scheme will be analysed in the result section below.
After convergence of the FPI scheme, Eqs. 8 and 9 hold for all active slip systems. As such, the stress increment can be computed using Eq. 8 as

$$ \Delta \underline{\sigma} = \underline{\mathcal{C}}_e \Delta \underline{\mathbf{D}} - \underline{\mathcal{C}}_e \hat{\mathbf{P}} \Delta \underline{\mathbf{\gamma}}, \quad (16) $$

and the consistent algorithmic tangent $\underline{\mathcal{C}}_{ep}$ can be derived as [12]

$$ \underline{\mathcal{C}}_{ep} = (\underline{\mathcal{C}}_e^{-1} + \hat{\mathbf{P}}\mathbf{G}^{-1}\hat{\mathbf{P}}^{\mathrm{T}})^{-1}. \quad (17) $$

Here, $G^{(\alpha\beta)} = \partial\phi^{(\alpha)}/\partial\gamma^{(\beta)}$.

<table><caption>Table 3: Fixed-point iteration based stress update algorithm for rate-independent CP</caption>
  <tbody>
    <tr>
      <td colspan="2"><b>1. Initialisation:</b></td>
    </tr>
    <tr>
      <td>1</td>
      <td>k &larr; 0</td>
    </tr>
    <tr>
      <td>2</td>
      <td>&Delta;&mathbf{y}<sub>k</sub> &larr; 0</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Compute lattice rotation, update &Ccirc;<sub>e</sub> and Schmid tensors &hat{P}<sup>(&alpha;)</sup></td>
    </tr>
    <tr>
      <td>4</td>
      <td>Set search directions: d<sup>(&alpha;)</sup> &larr; normalised &alpha;-th column of &hat{P}<sup>T</sup>&Ccirc;<sub>e</sub> &hat{P}</td>
    </tr>
    <tr>
      <td>5</td>
      <td>&sigma;<sup>*</sup><sub>n+1</sub> &larr; &sigma;<sub>n</sub> + &Ccirc;<sub>e</sub> &Delta;D</td>
    </tr>
    <tr>
      <td>6</td>
      <td>r<sub>0</sub> = min(1.0, &radic;[&Sigma;<sub>&alpha; &isin; A<sub>k</sub></sub> (&vert; &hat{P}<sup>(&alpha;)</sup> &sigma;<sup>*</sup><sub>n+1</sub> &vert; - &tau;<sub>f</sub><sup>(&alpha;)</sup>)<sup>2</sup>), A<sub>k</sub> = {&alpha; &isin; T: &vert;&tau;<sub>r</sub><sup>(&alpha;)</sup>&vert; &gt; &tau;<sub>f</sub><sup>(&alpha;)</sup>}</td>
    </tr>
    <tr>
      <td>7</td>
      <td>r &larr; r<sub>0</sub></td>
    </tr>
    <tr>
      <td colspan="2"><b>2. FPI-iterations:</b></td>
    </tr>
    <tr>
      <td>8</td>
      <td>do while r &gt; 10<sup>-4</sup></td>
    </tr>
    <tr>
      <td>9</td>
      <td>k &larr; k + 1</td>
    </tr>
    <tr>
      <td>10</td>
      <td>for &alpha; = 1, 2, ..., n<sub>s</sub> do</td>
    </tr>
    <tr>
      <td>11</td>
      <td>find t<sup>(&alpha;)</sup>: solve &phi;&#771;<sup>(&alpha;)</sup>(&mathbf{y}<sub>n</sub>, &Delta;&mathbf{y}<sub>k-1</sub> + t<sup>(&alpha;)</sup>d<sup>(&alpha;)</sup>) = 0</td>
    </tr>
    <tr>
      <td>12</td>
      <td>&Delta;&Delta;&mathbf{y}<sub>k</sub><sup>(&alpha;)</sup> &larr; t<sup>(&alpha;)</sup>d<sup>(&alpha;&alpha;)</sup></td>
    </tr>
    <tr>
      <td>13</td>
      <td>for &alpha; = 1, 2, ..., n<sub>s</sub> do</td>
    </tr>
    <tr>
      <td>14</td>
      <td>if &Delta;&mathbf{y}<sub>k-1</sub><sup>(&alpha;)</sup> + &Delta;&Delta;&mathbf{y}<sub>k</sub><sup>(&alpha;)</sup> = 0:</td>
    </tr>
    <tr>
      <td>15</td>
      <td>&Delta;&mathbf{y}<sub>k-1</sub><sup>&alpha;</sup> &larr; 0</td>
    </tr>
    <tr>
      <td>16</td>
      <td>&Delta;&Delta;&mathbf{y}<sub>k</sub><sup>(&alpha;)</sup> &larr; 0</td>
    </tr>
    <tr>
      <td>17</td>
      <td>find q: minimize &radic;[&Sigma;<sub>&alpha; &isin; A<sub>k</sub></sub> &phi;&#771;<sup>(&alpha;)</sup>(&mathbf{y}<sub>n</sub>, &Delta;&mathbf{y}<sub>k-1</sub> + q&Delta;&Delta;&mathbf{y}<sub>k</sub>)<sup>2</sup>], A<sub>k</sub> = {&alpha; &isin; T: &vert;&tau;<sub>r</sub><sup>(&alpha;)</sup>&vert; &gt; &tau;<sub>f</sub><sup>(&alpha;)</sup> &or; &vert;&Delta;&mathbf{y}<sub>k-1</sub> + q&Delta;&Delta;&mathbf{y}<sub>k</sub>&vert; &gt; 0}</td>
    </tr>
    <tr>
      <td>18</td>
      <td>if q &gt; 0:</td>
    </tr>
    <tr>
      <td>19</td>
      <td>&Delta;&Delta;&mathbf{y}<sub>k</sub><sup>(&alpha;)</sup> &larr; q&Delta;&Delta;&mathbf{y}<sub>k</sub></td>
    </tr>
    <tr>
      <td>20</td>
      <td>&Delta;&mathbf{y}<sub>k</sub> &larr; &Delta;&mathbf{y}<sub>k-1</sub> + &Delta;&Delta;&mathbf{y}<sub>k</sub></td>
    </tr>
    <tr>
      <td>21</td>
      <td>Compute lattice rotation, update &Ccirc;<sub>e</sub> and Schmid tensors &hat{P}<sup>(&alpha;)</sup></td>
    </tr>
    <tr>
      <td>22</td>
      <td>Set search directions: d<sup>(&alpha;)</sup> &larr; normalised &alpha;-th column of &hat{P}<sup>T</sup>&Ccirc;<sub>e</sub> &hat{P}</td>
    </tr>
    <tr>
      <td>23</td>
      <td>&sigma;<sup>*</sup><sub>n+1</sub> &larr; &sigma;<sub>n</sub> + &Ccirc;<sub>e</sub> &Delta;D</td>
    </tr>
    <tr>
      <td>24</td>
      <td>r &larr; &radic;[&Sigma;<sub>&alpha; &isin; A<sub>k</sub></sub> (&vert; &hat{P}<sup>(&alpha;)</sup> &sigma;<sup>*</sup><sub>n+1</sub> - &hat{P}<sup>(&alpha;)</sup>&Ccirc;<sub>e</sub> &hat{P}&Delta;&mathbf{y}&vert; - &tau;<sub>f</sub><sup>(&alpha;)</sup>)<sup>2</sup>], A<sub>k</sub> = {&alpha; &isin; T: &vert;&tau;<sub>r</sub><sup>(&alpha;)</sup>&vert; &gt; &tau;<sub>f</sub><sup>(&alpha;)</sup> &or; &vert;&Delta;&mathbf{y}<sub>k</sub>&vert; &gt; 0}</td>
    </tr>
    <tr>
      <td colspan="2"><b>3. Stress update and algorithmic tangent computation:</b></td>
    </tr>
    <tr>
      <td>25</td>
      <td>&sigma;<sub>n+1</sub> = &sigma;<sub>n</sub> + &Ccirc;<sub>e</sub> &Delta;D - &Ccirc;<sub>e</sub> &hat{P}&Delta;&mathbf{y}</td>
    </tr>
    <tr>
      <td>26</td>
      <td>&Ccirc;<sub>ep</sub> = (&Ccirc;<sub>e</sub><sup>-1</sup> + &hat{P}&mathbf{G}^{-1}&hat{P}<sup>T</sup>)<sup>-1</sup></td>
    </tr>
  </tbody>
</table>
After convergence of the FPI scheme, Eqs. 8 and 9 hold for all active slip systems. As such, the stress increment can be computed using Eq. 8 as

$$ \Delta \underline{\sigma} = \underline{\mathcal{C}}_e \Delta \underline{\mathbf{D}} - \underline{\mathcal{C}}_e \hat{\mathbf{P}} \Delta \underline{\mathbf{\gamma}}, \qquad (16) $$

and the consistent algorithmic tangent $\underline{\mathcal{C}}_{ep}$ can be derived as [12]

$$ \underline{\mathcal{C}}_{ep} = (\underline{\mathcal{C}}_e^{-1} + \hat{\mathbf{P}}\mathbf{G}^{-1}\hat{\mathbf{P}}^{\mathrm{T}})^{-1}. \qquad (17) $$

Here, $G^{(\alpha\beta)} = \partial\phi^{(\alpha)}/\partial\gamma^{(\beta)}$.

<table><caption>Table 3: Fixed-point iteration based stress update algorithm for rate-independent CP</caption>
  <tbody>
    <tr>
      <td colspan="2"><b>1. Initialisation:</b></td>
    </tr>
    <tr>
      <td>1</td>
      <td>k &larr; 0</td>
    </tr>
    <tr>
      <td>2</td>
      <td>&Delta;&mathbf{y}<sub>k</sub> &larr; 0</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Compute lattice rotation, update &Ccirc;<sub>e</sub> and Schmid tensors &hat{P}<sup>(&alpha;)</sup></td>
    </tr>
    <tr>
      <td>4</td>
      <td>Set search directions: d<sup>(&alpha;)</sup> &larr; normalised &alpha;-th column of &hat{P}<sup>T</sup>&Ccirc;<sub>e</sub> &hat{P}</td>
    </tr>
    <tr>
      <td>5</td>
      <td>&sigma;<sup>*</sup><sub>n+1</sub> &larr; &sigma;<sub>n</sub> + &Ccirc;<sub>e</sub> &Delta;D</td>
    </tr>
    <tr>
      <td>6</td>
      <td>r<sub>0</sub> = min(1.0, &radic;[&Sigma;<sub>&alpha; &isin; A<sub>k</sub></sub> (&vert; &hat{P}<sup>(&alpha;)</sup><sup>T</sup> &sigma;<sup>*</sup><sub>n+1</sub> &vert; &minus; &tau;<sub>f</sub><sup>(&alpha;)</sup>)<sup>2</sup>), A<sub>k</sub> = {&alpha; &isin; T: &vert;&tau;<sub>r</sub><sup>(&alpha;)</sup>&vert; &gt; &tau;<sub>f</sub><sup>(&alpha;)</sup>}</td>
    </tr>
    <tr>
      <td>7</td>
      <td>r &larr; r<sub>0</sub></td>
    </tr>
    <tr>
      <td colspan="2"><b>2. FPI-iterations:</b></td>
    </tr>
    <tr>
      <td>8</td>
      <td>do while r &gt; 10<sup>-4</sup></td>
    </tr>
    <tr>
      <td>9</td>
      <td>k &larr; k + 1</td>
    </tr>
    <tr>
      <td>10</td>
      <td>for &alpha; = 1, 2, ..., n<sub>s</sub> do</td>
    </tr>
    <tr>
      <td>11</td>
      <td>find t<sup>(&alpha;)</sup>: solve &phi;&#773;<sup>(&alpha;)</sup>(&mathbf{y}<sub>n</sub>, &Delta;&mathbf{y}<sub>k&minus;1</sub> + t<sup>(&alpha;)</sup>d<sup>(&alpha;)</sup>) = 0</td>
    </tr>
    <tr>
      <td>12</td>
      <td>&Delta;&Delta;&mathbf{y}<sub>k</sub><sup>(&alpha;)</sup> &larr; t<sup>(&alpha;)</sup>d<sup>(&alpha;&alpha;)</sup></td>
    </tr>
    <tr>
      <td>13</td>
      <td>for &alpha; = 1, 2, ..., n<sub>s</sub> do</td>
    </tr>
    <tr>
      <td>14</td>
      <td>if &Delta;&mathbf{y}<sub>k&minus;1</sub><sup>(&alpha;)</sup> + &Delta;&Delta;&mathbf{y}<sub>k</sub><sup>(&alpha;)</sup> = 0:</td>
    </tr>
    <tr>
      <td>15</td>
      <td>&Delta;&mathbf{y}<sub>k&minus;1</sub><sup>&alpha;</sup> &larr; 0</td>
    </tr>
    <tr>
      <td>16</td>
      <td>&Delta;&Delta;&mathbf{y}<sub>k</sub><sup>(&alpha;)</sup> &larr; 0</td>
    </tr>
    <tr>
      <td>17</td>
      <td>find q: minimize &radic;[&Sigma;<sub>&alpha; &isin; A<sub>k</sub></sub> &phi;&#773;<sup>(&alpha;)</sup>(&mathbf{y}<sub>n</sub>, &Delta;&mathbf{y}<sub>k&minus;1</sub> + q&Delta;&Delta;&mathbf{y}<sub>k</sub>)<sup>2</sup>], A<sub>k</sub> = {&alpha; &isin; T: &vert;&tau;<sub>r</sub><sup>(&alpha;)</sup>&vert; &gt; &tau;<sub>f</sub><sup>(&alpha;)</sup> &or; &vert;&Delta;&mathbf{y}<sub>k&minus;1</sub> + q&Delta;&Delta;&mathbf{y}<sub>k</sub>&vert; &gt; 0}</td>
    </tr>
    <tr>
      <td>18</td>
      <td>if q &gt; 0:</td>
    </tr>
    <tr>
      <td>19</td>
      <td>&Delta;&Delta;&mathbf{y}<sub>k</sub><sup>(&alpha;)</sup> &larr; q&Delta;&Delta;&mathbf{y}<sub>k</sub></td>
    </tr>
    <tr>
      <td>20</td>
      <td>&Delta;&mathbf{y}<sub>k</sub> &larr; &Delta;&mathbf{y}<sub>k&minus;1</sub> + &Delta;&Delta;&mathbf{y}<sub>k</sub></td>
    </tr>
    <tr>
      <td>21</td>
      <td>Compute lattice rotation, update &Ccirc;<sub>e</sub> and Schmid tensors &hat{P}<sup>(&alpha;)</sup></td>
    </tr>
    <tr>
      <td>22</td>
      <td>Set search directions: d<sup>(&alpha;)</sup> &larr; normalised &alpha;-th column of &hat{P}<sup>T</sup>&Ccirc;<sub>e</sub> &hat{P}</td>
    </tr>
    <tr>
      <td>23</td>
      <td>&sigma;<sup>*</sup><sub>n+1</sub> &larr; &sigma;<sub>n</sub> + &Ccirc;<sub>e</sub> &Delta;D</td>
    </tr>
    <tr>
      <td>24</td>
      <td>r &larr; &radic;[&Sigma;<sub>&alpha; &isin; A<sub>k</sub></sub> (&vert; &hat{P}<sup>(&alpha;)</sup><sup>T</sup> &sigma;<sup>*</sup><sub>n+1</sub> &minus; &hat{P}<sup>(&alpha;)</sup><sup>T</sup>&Ccirc;<sub>e</sub> &hat{P}&Delta;&mathbf{y}&vert; &minus; &tau;<sub>f</sub><sup>(&alpha;)</sup>)<sup>2</sup>], A<sub>k</sub> = {&alpha; &isin; T: &vert;&tau;<sub>r</sub><sup>(&alpha;)</sup>&vert; &gt; &tau;<sub>f</sub><sup>(&alpha;)</sup> &or; &vert;&Delta;&mathbf{y}<sub>k</sub>&vert; &gt; 0}</td>
    </tr>
    <tr>
      <td colspan="2"><b>3. Stress update and algorithmic tangent computation:</b></td>
    </tr>
    <tr>
      <td>25</td>
      <td>&sigma;<sub>n+1</sub> = &sigma;<sub>n</sub> + &Ccirc;<sub>e</sub> &Delta;D - &Ccirc;<sub>e</sub> &hat{P}&Delta;&mathbf{y}</td>
    </tr>
    <tr>
      <td>26</td>
      <td>&Ccirc;<sub>ep</sub> = (&Ccirc;<sub>e</sub><sup>-1</sup> + &hat{P}&mathbf{G}^{-1}&hat{P}<sup>T</sup>)<sup>-1</sup></td>
    </tr>
  </tbody>
</table>
## Results and discussion

In this section, the stability of the novel FPI-based stress update algorithm will be assessed. Its behaviour is compared against that of a conventional Newton Raphson (NR)-based CPP algorithm with an iterative active set update procedure. The NR-based algorithm, implemented as described in [12], has two options to initialise the active set at the start of the increment. It could (i) be initialised as the active set of the previous increment or (ii) include all systems for which the trial stress renders $\tau_r^{(\alpha)} > \tau_f^{(\alpha)}$. Analyses with the NR-based algorithm are aborted if an update to the active set results in a set that has already been tried (and failed) during an increment. Both algorithms are implemented in a FORTRAN module that returns the increments of stress and plastic slip of a material point for a given deformation increment. All analyses below are conducted for an initially stress-free FCC single crystal of which the lattice coordinate system is aligned with the global coordinate system at the start of the analysis.

Stability of the solution scheme for a deformation-driven problem. We first consider a simple shear test in which a material point is deformed by prescribing a total deformation gradient tensor with $F_{12} = 4.0$ in 10, 100 and 1000 increments. Solutions are attempted with both the FPI-based and the NR-based algorithm. For the latter, a converged solution was only reached for the test with 1000 increments using option (ii) to initialise the active set. The FPI-based algorithm, on the other hand, converged for all increment counts considered.

Figure 2 shows the resulting stress-strain curves. The solutions obtained with the FPI scheme stably converge with an increasing number of increments. For the simulation with 1000 increments, the FPI-based algorithm yields the same result as the NR-based algorithm. The fact that the FPI-based algorithm correctly predicted the material response in as little as 10 increments shows that it converges to a solution even in the presence of large rotations. The FPI scheme thus allows for a large increment size.

Figure 2 furthermore shows that only very few fixed-point iterations per increment are required during the simulations with 100 and 1000 increments. At 1000 increments, the average number of fixed-point iterations is comparable to the average number of NR-iterations. Recall that each fixed-point iteration requires only function evaluations and simple bracketing operations to compute the solution to the 12 (for FCC) one-dimensional root-finding problems (line 11 of Table 3) and the scalar optimisation problem (line 17 of Table 3), whereas an iteration with the NR-based algorithm requires to solve a system of equations of which the size equals the number of active slip systems. For the simple shear problem solved in 1000 increments, the average number of iterations in solving each root-finding problem is only 11, that to solve the scalar optimisation problem is 23, while the number of active system in the majority of the increments equals 8. These numbers and the iteration counts in Figure 2 indicate that the novel FPI scheme can be competitive with a conventional NR-based stress update scheme, especially when the latter requires multiple restarts due to the iterative active set search. A more detailed analysis of the computational time is, however, required to further asses the performance of the FPI scheme.

Influence of lattice rotations. To investigate the effect of lattice rotations on the convergence behaviour of the FPI scheme, we reran the simple shear test with 10 and 100 increments while turning off the lattice rotations during the fixed-point iterations (Lines 21 – 23 in Table 3). Rotations are now updated in a staggered manner only after convergence of the FPI scheme. For the simulation with 100 increments, the resulting stress-strain curve in Figure 3 lags only slightly behind that of the corresponding analysis in which lattice rotations during the fixed-point iterations were enabled. For the simulation with 100 increments, lattice rotations during each increment are thus small. Figure 3 shows that in such a case, enabling lattice rotations during the fixed-point iterations does not significantly affect the number of fixed-point iterations.
## Results and discussion

In this section, the stability of the novel FPI-based stress update algorithm will be assessed. Its behaviour is compared against that of a conventional Newton Raphson (NR)-based CPP algorithm with an iterative active set update procedure. The NR-based algorithm, implemented as described in [12], has two options to initialise the active set at the start of the increment. It could (i) be initialised as the active set of the previous increment or (ii) include all systems for which the trial stress renders $\tau_r^{(\alpha)} > \tau_f^{(\alpha)}$. Analyses with the NR-based algorithm are aborted if an update to the active set results in a set that has already been tried (and failed) during an increment. Both algorithms are implemented in a FORTRAN module that returns the increments of stress and plastic slip of a material point for a given deformation increment. All analyses below are conducted for an initially stress-free FCC single crystal of which the lattice coordinate system is aligned with the global coordinate system at the start of the analysis.

Stability of the solution scheme for a deformation-driven problem. We first consider a simple shear test in which a material point is deformed by prescribing a total deformation gradient tensor with $F_{12} = 4.0$ in 10, 100 and 1000 increments. Solutions are attempted with both the FPI-based and the NR-based algorithm. For the latter, a converged solution was only reached for the test with 1000 increments using option (ii) to initialise the active set. The FPI-based algorithm, on the other hand, converged for all increment counts considered.

Figure 2 shows the resulting stress-strain curves. The solutions obtained with the FPI scheme stably converge with an increasing number of increments. For the simulation with 1000 increments, the FPI-based algorithm yields the same result as the NR-based algorithm. The fact that the FPI-based algorithm correctly predicted the material response in as little as 10 increments shows that it converges to a solution even in the presence of large rotations. The FPI scheme thus allows for a large increment size.

Figure 2 furthermore shows that only very few fixed-point iterations per increment are required during the simulations with 100 and 1000 increments. At 1000 increments, the average number of fixed-point iterations is comparable to the average number of NR-iterations. Recall that each fixed-point iteration requires only function evaluations and simple bracketing operations to compute the solution to the 12 (for FCC) one-dimensional root-finding problems (line 11 of Table 3) and the scalar optimisation problem (line 17 of Table 3), whereas an iteration with the NR-based algorithm requires to solve a system of equations of which the size equals the number of active slip systems. For the simple shear problem solved in 1000 increments, the average number of iterations in solving each root-finding problem is only 11, that to solve the scalar optimisation problem is 23, while the number of active system in the majority of the increments equals 8. These numbers and the iteration counts in Figure 2 indicate that the novel FPI scheme can be competitive with a conventional NR-based stress update scheme, especially when the latter requires multiple restarts due to the iterative active set search. A more detailed analysis of the computational time is, however, required to further asses the performance of the FPI scheme.

Influence of lattice rotations. To investigate the effect of lattice rotations on the convergence behaviour of the FPI scheme, we reran the simple shear test with 10 and 100 increments while turning off the lattice rotations during the fixed-point iterations (Lines 21 – 23 in Table 3). Rotations are now updated in a staggered manner only after convergence of the FPI scheme. For the simulation with 100 increments, the resulting stress-strain curve in Figure 3 lags only slightly behind that of the corresponding analysis in which lattice rotations during the fixed-point iterations were enabled. For the simulation with 100 increments, lattice rotations during each increment are thus small. Figure 3 shows that in such a case, enabling lattice rotations during the fixed-point iterations does not significantly affect the number of fixed-point iterations.
<figure><img src="image_2028.png" /><figcaption>Figure 2: Stress – strain curves for a single crystal under simple shear deformation, prescribed in 10, 100 and 1000 increments, as predicted by a conventional NR-based [12] and the proposed FPI-based stress update algorithm, along with the average number of NR or FPI iterations per increment.</figcaption></figure>

<figure><img src="image_2029.png" /><figcaption>Figure 3: Stress – strain curves for a single crystal under simple shear deformation, prescribed in 10 and 100 increments, as predicted by the proposed FPI-based stress update algorithm with (Mode = cont) and without (Mode = stag) rotation updates during the FPI-iterations, along with the average number of FPI iterations per increment.</figcaption></figure>

This is different for the simulation with 10 increments. The resulting stress-strain curve in Figure 3 differs significantly from that of the corresponding analysis with lattice rotations during the fixed-point iterations enabled. This indicates that the fixed-point now moves considerably during the iterations. Clearly, staggered rotation updates do not yield accurate results for such a limited number of increments. Nonetheless, for this test, enabling rotation updates during fixed-point iterations increases the average number of fixed-point iterations by a factor four. Severe lattice rotations within an increment could thus slow down the convergence of the FPI scheme.

Stability of the solution scheme for a stress-driven problem. The tests considered so far were deformation-driven. As such, no global equilibrium iterations were needed to obtain a converged solution. To assess the stability of the proposed FPI-based stress update algorithm within such equilibrium iterations, we now consider a stress-driven problem.

For this problem, the incremental deformation rate $\Delta\underline{D}_n$ that produces the prescribed stress $\underline{s}_{n+1}$ at the end of an increment $n$ has to be iteratively computed. At the end of the $k^{th}$ equilibrium iteration, $\Delta\underline{D}_n$ is updated as $\Delta\underline{D}_n^{k+1} = \Delta\underline{D}_n^k + \Delta\Delta\underline{D}_n^k$, with $\Delta\Delta\underline{D}_n^k = \underline{C}_{ep}^{-1}(\underline{s}_{n+1} - \underline{s}^k)$. The bracketed term in the last expression is the difference between the prescribed stress state and the stress as currently predicted by the stress update algorithm. Equilibrium iterations continue while the magnitude of the largest element in $(\underline{s}_{n+1} - \underline{s}^k)$ is larger than $10^{-3}$ times the magnitude of the largest element in the initial difference vector $(\underline{s}_{n+1} - \underline{s}^1)$.

We considered a total of 968 different load cases in which the normal components of the stress tensor are prescribed. These components are given by 968 points that are evenly distributed on a sphere with a radius of 400 MPa to cover a large space of possible load directions. The Fibonacci lattice [17] is used to distribute the points. Each load case is applied in 10 increments and solved both with the FPI-based and the NR-based stress update algorithms. The FPI-based algorithm converged in 96.3% of the cases, whereas the NR-based algorithm converged in only 76.4% of the cases. For this stress-drive problem, the NR-based algorithm was found to converge best when using option (i) to initialise the active set.
<figure><img src="image_2020.png" /><figcaption>Figure 2: Stress – strain curves for a single crystal under simple shear deformation, prescribed in 10, 100 and 1000 increments, as predicted by a conventional NR-based [12] and the proposed FPI-based stress update algorithm, along with the average number of NR or FPI iterations per increment.</figcaption></figure>

<figure><img src="image_2021.png" /><figcaption>Figure 3: Stress – strain curves for a single crystal under simple shear deformation, prescribed in 10 and 100 increments, as predicted by the proposed FPI-based stress update algorithm with (Mode = cont) and without (Mode = stag) rotation updates during the FPI-iterations, along with the average number of FPI iterations per increment.</figcaption></figure>

This is different for the simulation with 10 increments. The resulting stress-strain curve in Figure 3 differs significantly from that of the corresponding analysis with lattice rotations during the fixed-point iterations enabled. This indicates that the fixed-point now moves considerably during the iterations. Clearly, staggered rotation updates do not yield accurate results for such a limited number of increments. Nonetheless, for this test, enabling rotation updates during fixed-point iterations increases the average number of fixed-point iterations by a factor four. Severe lattice rotations within an increment could thus slow down the convergence of the FPI scheme.

Stability of the solution scheme for a stress-driven problem. The tests considered so far were deformation-driven. As such, no global equilibrium iterations were needed to obtain a converged solution. To assess the stability of the proposed FPI-based stress update algorithm within such equilibrium iterations, we now consider a stress-driven problem.

For this problem, the incremental deformation rate $\Delta\underline{D}_n$ that produces the prescribed stress $\underline{s}_{n+1}$ at the end of an increment $n$ has to be iteratively computed. At the end of the $k^{th}$ equilibrium iteration, $\Delta\underline{D}_n$ is updated as $\Delta\underline{D}_n^{k+1} = \Delta\underline{D}_n^k + \Delta\Delta\underline{D}_n^k$, with $\Delta\Delta\underline{D}_n^k = \underline{C}_{ep}^{-1}(\underline{s}_{n+1} - \underline{s}^k)$. The bracketed term in the last expression is the difference between the prescribed stress state and the stress as currently predicted by the stress update algorithm. Equilibrium iterations continue while the magnitude of the largest element in $(\underline{s}_{n+1} - \underline{s}^k)$ is larger than $10^{-3}$ times the magnitude of the largest element in the initial difference vector $(\underline{s}_{n+1} - \underline{s}^1)$.

We considered a total of 968 different load cases in which the normal components of the stress tensor are prescribed. These components are given by 968 points that are evenly distributed on a sphere with a radius of 400 MPa to cover a large space of possible load directions. The Fibonacci lattice [17] is used to distribute the points. Each load case is applied in 10 increments and solved both with the FPI-based and the NR-based stress update algorithms. The FPI-based algorithm converged in 96.3% of the cases, whereas the NR-based algorithm converged in only 76.4% of the cases. For this stress-drive problem, the NR-based algorithm was found to converge best when using option (i) to initialise the active set.
The fact that not all load cases have converged is caused by a well-known property of rate-independent CPP algorithms. If the active set changes during equilibrium iterations, the resulting changes in the algorithmic tangent can lead to non-convergence of these iterations [18]. The FPI scheme still uses an active set and hence is not invulnerable to this behaviour. Nonetheless, it does not suffer from breakdowns of the iterative active set search procedure that cause the NR-based algorithm to abort prematurely, and therefore allows more of the load cases to converge. Given that the FPI scheme is able to resolve the vast majority of load cases, it seems a promising alternative to CPP algorithms with iterative active set search.

## Conclusion

A novel stress update algorithm for rate-independent crystal plasticity was presented in this work. It was developed based on the observation that the set of active slip systems does not have to be explicitly initialised and updated if, for each of the crystal's $n_s$ slip systems, a single condition is formulated that has to be zero both when the system is active as well as when it is inactive. The slip state for which each of these conditions is zero is defined by a hypersurface in $n_s$-dimensional space. We proposed a fixed-point iteration (FPI) scheme to trace these hypersurfaces until they all intersect. The resulting stress update algorithm simultaneously provides the set of active slip systems and their slip activity, without inducing spurious stresses on nonactive systems. It stably converges to a solution even when the crystal lattice rotates significantly during iterations, allowing for large increment sizes. The proposed stress update algorithm is furthermore shown to be efficient and significantly more stable than a conventional Newton-Raphson-based stress update scheme featuring an iterative active set search procedure.

## References

[1] H. K. Akpama, M. B. Bettaieb, F. Abed-Meraim, Numerical integration of rate-independent BCC single crystal plasticity models: comparative study of two classes of numerical algorithms, Int J Numerer Meth Eng, 108:5 (2016) 363-422. https://doi.org/10.1002/nme.5215

[2] C. Miehe, J. Schröder, A comparative study of stress update algorithms for rate-independent and rate-dependent crystal plasticity, Int J Numer Meth Eng, 50:2 (2001) 273-298. https://doi.org/10.1002/1097-0207(20010120)50:2<273::AID-NME17>3.3.CO;2-H

[3] C. Miehe, J. Schotte, Crystal Plasticity and Evolution of Polycrystalline Microstructure, in: E. Stein, R. De Borst, and T. J. R. Hughes (Eds.), Encyclopedia of Computational Mechanics Second Edition, John Wiley & Sons, Ltd. 2004, pp. 267-289.

[4] R. McGinty, D. McDowell, A semi-implicit integration scheme for rate independent finite crystal plasticity, Int J Plasticity, 22:6 (2006) 996-1025. https://doi.org/10.1016/j.ijplas.2005.06.002

[5] C. Miehe, D. Rosato, Fast texture updates in fcc polycrystal plasticity based on a linear active-set-estimate of the lattice spin, J Mech Phys Solids, 55:12 (2007) 2687-2716. https://doi.org/10.1016/j.jmps.2007.04.002

[6] M. Ben Bettaieb, O. Débordes, A. Dogui, L. Duchêne, C. Keller, On the numerical integration of rate independent single crystal behavior at large strain, Int J Plasticity, 32-33 (2012) 184-217. https://doi.org/10.1016/j.ijplas.2011.10.010

[7] T. Mánik, H. M. Asadkandi, B. Holmedal, A robust algorithm for rate-independent crystal plasticity, Comput Meth Appl Mech Eng, 393 (2022) 114831. https://doi.org/10.1016/j.cma.2022.114831
The fact that not all load cases have converged is caused by a well-known property of rate-independent CPP algorithms. If the active set changes during equilibrium iterations, the resulting changes in the algorithmic tangent can lead to non-convergence of these iterations [18]. The FPI scheme still uses an active set and hence is not invulnerable to this behaviour. Nonetheless, it does not suffer from breakdowns of the iterative active set search procedure that cause the NR-based algorithm to abort prematurely, and therefore allows more of the load cases to converge. Given that the FPI scheme is able to resolve the vast majority of load cases, it seems a promising alternative to CPP algorithms with iterative active set search.

## Conclusion

A novel stress update algorithm for rate-independent crystal plasticity was presented in this work. It was developed based on the observation that the set of active slip systems does not have to be explicitly initialised and updated if, for each of the crystal's $n_s$ slip systems, a single condition is formulated that has to be zero both when the system is active as well as when it is inactive. The slip state for which each of these conditions is zero is defined by a hypersurface in $n_s$-dimensional space. We proposed a fixed-point iteration (FPI) scheme to trace these hypersurfaces until they all intersect. The resulting stress update algorithm simultaneously provides the set of active slip systems and their slip activity, without inducing spurious stresses on nonactive systems. It stably converges to a solution even when the crystal lattice rotates significantly during iterations, allowing for large increment sizes. The proposed stress update algorithm is furthermore shown to be efficient and significantly more stable than a conventional Newton-Raphson-based stress update scheme featuring an iterative active set search procedure.

## References

[1] H. K. Akpama, M. B. Bettaieb, F. Abed-Meraim, Numerical integration of rate-independent BCC single crystal plasticity models: comparative study of two classes of numerical algorithms, Int J Numerer Meth Eng, 108:5 (2016) 363-422. https://doi.org/10.1002/nme.5215

[2] C. Miehe, J. Schröder, A comparative study of stress update algorithms for rate-independent and rate-dependent crystal plasticity, Int J Numer Meth Eng, 50:2 (2001) 273-298. https://doi.org/10.1002/1097-0207(20010120)50:2<273::AID-NME17>3.3.CO;2-H

[3] C. Miehe, J. Schotte, Crystal Plasticity and Evolution of Polycrystalline Microstructure, in: E. Stein, R. De Borst, and T. J. R. Hughes (Eds.), Encyclopedia of Computational Mechanics Second Edition, John Wiley & Sons, Ltd. 2004, pp. 267-289.

[4] R. McGinty, D. McDowell, A semi-implicit integration scheme for rate independent finite crystal plasticity, Int J Plasticity, 22:6 (2006) 996-1025. https://doi.org/10.1016/j.ijplas.2005.06.002

[5] C. Miehe, D. Rosato, Fast texture updates in fcc polycrystal plasticity based on a linear active-set-estimate of the lattice spin, J Mech Phys Solids, 55:12 (2007) 2687-2716. https://doi.org/10.1016/j.jmps.2007.04.002

[6] M. Ben Bettaieb, O. Débordes, A. Dogui, L. Duchêne, C. Keller, On the numerical integration of rate independent single crystal behavior at large strain, Int J Plasticity, 32-33 (2012) 184-217. https://doi.org/10.1016/j.ijplas.2011.10.010

[7] T. Mánik, H. M. Asadkandi, B. Holmedal, A robust algorithm for rate-independent crystal plasticity, Comput Meth Appl Mech Eng, 393 (2022) 114831. https://doi.org/10.1016/j.cma.2022.114831
[8] S. Forest, M. B. Rubin, A rate-independent crystal plasticity model with a smooth elastic-plastic transition and no slip indeterminacy, Eur J Mech A-Solid, 55 (2016) 278-288.
https://doi.org/10.1016/j.euromechsol.2015.08.012

[9] F. Roters, P. Eisenlohr, T. R. Bieler, D. Raabe, Crystal Plasticity Finite Element Methods in Materials Science and Engineering. Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, 2010.
https://doi.org/10.1002/9783527631483

[10] M. Schmidt-Baldassari, Numerical concepts for rate-independent single crystal plasticity, Comput Meth Appl Mech Eng, 192:11-12 (2003) 1261-1280. https://doi.org/10.1016/S0045-7825(02)00563-7

[11] L. Scheunemann, P. S. B. Nigro, J. Schröder, P. M. Pimenta, A novel algorithm for rate independent small strain crystal plasticity based on the infeasible primal-dual interior point method, Int J Plasticity, 124 (2020) 1-19. https://doi.org/10.1016/j.ijplas.2019.07.020

[12] E. S. Perdahcıoğlu, A rate-independent crystal plasticity algorithm based on the interior point method, Comput Meth Appl Mech Eng, 418 (2024) 116533.
https://doi.org/10.1016/j.cma.2023.116533

[13] L. P. Kubin, B. Devincre, M. Tang, Mesoscopic modelling and simulation of plasticity in fcc and bcc crystals: Dislocation intersections and mobility, J Comput-Aided Mater Des, 5:1 (1998) 31-54. https://doi.org/10.1023/A:1008648120261

[14] M. Becker, "Incompatibility and Instability Based Size Effects in Crystals and Composites at Finite Elastoplastic Strains," PhD-thesis, Universität Stuttgart, Stuttgart, 2006.

[15] J. Mandel, Generalisation de la theorie de plasticite de W. T. Koiter, Int J Solids Struct, 1:3 (1965) 273-295. https://doi.org/10.1016/0020-7683(65)90034-X

[16] R. P. Brent, Algorithms for Minimization without Derivatives. Prentice-Hall, New Jersey, 1973.

[17] Á. González, Measurement of Areas on a Sphere Using Fibonacci and Latitude-Longitude Lattices, Math Geosci, 42:1 (2009) 49-64. https://doi.org/10.1007/s11004-009-9257-x

[18] T. Mánik, B. Holmedal, Review of the Taylor ambiguity and the relationship between rate-independent and rate-dependent full-constraints Taylor models, Int J Plasticity, 55 (2014) 152-181. https://doi.org/10.1016/j.ijplas.2013.10.002
[8] S. Forest, M. B. Rubin, A rate-independent crystal plasticity model with a smooth elastic-plastic transition and no slip indeterminacy, Eur J Mech A-Solid, 55 (2016) 278-288.
https://doi.org/10.1016/j.euromechsol.2015.08.012

[9] F. Roters, P. Eisenlohr, T. R. Bieler, D. Raabe, Crystal Plasticity Finite Element Methods in Materials Science and Engineering. Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, 2010.
https://doi.org/10.1002/9783527631483

[10] M. Schmidt-Baldassari, Numerical concepts for rate-independent single crystal plasticity, Comput Meth Appl Mech Eng, 192:11-12 (2003) 1261-1280. https://doi.org/10.1016/S0045-7825(02)00563-7

[11] L. Scheunemann, P. S. B. Nigro, J. Schröder, P. M. Pimenta, A novel algorithm for rate independent small strain crystal plasticity based on the infeasible primal-dual interior point method, Int J Plasticity, 124 (2020) 1-19. https://doi.org/10.1016/j.ijplas.2019.07.020

[12] E. S. Perdahcıoğlu, A rate-independent crystal plasticity algorithm based on the interior point method, Comput Meth Appl Mech Eng, 418 (2024) 116533.
https://doi.org/10.1016/j.cma.2023.116533

[13] L. P. Kubin, B. Devincre, M. Tang, Mesoscopic modelling and simulation of plasticity in fcc and bcc crystals: Dislocation intersections and mobility, J Comput-Aided Mater Des, 5:1 (1998) 31-54. https://doi.org/10.1023/A:1008648120261

[14] M. Becker, "Incompatibility and Instability Based Size Effects in Crystals and Composites at Finite Elastoplastic Strains," PhD-thesis, Universität Stuttgart, Stuttgart, 2006.

[15] J. Mandel, Generalisation de la theorie de plasticite de W. T. Koiter, Int J Solids Struct, 1:3 (1965) 273-295. https://doi.org/10.1016/0020-7683(65)90034-X

[16] R. P. Brent, Algorithms for Minimization without Derivatives. Prentice-Hall, New Jersey, 1973.

[17] Á. González, Measurement of Areas on a Sphere Using Fibonacci and Latitude-Longitude Lattices, Math Geosci, 42:1 (2009) 49-64. https://doi.org/10.1007/s11004-009-9257-x

[18] T. Mánik, B. Holmedal, Review of the Taylor ambiguity and the relationship between rate-independent and rate-dependent full-constraints Taylor models, Int J Plasticity, 55 (2014) 152-181. https://doi.org/10.1016/j.ijplas.2013.10.002
# Fatigue life analysis of hot forming dies produced by the L-PBF and WA-DED additive technologies

ALIMOV Artem¹,ᵃ*, SVIRIDOV Alexander¹,ᵇ and HÄRTEL Sebastian¹,ᶜ

¹Chair of Hybrid Manufacturing, Brandenburg University of Technology, Cottbus, Germany

ᵃalimov@b-tu.de, ᵇsviridov@b-tu.de, ᶜhaertel@b-tu.de

**Keywords:** Fatigue, LCF, Additive Manufacturing, Hot Forming Dies, L-PBF, WA-DED, WAAM

**Abstract.** Additive manufacturing (AM) technologies are widely used for the production of individual complex precise parts and small-scale production. AM has been used in many applications for the manufacturing of sheet metal forming dies. The production of hot bulk forming dies by additive manufacturing technology has not been extensively studied in the literature. Traditional hot working steels exhibit many problems during the AM process, such as crack formation due to complex thermal processing cycles. Maraging/precipitation-hardened (PH) steels may be a good alternative because of their ability to undergo additive manufacturing. One of the main causes of hot forming die failures is thermomechanical fatigue. The fatigue performance of 17-4PH was investigated primarily during high-cycle fatigue (HCF) and low-cycle fatigue (LCF) modes without heat treatment as well as with H1100 heat treatment. The LCF of dies with the heat treatment designation H900, which is most appropriate for die production, has not been intensively investigated. In this work, the mechanical properties and fatigue life of 17-4PH precipitation-hardened steel produced by laser powder bed fusion (L-PBF) and wire arc direct energy deposition (WA-DED) additive technologies were experimentally investigated, and numerical models for high-temperature low cycle fatigue life were developed for 17-4PH steel. Based on the FEM simulation in the commercial program QForm UK, the fatigue life of the produced dies was predicted.

**Introduction**

Additive manufacturing (AM) technologies represent methods of layer-by-layer production of three-dimensional objects from a CAD model [1]. The main application fields of AM include the aerospace [2], energy [3] and automotive industries [4] as well as biomedicine [5]. Previously reported AMs were used primarily for rapid prototyping, for the production of individual parts or for small batch production [6], but the rapid development of AMs in recent years has extended the application of these materials to the production of casting molds [7] and sheet metal forming dies [8]. The additive manufacturing of hot bulk forging dies has not been covered as extensively in the literature. Huskic et al. [9] and Junker et al. [10] reported the additive manufacturing of forming dies from 1.2709 and 1.2343 tool steels, respectively. Typically, tool steels tend to develop hot cracks in additive manufacturing [11] and require special measures such as preheating the build plate to high temperatures of up to 800°C for defect-free production.

Currently, the most commonly used additive manufacturing methods are powder bed fusion (PBF) and direct energy deposition (DED) methods [12]. The energy source can be an electric arc, laser, plasma or electron beam. Laser Powder Bed Fusion (L-PBF) technology fabricates complex 3D components from powder by layer-by-layer selective laser melting. It offers one of the highest accuracies among AM technologies [13]. The main disadvantage of L-PBF is the limited workspace and relatively low performance. Wire-arc DED (WA-DED) has the highest productivity among the AM methods. Due to the use of an electric arc, WA-DED is believed to

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