specimen significantly affected the measured stress. More accurate processing of specimens is required in the future.

<figure><img src="image_1022.png" /><figcaption>Fig. 5 Experimental results using proposed apparatus, (a) validation of stress ratio, (b) comparison with yield loci modeled using standard biaxial tensile testing machine.</figcaption></figure>

## Conclusion

We developed a novel biaxial tensile testing apparatus that uses the transmission of a tensile load by a looped rope and conducted a biaxial tensile test on a mild steel sheet for preliminary testing. Although certain accuracy problems were observed, it was confirmed that the biaxial test could be performed at the intended load ratio. In the future, improvements will be made for its practical use as an efficient and reasonably priced apparatus for biaxial tensile testing.

## References

[1] ISO 16842: 2021, Metallic materials - Sheet and strip - Biaxial tensile testing method using a cruciform test piece.

[2] A. Makinde, L. Thibodeau, K. W. Neale, Development of an apparatus for biaxial testing for cruciform specimens, Exp. Mech. 32 (1992) 138-144. https://doi.org/10.1007/BF02324725

[3] J.P. Boehler, S. Demmerle, S. Koss, A new direct biaxial testing machine for anisotropic materials, Exp. Mech. 34 (1994) 1-9. https://doi.org/10.1007/BF02328435

[4] T. Kuwabara, S. Ikeda, T. Kuroda, Measurement and analysis of differential work hardening in cold-rolled steel sheet under biaxial tension, J. Material Process. Technol. 80/81(1998) 517-523. https://doi.org/10.1016/S0924-0136(98)00155-1

[5] G. Ferron, A. Makinde, Design and development of a biaxial strength testing device, J. Test. Evaluat. 16 (1988) 253-256. https://doi.org/10.1520/JTE10375J
specimen significantly affected the measured stress. More accurate processing of specimens is required in the future.

<figure><img src="image_1027.png" /><figcaption>Fig. 5 Experimental results using proposed apparatus, (a) validation of stress ratio, (b) comparison with yield loci modeled using standard biaxial tensile testing machine.</figcaption></figure>

## Conclusion

We developed a novel biaxial tensile testing apparatus that uses the transmission of a tensile load by a looped rope and conducted a biaxial tensile test on a mild steel sheet for preliminary testing. Although certain accuracy problems were observed, it was confirmed that the biaxial test could be performed at the intended load ratio. In the future, improvements will be made for its practical use as an efficient and reasonably priced apparatus for biaxial tensile testing.

## References

[1] ISO 16842: 2021, Metallic materials - Sheet and strip - Biaxial tensile testing method using a cruciform test piece.

[2] A. Makinde, L. Thibodeau, K. W. Neale, Development of an apparatus for biaxial testing for cruciform specimens, Exp. Mech. 32 (1992) 138-144. https://doi.org/10.1007/BF02324725

[3] J.P. Boehler, S. Demmerle, S. Koss, A new direct biaxial testing machine for anisotropic materials, Exp. Mech. 34 (1994) 1-9. https://doi.org/10.1007/BF02328435

[4] T. Kuwabara, S. Ikeda, T. Kuroda, Measurement and analysis of differential work hardening in cold-rolled steel sheet under biaxial tension, J. Material Process. Technol. 80/81(1998) 517-523. https://doi.org/10.1016/S0924-0136(98)00155-1

[5] G. Ferron, A. Makinde, Design and development of a biaxial strength testing device, J. Test. Evaluat. 16 (1988) 253-256. https://doi.org/10.1520/JTE10375J
[6] T. Nagayasu, S. Takahashi, T. Kuwabara, Development of compact biaxial tensile testing apparatus using conventional compression testing machine and evaluation of the test results, Proc. IDDRG 2010 May 31-June 2, 2010, Graz, Austria, 593-602.

[7] M. Merklein, M. Biasutti, Development of a biaxial tensile machine for characterization of sheet metals, J. Mater. Processing Technol. 213-6 (2013) 939-946. https://doi.org/10.1016/j.jmatprotec.2012.12.005

[8] T. Kuwabara, H. Nakano, 4th meeting meterial, Research committee on Advancement of steel sheet forming by advanced multiaxial stress tests, 2017, The iron and steel institute of Japan.

[9] Y. Hanabusa, H. Takizawa, T. Kuwabara, Numerical verification of a biaxial tensile test method using a cruciform specimen, J. Mater. Processing Technol. 213-6 (2013) 961-970. https://doi.org/10.1016/j.jmatprotec.2012.12.007
[6] T. Nagayasu, S. Takahashi, T. Kuwabara, Development of compact biaxial tensile testing apparatus using conventional compression testing machine and evaluation of the test results, Proc. IDDRG 2010 May 31-June 2, 2010, Graz, Austria, 593-602.

[7] M. Merklein, M. Biasutti, Development of a biaxial tensile machine for characterization of sheet metals, J. Mater. Processing Technol. 213-6 (2013) 939-946. https://doi.org/10.1016/j.jmatprotec.2012.12.005

[8] T. Kuwabara, H. Nakano, 4th meeting meterial, Research committee on Advancement of steel sheet forming by advanced multiaxial stress tests, 2017, The iron and steel institute of Japan.

[9] Y. Hanabusa, H. Takizawa, T. Kuwabara, Numerical verification of a biaxial tensile test method using a cruciform specimen, J. Mater. Processing Technol. 213-6 (2013) 961-970. https://doi.org/10.1016/j.jmatprotec.2012.12.007
# Enhancement of fracture strain during abrupt orthogonal strain-path changes in ferrite/martensite dual phase steel

MATSUNO Takashi<fnref n="1" />,ᵃ*, KINOSHITA Nanami<fnref n="1" />,ᵇ, HOKIMOTO Keisuke<fnref n="1" />,ᶜ, HAMA Takayuki<fnref n="2" />,ᵈ and HONDA Yoshiaki<fnref n="3" />,ᵉ

**Keywords:** Ductile Fracture, Strain-Path Change, Post-Necking, Dual Phase Steel

**Abstract** In multistage press forming, the scrutiny of work-hardening behavior attributed to abrupt strain path changes in deformation paths has been pivotal in the context of press formability. However, applying this method to stretch flanging presents a formidable complexity. In such scenarios, substantial deformations resulting from shear-cut processing and subsequent stretch-flange formation significantly influence the local fracture strains, a critical determinant of process success. Therefore, our investigation examined the impact of deformation path variations on the fracture strain of dual-phase steel sheets. A crucial facet of our experimental approach involved manipulating the strain level during the initial tensile deformation. Tensile deformation was extended beyond the onset of necking, culminating in the emergence of microvoids within the material. Following this preliminary extensive strain, successive tensile tests were conducted on diminutive round-bar specimens, subject to a 90° change in the deformation path, to investigate their influence on work hardening and fracture strain. The outcomes of the diameter-measuring tensile tests unveiled a pronounced cross effect during the early stages of deformation, particularly concerning work-hardening behavior. However, as the deformation advanced, the results conformed to the stress-strain curve observed in the same direction of tension. Notably, the fracture strain exhibited an approximately 10% enhancement attributable to the abrupt strain path change, thereby having the potential to augment the practical stretch flangeability.

**Introduction**

In press forming, the deformation path of a material commonly changes partway through the process. Processes involving multiple dies, such as progressive presses, often entail abrupt alterations in the deformation paths. Forming processes like stretch flanging may incorporate a 90° abrupt strain path change, particularly when the initial deformation is characterized as shear-cut processing. The ductility of automotive crash-resistant components, which are manufactured via press forming, depends on how they respond to an abrupt strain path change between press forming and crushing deformations.

Dortance, despite its impw reports address the fracture strain associated with sudden changes in the deformation path (local plastic strain at fracture). In single-phase steels, an abrupt strain path change only alters work hardening, with minimal impact on elongation [1]. However, elongation increases due to an abrupt strain path change in multi-phase steels like ferrite-pearlite steels [1].

Therefore, this study investigates work hardening and fracture strain in response to pre-strain surpassing plastic instability in dual-phase (DP) steels. Tiny round bar tensile specimens were cut

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

<footnote n="1">¹Faculty of Engineering, Tottori University, 4-101 Koyama-cho-minami, Tottori, Tottori 680-8552, Japan</footnote>

<footnote n="2">²Graduate School of Energy Science, Kyoto University, Yoshida-Honmachi, Sakyo-ku, Kyoto 606-8501, Japan</footnote>

<footnote n="3">³Nippon Steel Corporation, 20-1 Shintomi, Futtsu-City, Chiba, 293-8511 Japan</footnote>

<footnote>ᵃmatsu@tottori-u.ac.jp, ᵇm23J3015z@edu.tottori-u.ac.jp, ᶜb20T1092z@edu.tottori-u.ac.jp,</footnote>

<footnote>ᵈhama@energy.kyoto-u.ac.jp, ᶜhonda.x66.yoshiaki@jp.nipponsteel.com</footnote>
# Enhancement of fracture strain during abrupt orthogonal strain-path changes in ferrite/martensite dual phase steel

MATSUNO Takashi<fnref n="1" />,ᵃ*, KINOSHITA Nanami<fnref n="1" />,ᵇ, HOKIMOTO Keisuke<fnref n="1" />,ᶜ, HAMA Takayuki<fnref n="2" />,ᵈ and HONDA Yoshiaki<fnref n="3" />,ᵉ

**Keywords:** Ductile Fracture, Strain-Path Change, Post-Necking, Dual Phase Steel

**Abstract** In multistage press forming, the scrutiny of work-hardening behavior attributed to abrupt strain path changes in deformation paths has been pivotal in the context of press formability. However, applying this method to stretch flanging presents a formidable complexity. In such scenarios, substantial deformations resulting from shear-cut processing and subsequent stretch-flange formation significantly influence the local fracture strains, a critical determinant of process success. Therefore, our investigation examined the impact of deformation path variations on the fracture strain of dual-phase steel sheets. A crucial facet of our experimental approach involved manipulating the strain level during the initial tensile deformation. Tensile deformation was extended beyond the onset of necking, culminating in the emergence of microvoids within the material. Following this preliminary extensive strain, successive tensile tests were conducted on diminutive round-bar specimens, subject to a 90° change in the deformation path, to investigate their influence on work hardening and fracture strain. The outcomes of the diameter-measuring tensile tests unveiled a pronounced cross effect during the early stages of deformation, particularly concerning work-hardening behavior. However, as the deformation advanced, the results conformed to the stress-strain curve observed in the same direction of tension. Notably, the fracture strain exhibited an approximately 10% enhancement attributable to the abrupt strain path change, thereby having the potential to augment the practical stretch flangeability.

**Introduction**

In press forming, the deformation path of a material commonly changes partway through the process. Processes involving multiple dies, such as progressive presses, often entail abrupt alterations in the deformation paths. Forming processes like stretch flanging may incorporate a 90° abrupt strain path change, particularly when the initial deformation is characterized as shear-cut processing. The ductility of automotive crash-resistant components, which are manufactured via press forming, depends on how they respond to an abrupt strain path change between press forming and crushing deformations.

Dortance, despite its impw reports address the fracture strain associated with sudden changes in the deformation path (local plastic strain at fracture). In single-phase steels, an abrupt strain path change only alters work hardening, with minimal impact on elongation [1]. However, elongation increases due to an abrupt strain path change in multi-phase steels like ferrite-pearlite steels [1].

Therefore, this study investigates work hardening and fracture strain in response to pre-strain surpassing plastic instability in dual-phase (DP) steels. Tiny round bar tensile specimens were cut

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

<footnote n="1">¹Faculty of Engineering, Tottori University, 4-101 Koyama-cho-minami, Tottori, Tottori 680-8552, Japan</footnote>

<footnote n="2">²Graduate School of Energy Science, Kyoto University, Yoshida-Honmachi, Sakyo-ku, Kyoto 606-8501, Japan</footnote>

<footnote n="3">³Nippon Steel Corporation, 20-1 Shintomi, Futtsu-City, Chiba, 293-8511 Japan</footnote>

<footnote>ᵃmatsu@tottori-u.ac.jp, ᵇm23J3015z@edu.tottori-u.ac.jp, ᶜb20T1092z@edu.tottori-u.ac.jp,</footnote>

<footnote>ᵈhama@energy.kyoto-u.ac.jp, ᶜhonda.x66.yoshiaki@jp.nipponsteel.com</footnote>
from sheet-type tensile tested specimens that had undergone pre-straining exceeding the onset of necking. Subsequent diameter measurement-type tensile tests were executed [2-4] in the same direction and at a 90° angle to the pre-strain.

## Experimental Methods

### Material

Two DP steels, DP-A and DP-B, with tensile strengths of 560 and 650 MPa, respectively, were used for the tests. The chemical compositions and mechanical properties are detailed in Tables 1 and 2, respectively. These materials exhibited hard martensite islands dispersed within a soft ferrite matrix. Figure 1 shows the micrographs of these microstructures. Specifically, DP-A and DP-B comprised 4.3% and 10.5% martensite volume fraction, respectively.

<table><caption>Table 1 Chemical components of tested steels (mass%)</caption><thead><tr><th></th><th>C</th><th>Si</th><th>Mn</th></tr></thead><tbody><tr><th>DP-A</th><td>0.038</td><td>0.50</td><td>1.50</td></tr><tr><th>DP-B</th><td>0.067</td><td>0.50</td><td>1.50</td></tr></tbody></table>

<table><caption>Table 2 Mechanical properties of tested steels</caption><thead><tr><th></th><th>Yield stress<br>[MPa]</th><th>Tensile strength<br>[MPa]</th><th>Uniform<br>elongation [%]</th><th>Total elongation<br>[%]</th></tr></thead><tbody><tr><th>DP-A</th><td>335</td><td>561</td><td>19</td><td>31</td></tr><tr><th>DP-B</th><td>353</td><td>650</td><td>16</td><td>24</td></tr></tbody></table>

<figure><img src="image_1024.png" /><figcaption>Fig. 1 Micrographs of the DP steels used in this study: (a) DP-A containing 4.3% martensite and (b) DP-B containing 10.5% martensite. Grains colored in white are martensite, and remaining brown grains are ferrite.</figcaption></figure>

### Preparing tiny round-bar specimen

We subjected No. 5 type specimens, 2.5-mm thick and standardized by Japanese Industrial Standards (JIS) Z2241, to preliminary tensile testing before subsequent compact tensile tests. These specimens were elongated by up to 23.14% and 19.44% for DP-A and DP-B, respectively,
from sheet-type tensile tested specimens that had undergone pre-straining exceeding the onset of necking. Subsequent diameter measurement-type tensile tests were executed [2-4] in the same direction and at a 90° angle to the pre-strain.

## Experimental Methods

### Material

Two DP steels, DP-A and DP-B, with tensile strengths of 560 and 650 MPa, respectively, were used for the tests. The chemical compositions and mechanical properties are detailed in Tables 1 and 2, respectively. These materials exhibited hard martensite islands dispersed within a soft ferrite matrix. Figure 1 shows the micrographs of these microstructures. Specifically, DP-A and DP-B comprised 4.3% and 10.5% martensite volume fraction, respectively.

<table><caption>Table 1 Chemical components of tested steels (mass%)</caption><thead><tr><th></th><th>C</th><th>Si</th><th>Mn</th></tr></thead><tbody><tr><th>DP-A</th><td>0.038</td><td>0.50</td><td>1.50</td></tr><tr><th>DP-B</th><td>0.067</td><td>0.50</td><td>1.50</td></tr></tbody></table>

<table><caption>Table 2 Mechanical properties of tested steels</caption><thead><tr><th></th><th>Yield stress<br>[MPa]</th><th>Tensile strength<br>[MPa]</th><th>Uniform<br>elongation [%]</th><th>Total elongation<br>[%]</th></tr></thead><tbody><tr><th>DP-A</th><td>335</td><td>561</td><td>19</td><td>31</td></tr><tr><th>DP-B</th><td>353</td><td>650</td><td>16</td><td>24</td></tr></tbody></table>

<figure><img src="image_1029.png" /><figcaption>Fig. 1 Micrographs of the DP steels used in this study: (a) DP-A containing 4.3% martensite and (b) DP-B containing 10.5% martensite. Grains colored in white are martensite, and remaining brown grains are ferrite.</figcaption></figure>

### Preparing tiny round-bar specimen

We subjected No. 5 type specimens, 2.5-mm thick and standardized by Japanese Industrial Standards (JIS) Z2241, to preliminary tensile testing before subsequent compact tensile tests. These specimens were elongated by up to 23.14% and 19.44% for DP-A and DP-B, respectively,
surpassing their uniform elongations of 19% and 16%, as indicated in **Table 2**. Prior to the compact tensile tests, which involved real-time diameter measurements, the specimens experienced necking due to deformations exceeding their uniform elongation.

As depicted in **Fig. 2**, the preliminary tensile test was aligned such that the tensile direction of the tiny round-bar specimens corresponded to the transverse direction of the cold-rolling process (referred to as the transverse direction). Accordingly, the JIS No. 5 specimens were tensioned in the transverse direction without a strain path change (**Fig. 2a**) and in the rolling direction in the presence of an orthogonal strain path change (**Fig. 2b**).

The design of the tiny round-bar specimens is presented in **Fig. 3**. The thinnest part of these specimens had a diameter of 1.0 mm, enabling their extraction from the JIS No. 5 specimens that had undergone pre-straining exceeding the onset of necking.

<figure><img src="image_1025.png" /><figcaption>Fig. 2 Sampling position of the tiny round-bar specimens and the tensile direction of the JIS No. 5 tensile test prior to compact tensile testing: (a) without subsequent strain path change and (b) with a subsequent orthogonal strain path change.</figcaption></figure>

<figure><img src="image_1026.png" /><figcaption>Fig. 3 Geometry of the tiny round-bar specimen depicted in an axisymmetric two-dimensional drawing.</figcaption></figure>

### Compact tensile test with real-time diameter measurements

After extracting the tiny round-bar specimens, compact tensile tests with real-time diameter measurements were executed, generating true stress versus cross-sectional reduction area (CSRA,
surpassing their uniform elongations of 19% and 16%, as indicated in **Table 2**. Prior to the compact tensile tests, which involved real-time diameter measurements, the specimens experienced necking due to deformations exceeding their uniform elongation.

As depicted in **Fig. 2**, the preliminary tensile test was aligned such that the tensile direction of the tiny round-bar specimens corresponded to the transverse direction of the cold-rolling process (referred to as the transverse direction). Accordingly, the JIS No. 5 specimens were tensioned in the transverse direction without a strain path change (**Fig. 2a**) and in the rolling direction in the presence of an orthogonal strain path change (**Fig. 2b**).

The design of the tiny round-bar specimens is presented in **Fig. 3**. The thinnest part of these specimens had a diameter of 1.0 mm, enabling their extraction from the JIS No. 5 specimens that had undergone pre-straining exceeding the onset of necking.

<figure><img src="image_1030.png" /><figcaption>Fig. 2 Sampling position of the tiny round-bar specimens and the tensile direction of the JIS No. 5 tensile test prior to compact tensile testing: (a) without subsequent strain path change and (b) with a subsequent orthogonal strain path change.</figcaption></figure>

<figure><img src="image_1031.png" /><figcaption>Fig. 3 Geometry of the tiny round-bar specimen depicted in an axisymmetric two-dimensional drawing.</figcaption></figure>

### Compact tensile test with real-time diameter measurements

After extracting the tiny round-bar specimens, compact tensile tests with real-time diameter measurements were executed, generating true stress versus cross-sectional reduction area (CSRA,
$\rho$) curves. In these compact tensile tests, the diameters, $D_a$ and $D_b$, of the thinnest cross-sectional area under tension were continuously measured in two directions. The specimen was oriented such that $D_b$ corresponded to the thickness of the original steel material. The true stress $\sigma_t$ and $\rho$ in the cross-section were calculated as follows:

$$ \sigma_t = \frac{4F}{\pi D_a D_b}, \qquad (1) $$

$$ \rho = \ln \left( \frac{D_0^2}{D_a D_b} \right), \qquad (2) $$

where $D_0$ corresponds to the initial diameter of the specimen (1.0 mm).

## Results and discussion

**Figure 4** illustrates the $\sigma_t - \rho$ curves measured for DP-A and DP-B. In the materials subjected to orthogonal abrupt strain path changes (referred to as the “orthogonal-loading case”), an increased initial yield stress was observed compared to those deformed in the same direction as the pre-strain (referred to as the “same-direction case”), although their subsequent $\sigma_t$ decreased. Notably, the same-direction and orthogonal-loading cases reached peak loads at $\rho = 0.01$ and $\rho = 0.03$, respectively, marking the onset of necking in the specimens.

In DP-A (Fig. 4a), at the marked final point on the plot, the cross-sectional reduction area (CSRA) $\rho$ was 0.92 for the orthogonal-loading case, compared to $\rho = 0.87$ for the same-direction case, indicating a marginal enhancement in the fracture strain. Given the challenges in precisely determining the fracture timing in the $\sigma_t - \rho$ curves, the definitive $\rho$ was deduced from scanning electron microscope (SEM) images of the fracture surface (Fig. 5). The fracture surface areas were quantified from these images. As a result, a comparative analysis between the same direction (Fig. 5a) and orthogonal-loading cases (Fig. 5b) also exhibited a 0.09 increase in $\rho$. Previous studies have suggested that in multi-phase steels, an abrupt strain path change can relocate sites of microscopic plastic deformation, which is considered to alleviate deformation localization during pre-straining and thereby augment ductility. This study hypothesized a similar phenomenon in the DP steels, albeit with a relatively minor effect. Variations in work hardening and resultant stress states may obscure microscopic enhancement in ductility.

In the case of DP-B, the impact of a strain path change on the $\sigma_t - \rho$ curves was less significant than DP-A (Fig. 4b). The improvement stemming from the strain-path change was modest but contributed to increased fracture strain. Final $\rho$ values on the $\sigma_t - \rho$ curve were 0.55 and 0.58 for the same-direction and orthogonal-loading cases, respectively. SEM images of the fracture surfaces (Figs. 5c and 5d) unveiled a 0.1 increase in the actual final $\rho$, mirroring the observation in DP-A. It suggests that the martensite volume fraction does not markedly influence the effect of strain path change on the increase in fracture strain.

An additional consideration pertains to the disparity in anisotropy. Fig. 6 illustrates the history of $D_b/D_a$ relative to $\rho$. Lower values signify a flattened elliptical shape in the thickness direction. The behaviors observed in DP-A (Fig. 6a) and DP-B (Fig. 6b) are akin. Therefore, we focus on DP-A to discuss its anisotropy: For $\rho$ lower than approximately 0.45, $D_b/D_a$ was marginally smaller for the same-direction cases than for the orthogonal-loading cases. The relationship started to invert, commencing from approximately $\rho = 0.45$ with an intensified flattening in the thickness direction in the abrupt strain path change of materials. In orthogonal-loading cases, no flattening occurred until $\rho$ reached approximately 0.08, after which $D_b/D_a$ decreased sharply. This flattening of the cross-section indicates more pronounced necking in the section involving $D_b$, potentially leading to a transformation into a geometry with higher triaxiality, which might diminish the impact on the enhanced fracture strain.
$\rho$) curves. In these compact tensile tests, the diameters, $D_a$ and $D_b$, of the thinnest cross-sectional area under tension were continuously measured in two directions. The specimen was oriented such that $D_b$ corresponded to the thickness of the original steel material. The true stress $\sigma_t$ and $\rho$ in the cross-section were calculated as follows:

$$ \sigma_t = \frac{4F}{\pi D_a D_b}, \qquad (1) $$

$$ \rho = \ln \left( \frac{D_0^2}{D_a D_b} \right), \qquad (2) $$

where $D_0$ corresponds to the initial diameter of the specimen (1.0 mm).

## Results and discussion

**Figure 4** illustrates the $\sigma_t - \rho$ curves measured for DP-A and DP-B. In the materials subjected to orthogonal abrupt strain path changes (referred to as the “orthogonal-loading case”), an increased initial yield stress was observed compared to those deformed in the same direction as the pre-strain (referred to as the “same-direction case”), although their subsequent $\sigma_t$ decreased. Notably, the same-direction and orthogonal-loading cases reached peak loads at $\rho = 0.01$ and $\rho = 0.03$, respectively, marking the onset of necking in the specimens.

In DP-A (Fig. 4a), at the marked final point on the plot, the cross-sectional reduction area (CSRA) $\rho$ was 0.92 for the orthogonal-loading case, compared to $\rho = 0.87$ for the same-direction case, indicating a marginal enhancement in the fracture strain. Given the challenges in precisely determining the fracture timing in the $\sigma_t - \rho$ curves, the definitive $\rho$ was deduced from scanning electron microscope (SEM) images of the fracture surface (Fig. 5). The fracture surface areas were quantified from these images. As a result, a comparative analysis between the same direction (Fig. 5a) and orthogonal-loading cases (Fig. 5b) also exhibited a 0.09 increase in $\rho$. Previous studies have suggested that in multi-phase steels, an abrupt strain path change can relocate sites of microscopic plastic deformation, which is considered to alleviate deformation localization during pre-straining and thereby augment ductility. This study hypothesized a similar phenomenon in the DP steels, albeit with a relatively minor effect. Variations in work hardening and resultant stress states may obscure microscopic enhancement in ductility.

In the case of DP-B, the impact of a strain path change on the $\sigma_t - \rho$ curves was less significant than DP-A (Fig. 4b). The improvement stemming from the strain-path change was modest but contributed to increased fracture strain. Final $\rho$ values on the $\sigma_t - \rho$ curve were 0.55 and 0.58 for the same-direction and orthogonal-loading cases, respectively. SEM images of the fracture surfaces (Figs. 5c and 5d) unveiled a 0.1 increase in the actual final $\rho$, mirroring the observation in DP-A. It suggests that the martensite volume fraction does not markedly influence the effect of strain path change on the increase in fracture strain.

An additional consideration pertains to the disparity in anisotropy. Fig. 6 illustrates the history of $D_b/D_a$ relative to $\rho$. Lower values signify a flattened elliptical shape in the thickness direction. The behaviors observed in DP-A (Fig. 6a) and DP-B (Fig. 6b) are akin. Therefore, we focus on DP-A to discuss its anisotropy: For $\rho$ lower than approximately 0.45, $D_b/D_a$ was marginally smaller for the same-direction cases than for the orthogonal-loading cases. The relationship started to invert, commencing from approximately $\rho = 0.45$ with an intensified flattening in the thickness direction in the abrupt strain path change of materials. In orthogonal-loading cases, no flattening occurred until $\rho$ reached approximately 0.08, after which $D_b/D_a$ decreased sharply. This flattening of the cross-section indicates more pronounced necking in the section involving $D_b$, potentially leading to a transformation into a geometry with higher triaxiality, which might diminish the impact on the enhanced fracture strain.
<figure><img src="image_1027.png" /></figure>

<figure><img src="image_1028.png" /><figcaption>Fig. 4 True stress $\sigma_t$ versus cross-sectional reduction area (CSRA) $\rho$: (a) DP-A and (b) DP-B.</figcaption></figure>

<figure><img src="image_1029.png" /><figcaption>Fig. 5 SEM images of fracture surface: (a) same-direction case for DP-A ($\rho$ = 0.95), (b) orthogonal-loading case for DP-A ($\rho$ = 1.04), (c) same-direction case for DP-B ($\rho$ = 0.60), and (d) orthogonal-loading case for DP-B($\rho$ = 0.70).</figcaption></figure>
<figure><img src="image_1032.png" /></figure>

<figure><img src="image_1033.png" /><figcaption>Fig. 4 True stress $\sigma_t$ versus cross-sectional reduction area (CSRA) $\rho$: (a) DP-A and (b) DP-B.</figcaption></figure>

<figure><img src="image_1034.png" /><figcaption>Fig. 5 SEM images of fracture surface: (a) same-direction case for DP-A ($\rho$ = 0.95), (b) orthogonal-loading case for DP-A ($\rho$ = 1.04), (c) same-direction case for DP-B ($\rho$ = 0.60), and (d) orthogonal-loading case for DP-B($\rho$ = 0.70).</figcaption></figure>
<figure><img src="image_1030.png" /><figcaption>Fig. 6 D<sub>b</sub>/D<sub>a</sub> versus CSRA ρ: (a) DP-A and (b) DP-B.</figcaption></figure>

**Summary**

This study investigated the impact of enhancing fracture strain during an abrupt orthogonal strain path change in two types of DP steels with different martensite volume fractions. The results validated an improvement in fracture strains for both types of DP steels, with enhancements of 0.1. While this quantitative behavior differs from that observed in ferrite-pearlite steels, it may significantly contribute to press-formability in processes like stretch-flange forming, where a difference of 0.1 in fracture strain holds substantial importance. We postulated that the degree of improvement in fracture strain due to the abrupt strain path change was influenced not only by the previously reported microscopic deformation concentration mechanisms but also by variations in anisotropy and work-hardening behavior.

**References**

[1] T. Yasutomi, T. Matsuno, E. Sakurada, S. Yonemura, H. Shoji, M. Ohata, Ductile Fracture Behavior for Two-phase Steel under Tensile Deformation with Strain Path Change, Tetsu-to-Hagane, 108(2022) 772-783(In Japanese). https://doi.org/10.2355/tetsutohagane.TETSU-2022-020

[2] T. Matsuno, D. Kondo, T. Hama, T. Naito, Y. Okitsu, S. Hayashi, K. Takada, Flow stress curves for 980MPa- and 1.5GPa-class ultra-high-strength steel sheets weakened under high-stress triaxiality, Int. J. Mech. Sci., 261(2024), 108671. https://doi.org/10.1016/j.ijmecsci.2023.108671

[3] T. Matsuno, T. Fujita, T. Matsuda, T. Shibayama, T. Hojo, I. Watanabe, Unstable stress-triaxiality development and contrasting weakening in two types of high-strength transformation-induced plasticity (TRIP) steels: Insights from a new compact tensile testing method, J. Mater. Process. Technol., 322(2023), 118174. https://doi.org/10.1016/j.jmatprotec.2023.118174

[4] T. Matsuno, T. Hojo, I. Watanabe, A. Shiro, T. Shobu, K. Kajiwara, Tensile deformation behavior of TRIP-aided bainitic ferrite steel in the post-necking strain region, Sci. Tech. Adv. Mater.: Methods, 1(2021), 56-74. https://doi.org/10.1080/27660400.2021.1922207
<figure><img src="image_1035.png" /><figcaption>Fig. 6 D<sub>b</sub>/D<sub>a</sub> versus CSRA ρ: (a) DP-A and (b) DP-B.</figcaption></figure>

**Summary**

This study investigated the impact of enhancing fracture strain during an abrupt orthogonal strain path change in two types of DP steels with different martensite volume fractions. The results validated an improvement in fracture strains for both types of DP steels, with enhancements of 0.1. While this quantitative behavior differs from that observed in ferrite-pearlite steels, it may significantly contribute to press-formability in processes like stretch-flange forming, where a difference of 0.1 in fracture strain holds substantial importance. We postulated that the degree of improvement in fracture strain due to the abrupt strain path change was influenced not only by the previously reported microscopic deformation concentration mechanisms but also by variations in anisotropy and work-hardening behavior.

**References**

[1] T. Yasutomi, T. Matsuno, E. Sakurada, S. Yonemura, H. Shoji, M. Ohata, Ductile Fracture Behavior for Two-phase Steel under Tensile Deformation with Strain Path Change, Tetsu-to-Hagane, 108(2022) 772-783 (In Japanese). https://doi.org/10.2355/tetsutohagane.TETSU-2022-020

[2] T. Matsuno, D. Kondo, T. Hama, T. Naito, Y. Okitsu, S. Hayashi, K. Takada, Flow stress curves for 980MPa- and 1.5GPa-class ultra-high-strength steel sheets weakened under high-stress triaxiality, Int. J. Mech. Sci., 261(2024), 108671. https://doi.org/10.1016/j.ijmecsci.2023.108671

[3] T. Matsuno, T. Fujita, T. Matsuda, T. Shibayama, T. Hojo, I. Watanabe, Unstable stress-triaxiality development and contrasting weakening in two types of high-strength transformation-induced plasticity (TRIP) steels: Insights from a new compact tensile testing method, J. Mater. Process. Technol., 322(2023), 118174. https://doi.org/10.1016/j.jmatprotec.2023.118174

[4] T. Matsuno, T. Hojo, I. Watanabe, A. Shiro, T. Shobu, K. Kajiwara, Tensile deformation behavior of TRIP-aided bainitic ferrite steel in the post-necking strain region, Sci. Tech. Adv. Mater.: Methods, 1(2021), 56-74. https://doi.org/10.1080/27660400.2021.1922207
# Numerical assessment of the role of anisotropy on strain localization in uniaxial tension

CORALLO Luca<fnref n="1" />,ᵃ*, CAZACU Oana<fnref n="2" />,ᵇ and VERLEYSEN Patricia<fnref n="1" />,ᶜ

**Keywords:** Uniaxial Tensile Test, Necking Bands, Orthotropy

**Abstract.** Despite the experimental evidence of the effect of loading direction on the inclination of necking bands in uniaxial tension, theoretical and numerical studies are mainly devoted to isotropic materials. Recently, theoretical, numerical, and experimental works have put into evidence the key role of the material anisotropy on the formation of localized necking bands. In particular, Cazacu and Rodriguez (2019) [1] provided analytical expressions for the orientations of the necking bands that develop under uniaxial tension in flat specimens. It was shown that there is a switch in the orientation of main necking band from acute to obtuse, which is correlated with the anisotropy in yield stresses. In this paper, we conduct a FE study on virtual materials with the same anisotropy in Lankford coefficients but different anisotropy in yield stresses. We show that, although both materials have a very slight anisotropy in yield stresses, it strongly affects the localization behavior.

**Introduction**

Early experimental studies have put into evidence strain localization in flat tensile metallic specimens subjected to quasi-static uniaxial tension [2,3]. For isotropic elastic-plastic metallic materials, Neilsen and Schreyer [4] derived analytical expressions showing that under uniaxial tension there are two localization bands which are equally inclined with respect to the loading axis, the inclination being of ±54.74°. On the other hand, in an early work of 1928, Korber and Hoff [5] reported the influence of the loading direction on the inclinations of the necking bands in sheet metals subjected to uniaxial tension. In 1948, Hill [6] showed that bands of equal inclination to the loading axis develop in the rolling direction (RD) and transverse direction (TD). In the same work, Hill provided analytic expressions for the bands inclinations in terms of the anisotropy coefficients of his model showing that, it might exist an off-axis specimen orientation for which bands of equal inclination may develop. The recent work of Cazacu and co-workers [1,7,8] further elucidated the role of material anisotropy on the inclination of localization bands in flat tensile specimens. In [1], Cazacu and Rodríguez-Martínez conducted a theoretical investigation on the effect of the loading orientation on the inclinations of localization bands in flat tensile specimens. It was demonstrated that, if the yield function and plastic potential of the material are known, the bands angles can be easily obtained by solving a second-order algebraic equation. Moreover, the authors formulated an explicit equation expressed in terms of Lankford coefficient to calculate the inclination of the necking band for the orientation corresponding to a minimum or maximum of the function σ(θ) in the plane of the sheet. As concerns numerical modeling of strain localization, beginning with the pioneering finite-element (FE) study of Tvergaard (1993) [9], to trigger the development of necking bands, geometrical or material imperfections are often introduced (see also [10,11]). Moreover, the influence of voids present in a metallic material on the development of instabilities

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

<footnote n="1">¹Department of Electromechanical, Systems and Metal Engineering, MST-DyMaLab Research Group, Ghent University, Technologiepark 46, 9052 Zwijnaarde, Belgium</footnote>

<footnote n="2">²Department of Materials Sciences and Engineering, University of Arizona, 1235 E. James E. Rogers Ways, Tucson, AZ 85721, USA</footnote>

<footnote>ᵃLuca.Corallo@Unict.it, ᵇoanacazacu@arizona.edu, ᶜPatricia.Verleysen@UGent.be</footnote>
# Numerical assessment of the role of anisotropy on strain localization in uniaxial tension

CORALLO Luca<fnref n="1" />,ᵃ*, CAZACU Oana<fnref n="2" />,ᵇ and VERLEYSEN Patricia<fnref n="1" />,ᶜ

**Keywords:** Uniaxial Tensile Test, Necking Bands, Orthotropy

**Abstract.** Despite the experimental evidence of the effect of loading direction on the inclination of necking bands in uniaxial tension, theoretical and numerical studies are mainly devoted to isotropic materials. Recently, theoretical, numerical, and experimental works have put into evidence the key role of the material anisotropy on the formation of localized necking bands. In particular, Cazacu and Rodriguez (2019) [1] provided analytical expressions for the orientations of the necking bands that develop under uniaxial tension in flat specimens. It was shown that there is a switch in the orientation of main necking band from acute to obtuse, which is correlated with the anisotropy in yield stresses. In this paper, we conduct a FE study on virtual materials with the same anisotropy in Lankford coefficients but different anisotropy in yield stresses. We show that, although both materials have a very slight anisotropy in yield stresses, it strongly affects the localization behavior.

**Introduction**

Early experimental studies have put into evidence strain localization in flat tensile metallic specimens subjected to quasi-static uniaxial tension [2,3]. For isotropic elastic-plastic metallic materials, Neilsen and Schreyer [4] derived analytical expressions showing that under uniaxial tension there are two localization bands which are equally inclined with respect to the loading axis, the inclination being of ±54.74°. On the other hand, in an early work of 1928, Korber and Hoff [5] reported the influence of the loading direction on the inclinations of the necking bands in sheet metals subjected to uniaxial tension. In 1948, Hill [6] showed that bands of equal inclination to the loading axis develop in the rolling direction (RD) and transverse direction (TD). In the same work, Hill provided analytic expressions for the bands inclinations in terms of the anisotropy coefficients of his model showing that, it might exist an off-axis specimen orientation for which bands of equal inclination may develop. The recent work of Cazacu and co-workers [1,7,8] further elucidated the role of material anisotropy on the inclination of localization bands in flat tensile specimens. In [1], Cazacu and Rodríguez-Martínez conducted a theoretical investigation on the effect of the loading orientation on the inclinations of localization bands in flat tensile specimens. It was demonstrated that, if the yield function and plastic potential of the material are known, the bands angles can be easily obtained by solving a second-order algebraic equation. Moreover, the authors formulated an explicit equation expressed in terms of Lankford coefficient to calculate the inclination of the necking band for the orientation corresponding to a minimum or maximum of the function σ(θ) in the plane of the sheet. As concerns numerical modeling of strain localization, beginning with the pioneering finite-element (FE) study of Tvergaard (1993) [9], to trigger the development of necking bands, geometrical or material imperfections are often introduced (see also [10,11]). Moreover, the influence of voids present in a metallic material on the development of instabilities

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

<footnote n="1">¹Department of Electromechanical, Systems and Metal Engineering, MST-DyMaLab Research Group, Ghent University, Technologiepark 46, 9052 Zwijnaarde, Belgium</footnote>

<footnote n="2">²Department of Materials Sciences and Engineering, University of Arizona, 1235 E. James E. Rogers Ways, Tucson, AZ 85721, USA</footnote>

<footnote>ᵃLuca.Corallo@Unict.it, ᵇoanacazacu@arizona.edu, ᶜPatricia.Verleysen@UGent.be</footnote>
has been extensively studied using both micromechanical bifurcation analyses and models accounting for dilatational plasticity [12–14]. All the aforementioned studies concern experimental, theoretical, and numerical investigations of strain localization phenomena in isotropic materials. In [8], Corallo et al. conducted a comprehensive experimental, theoretical, and FE study of the effect of plastic anisotropy in Ti6Al4V deformed at high strain rates and its effects on strain localization. Full-field strain measurements revealed a strong influence of the loading direction on the inclination of the necking bands. The authors were able to capture the main experimental trends, i.e., sign reversal of the inclination of the necking bands with respect to the loading direction, strain fields at the sample surface, minimum sample width, by modelling the anisotropic response of Ti6Al4V using a constitutive model with yielding described by the CB2001 [15] orthotropic yield criterion. In this work, we further elucidate the role of both r-value and flow stress anisotropy on the FE predictions of necking bands for two virtual materials obtained varying the anisotropic parameters of Ti6Al4V.

## Theoretical analysis of necking bands

The main results of the theoretical analysis of necking bands in flat tensile specimens proposed by Cazacu and Rodríguez-Martínez [1] are presented in the following. Let us denote by θ the angle between the loading direction and the rolling direction, r(θ) the Lankford coefficient, β the necking band inclination. Starting from the fact that under uniaxial tension, the localization bands are along the directions of zero extension, the authors showed that irrespective of the mathematical form of the yield function and plastic potential used to model the behavior of an orthotropic material, bands of equal inclination with respect to the loading axis occur only for the specimens of orientation indicated with θ*, that corresponds to either maximum or minimum in the variation of the yield stress with the angle to the rolling direction and

$$ \beta(\theta^*) = \pm \tan^{-1} \left( \sqrt{\frac{1+r(\theta^*)}{r(\theta^*)}} \right) \quad (1) $$

For all the other specimen orientations, necking bands of different inclinations occur and these inclinations are the solutions of a second-order algebraic equation in u = tanβ. For associated plasticity, this equation is given by:

$$ A \left( u^2 - \frac{1+r(\theta)}{r(\theta)} \right) + 2Bu = 0 \quad (2) $$

$$ A = \frac{\partial \bar{\sigma}}{\partial \sigma_{xx}} \sin^2 \theta - \frac{\partial \bar{\sigma}}{\partial \sigma_{xy}} \sin(2\theta) + \frac{\partial \bar{\sigma}}{\partial \sigma_{yy}} \cos^2 \theta \quad (3) $$

$$ B = \sin\theta\cos\theta \left( \frac{\partial \bar{\sigma}}{\partial \sigma_{yy}} - \frac{\partial \bar{\sigma}}{\partial \sigma_{xx}} \right) + \cos(2\theta) \frac{\partial \bar{\sigma}}{\partial \sigma_{xy}} \quad (4) $$

where $\bar{\sigma}$ is the effective stress associated with the orthotropic yield function, and Oxyz is the coordinate system associated with the orthotropy axes, with x being along the rolling direction (RD) and y along the transverse direction (TD).

## Materials and constitutive model

The role of the material anisotropy on the formation of necking bands in uniaxial tension is investigated for two virtual materials having mild anisotropy in yield stresses and pronounced anisotropy in r-values, a behaviour which is typical for titanium materials [8,16–18]. Anisotropic
has been extensively studied using both micromechanical bifurcation analyses and models accounting for dilatational plasticity [12–14]. All the aforementioned studies concern experimental, theoretical, and numerical investigations of strain localization phenomena in isotropic materials. In [8], Corallo et al. conducted a comprehensive experimental, theoretical, and FE study of the effect of plastic anisotropy in Ti6Al4V deformed at high strain rates and its effects on strain localization. Full-field strain measurements revealed a strong influence of the loading direction on the inclination of the necking bands. The authors were able to capture the main experimental trends, i.e., sign reversal of the inclination of the necking bands with respect to the loading direction, strain fields at the sample surface, minimum sample width, by modelling the anisotropic response of Ti6Al4V using a constitutive model with yielding describe by the CB2001 [15] orthotropic yield criterion. In this work, we further elucidate the role of both r-value and flow stress anisotropy on the FE predictions of necking bands for two virtual materials obtained varying the anisotropic parameters of Ti6Al4V.

## Theoretical analysis of necking bands

The main results of the theoretical analysis of necking bands in flat tensile specimens proposed by Cazacu and Rodríguez-Martínez [1] are presented in the following. Let us denote by θ the angle between the loading direction and the rolling direction, r(θ) the Lankford coefficient, β the necking band inclination. Starting from the fact that under uniaxial tension, the localization bands are along the directions of zero extension, the authors showed that irrespective of the mathematical form of the yield function and plastic potential used to model the behavior of an orthotropic material, bands of equal inclination with respect to the loading axis occur only for the specimens of orientation indicated with θ*, that corresponds to either maximum or minimum in the variation of the yield stress with the angle to the rolling direction and

$$ \beta(\theta^*) = \pm \tan^{-1} \left( \sqrt{\frac{1+r(\theta^*)}{r(\theta^*)}} \right) \quad (1) $$

For all the other specimen orientations, necking bands of different inclinations occur and these inclinations are the solutions of a second-order algebraic equation in u = tanβ. For associated plasticity, this equation is given by:

$$ A \left( u^2 - \frac{1+r(\theta)}{r(\theta)} \right) + 2Bu = 0 \quad (2) $$

$$ A = \frac{\partial \bar{\sigma}}{\partial \sigma_{xx}} \sin^2 \theta - \frac{\partial \bar{\sigma}}{\partial \sigma_{xy}} \sin(2\theta) + \frac{\partial \bar{\sigma}}{\partial \sigma_{yy}} \cos^2 \theta \quad (3) $$

$$ B = \sin\theta\cos\theta \left( \frac{\partial \bar{\sigma}}{\partial \sigma_{yy}} - \frac{\partial \bar{\sigma}}{\partial \sigma_{xx}} \right) + \cos(2\theta) \frac{\partial \bar{\sigma}}{\partial \sigma_{xy}} \quad (4) $$

where $\bar{\sigma}$ is the effective stress associated with the orthotropic yield function, and Oxyz is the coordinate system associated with the orthotropy axes, with x being along the rolling direction (RD) and y along the transverse direction (TD).

## Materials and constitutive model

The role of the material anisotropy on the formation of necking bands in uniaxial tension is investigated for two virtual materials having mild anisotropy in yield stresses and pronounced anisotropy in r-values, a behaviour which is typical for titanium materials [8,16–18]. Anisotropic
properties, i.e., r-value and yield stress, as a function of θ for the two materials are compiled in Table 1.

<table><caption>Table 1. r-values and normalized flow stresses</caption><thead><tr><th></th><th></th><th>0°</th><th>15°</th><th>30°</th><th>45°</th><th>60°</th><th>75°</th><th>90°</th></tr></thead><tbody><tr><td rowspan="2">Material 1</td><td>r-value</td><td>2.50</td><td>2.30</td><td>1.85</td><td>1.27</td><td>1.01</td><td>1.16</td><td>1.23</td></tr><tr><td>σ<sub>θ</sub> / σ<sub>0</sub></td><td>1.00</td><td>0.99</td><td>0.95</td><td>0.92</td><td>0.97</td><td>0.98</td><td>1.01</td></tr><tr><td rowspan="2">Material 2</td><td>r-value</td><td>2.50</td><td>2.30</td><td>1.85</td><td>1.27</td><td>1.01</td><td>1.16</td><td>1.23</td></tr><tr><td>σ<sub>θ</sub> / σ<sub>0</sub></td><td>1.00</td><td>1.01</td><td>1.05</td><td>1.07</td><td>1.03</td><td>1.01</td><td>0.99</td></tr></tbody></table>

Note that both materials have same anisotropy in r-values with a minimum at $r(θ) = 60°$. However, they display different trends in sigma variation $σ_θ / σ_0$ vs. $θ$: for $θ = 45°$ the Material 1 presents a minimum, while the Material 2 has a maximum. To model the anisotropic behaviour of the materials listed in Table 1, a constitutive elastic-plastic model with yielding described by Cazacu and Barlat orthotropic criterion [15] is adopted. To simplify writing, this criterion will be designated as CB2001. The constitutive assumptions and model equations are presented in the following. The objective derivative of the Cauchy stress tensor $\sigma$ is calculated using:

$$ \dot{\sigma} = C^e: (D - D^p) \qquad (5) $$

in which $C^e$ is the isotropic elastic tensor, $D$ and $D^p$ are the total and plastic part of the rate of the deformation tensor, respectively. The double contracted product between two tensors is denoted with “:”. Moreover, an associated flow rule is assumed, i.e., the yield function and plastic potential coincide, resulting in:

$$ D^p = \dot{\lambda} \frac{\partial \varphi}{\partial} \qquad (6) $$

with $\varphi$ being the yield function. The orthotropic yield function CB2001 is expressed as [15]:

$$ \varphi(\cdot) = (J_2^0)^3 - c \cdot (J_3^0)^2 \qquad (7) $$

where $J_2^0$ and $J_3^0$ are the orthotropic second and third generalized invariants of the stress deviator tensor, respectively, and c is a material parameter (see also [15]). The effective stress $\bar{\sigma}$, associated to the orthotropic yield function $\varphi(\sigma)$ is:

$$ \bar{\sigma} = \tilde{A} \cdot [(J_2^0)^3 - c \cdot (J_3^0)^2]^{1/6} \qquad (8) $$

The constant $\tilde{A}$ is determined in such a way that for uniaxial tension $\bar{\sigma}$ reduces to the yield stress along RD. The anisotropic coefficients $a_k$ ($k = 1, 2...6$), $b_j$ ($j = 1, 2...11$) and c, associated to the anisotropic criterion are identified using the uniaxial tensile data listed in Table 1. The anisotropic constants are identified by minimizing the following cost function (e.g., see [19]):
properties, i.e., r-value and yield stress, as a function of θ for the two materials are compiled in Table 1.

<table><caption>Table 1. r-values and normalized flow stresses</caption><thead><tr><th></th><th></th><th>0°</th><th>15°</th><th>30°</th><th>45°</th><th>60°</th><th>75°</th><th>90°</th></tr></thead><tbody><tr><td rowspan="2">Material 1</td><td>r-value</td><td>2.50</td><td>2.30</td><td>1.85</td><td>1.27</td><td>1.01</td><td>1.16</td><td>1.23</td></tr><tr><td>σ<sub>θ</sub> / σ<sub>0</sub></td><td>1.00</td><td>0.99</td><td>0.95</td><td>0.92</td><td>0.97</td><td>0.98</td><td>1.01</td></tr><tr><td rowspan="2">Material 2</td><td>r-value</td><td>2.50</td><td>2.30</td><td>1.85</td><td>1.27</td><td>1.01</td><td>1.16</td><td>1.23</td></tr><tr><td>σ<sub>θ</sub> / σ<sub>0</sub></td><td>1.00</td><td>1.01</td><td>1.05</td><td>1.07</td><td>1.03</td><td>1.01</td><td>0.99</td></tr></tbody></table>

Note that both materials have same anisotropy in r-values with a minimum at $r(θ) = 60°$. However, they display different trends in sigma variation $σ_θ / σ_0$ vs. $θ$: for $θ = 45°$ the Material 1 presents a minimum, while the Material 2 has a maximum. To model the anisotropic behaviour of the materials listed in Table 1, a constitutive elastic-plastic model with yielding described by Cazacu and Barlat orthotropic criterion [15] is adopted. To simplify writing, this criterion will be designated as CB2001. The constitutive assumptions and model equations are presented in the following. The objective derivative of the Cauchy stress tensor $\sigma$ is calculated using:

$$ \dot{\sigma} = C^e: (D - D^p) \qquad (5) $$

in which $C^e$ is the isotropic elastic tensor, $D$ and $D^p$ are the total and plastic part of the rate of the deformation tensor, respectively. The double contracted product between two tensors is denoted with “:”. Moreover, an associated flow rule is assumed, i.e., the yield function and plastic potential coincide, resulting in:

$$ D^p = \dot{\lambda} \frac{\partial \varphi}{\partial} \qquad (6) $$

with $\varphi$ being the yield function. The orthotropic yield function CB2001 is expressed as [15]:

$$ \varphi(\cdot) = (J_2^0)^3 - c \cdot (J_3^0)^2 \qquad (7) $$

where $J_2^0$ and $J_3^0$ are the orthotropic second and third generalized invariants of the stress deviator tensor, respectively, and c is a material parameter (see also [15]). The effective stress $\bar{\sigma}$, associated to the orthotropic yield function $\varphi(\sigma)$ is:

$$ \bar{\sigma} = \tilde{A} \cdot [(J_2^0)^3 - c \cdot (J_3^0)^2]^{1/6} \qquad (8) $$

The constant $\tilde{A}$ is determined in such a way that for uniaxial tension $\bar{\sigma}$ reduces to the yield stress along RD. The anisotropic coefficients $a_k$ ($k = 1, 2...6$), $b_j$ ($j = 1, 2...11$) and c, associated to the anisotropic criterion are identified using the uniaxial tensile data listed in Table 1. The anisotropic constants are identified by minimizing the following cost function (e.g., see [19]):
$$E(a_k, b_j, c) = \sum_m w_m \left( \frac{\sigma_{\theta}^{CB2001}}{\sigma_{\theta}^{Data}} - 1 \right)^2 + \sum_n v_n \left( \frac{r_{\theta}^{CB2001}}{r_{\theta}^{Data}} - 1 \right)^2 \quad (9)$$

In Eq. (9) the range of the indices *m*, *n* is given by the number of tensile yield stresses and Lankford coefficients available while *w*<sub>*m*</sub>, *v*<sub>*n*</sub> are weights given to the respective experimental values. The superscripts *CB2001* and *Data* indicate the theoretical values of yield stress or Lankford coefficient obtained with the anisotropic model, and the data from Table 1, respectively. Since for sheet materials tests involving out-of-plane stress components are very difficult or non-feasible, in the 3-D FE calculations, isotropic values are assigned to the anisotropy coefficients associated to the out-of-plane stress components, namely: *a*<sub>5</sub>=*a*<sub>6</sub>=*b*<sub>6</sub>=*b*<sub>7</sub>=*b*<sub>8</sub>=*b*<sub>9</sub>=*b*<sub>11</sub>=1 (e.g., [16]). Table 2 compiles the numerical values of the anisotropic parameters obtained from the identification procedure.

<table><caption>Table 2. Parameters CB2001 for all the materials listed Table 3.</caption><thead><tr><th></th><th>a1</th><th>a2</th><th>a3</th><th>a4</th><th>c</th><th>b1</th><th>b2</th><th>b3</th><th>b4</th><th>b5</th><th>b10</th><th>c</th></tr></thead><tbody><tr><td>Material 1</td><td>0.87</td><td>0.87</td><td>0.36</td><td>0.74</td><td>0.03</td><td>1.06</td><td>-1.95</td><td>-10.0</td><td>-9.11</td><td>-4.01</td><td>-0.25</td><td>0.03</td></tr><tr><td>Material 2</td><td>0.58</td><td>0.46</td><td>0.61</td><td>0.48</td><td>0.38</td><td>0.20</td><td>2.13</td><td>-0.02</td><td>0.24</td><td>0.74</td><td>0.63</td><td>0.38</td></tr></tbody></table>

<figure><img src="image_1032.png" /><figcaption>Fig. 1. Prediction of (a) Lankford coefficients, (b) normalized flow stresses and (c) projection in the biaxial ($\sigma_{xx} - \sigma_{yy}$) plane (with $\sigma_{xy} = 0$) of the yield locus according to CB2001 for the two materials.</figcaption></figure>

Fig. 1a and Fig. 1b show the predicted Lankford coefficients and normalized flow stresses according to CB2001 for both materials in Table 1. Moreover, in Fig. 1c, the projections in the biaxial plane ($\sigma_{xx} - \sigma_{yy}$) plane (with $\sigma_{xy} = 0$) of the yield loci are given. The materials hardening
$$E(a_k, b_j, c) = \sum_m w_m \left( \frac{\sigma_{\theta}^{CB2001}}{\sigma_{\theta}^{Data}} - 1 \right)^2 + \sum_n v_n \left( \frac{r_{\theta}^{CB2001}}{r_{\theta}^{Data}} - 1 \right)^2 \quad (9)$$

In Eq. (9) the range of the indices *m*, *n* is given by the number of tensile yield stresses and Lankford coefficients available while *w*<sub>*m*</sub>, *v*<sub>*n*</sub> are weights given to the respective experimental values. The superscripts *CB2001* and *Data* indicate the theoretical values of yield stress or Lankford coefficient obtained with the anisotropic model, and the data from Table 1, respectively. Since for sheet materials tests involving out-of-plane stress components are very difficult or non-feasible, in the 3-D FE calculations, isotropic values are assigned to the anisotropy coefficients associated to the out-of-plane stress components, namely: *a*<sub>5</sub>=*a*<sub>6</sub>=*b*<sub>6</sub>=*b*<sub>7</sub>=*b*<sub>8</sub>=*b*<sub>9</sub>=*b*<sub>11</sub>=1 (e.g., [16]). Table 2 compiles the numerical values of the anisotropic parameters obtained from the identification procedure.

<table><caption>Table 2. Parameters CB2001 for all the materials listed Table 3.</caption><thead><tr><th></th><th>a1</th><th>a2</th><th>a3</th><th>a4</th><th>c</th><th>b1</th><th>b2</th><th>b3</th><th>b4</th><th>b5</th><th>b10</th><th>c</th></tr></thead><tbody><tr><td>Material 1</td><td>0.87</td><td>0.87</td><td>0.36</td><td>0.74</td><td>0.03</td><td>1.06</td><td>-1.95</td><td>-10.0</td><td>-9.11</td><td>-4.01</td><td>-0.25</td><td>0.03</td></tr><tr><td>Material 2</td><td>0.58</td><td>0.46</td><td>0.61</td><td>0.48</td><td>0.38</td><td>0.20</td><td>2.13</td><td>-0.02</td><td>0.24</td><td>0.74</td><td>0.63</td><td>0.38</td></tr></tbody></table>

<figure><img src="image_1037.png" /><figcaption>Fig. 1. Prediction of (a) Lankford coefficients, (b) normalized flow stresses and (c) projection in the biaxial ($\sigma_{xx} - \sigma_{yy}$) plane (with $\sigma_{xy} = 0$) of the yield locus according to CB2001 for the two materials.</figcaption></figure>

Fig. 1a and Fig. 1b show the predicted Lankford coefficients and normalized flow stresses according to CB2001 for both materials in Table 1. Moreover, in Fig. 1c, the projections in the biaxial plane ($\sigma_{xx} - \sigma_{yy}$) plane (with $\sigma_{xy} = 0$) of the yield loci are given. The materials hardening
behavior is considered isotropic and governed by the plastic strain. The power law in Eq. (10) is used to describe the flow curve for the RD sample.

$$ (\bar{\varepsilon}_p) = \varepsilon_0 \cdot [\varepsilon_{p0} + \bar{\varepsilon}_p]^n \quad (10) $$

The hardening parameters are identified based on the quasi-static static response of Ti4Al4V [8]. Where $\bar{\varepsilon}_p$ is the plastic strain, $\varepsilon_0 = 1472.15$ MPa, $\varepsilon_{p0} = 0.039$, $n = 0.127$ are hardening parameters. Standard values for the elastic properties are used, namely Young's modulus $E=117$ GPa, Poisson's ratio $\nu = 0.3$, and density $\rho = 4500$ kg/m³.

## Finite element model

Finite element simulations of tensile experiments are performed using the commercial finite element code Abaqus/Explicit. The CB2001 yield criterion previously presented is implemented as a user-defined material subroutine (VUMAT). In the simulation, the gauge length of the sample and part of the clamping region are modeled (see Fig. 2). The boundary conditions are imposed to the sample through control points. The control points are coupled to the right-hand and left-hand sides of the sample. A kinematic coupling type is used to couple all degrees of freedom of the control points, i.e., displacements (Ux,Uy,Uz), and rotations (URx, URy, URz), to the two sample end surfaces (see also Fig. 2). All degrees of freedom of the left-hand control point are fixed. A horizontal linear displacement of 3mm is imposed to the control point at the right-hand side of the sample. The tensile specimen is meshed with hexahedral elements with reduced integration (C3D8R in Abaqus) of size 0.1mm x 0.1mm x 0.1mm, resulting in a total of 60852 elements. FE simulations are performed for several in-plane orientations, namely $\theta = 0°, 15°, 30°, 45°, 60°, 75°$ and $90°$ to RD.

<figure><img src="image_1033.png" /><figcaption>Fig. 2. Schematic representation FE simulation uniaxial tensile experiment.</figcaption></figure>

## Results and discussion

In this section, FE results of uniaxial tensile simulations for all anisotropic materials listed in Table 1 are presented and discussed. First, the effect of both Lankford coefficient and normalized flow stress on the sign of the necking bands, is investigated. Next, the analytic predictions of the inclination of the necking band calculated using Eq. (2) are compared to the FE results. The angle $\beta$ between the loading direction and the localized necking band is considered to be positive for a counter-clockwise band orientation and negative for a clockwise band inclination (see also Fig. 3).
behavior is considered isotropic and governed by the plastic strain. The power law in Eq. (10) is used to describe the flow curve for the RD sample.

$$ (\bar{\varepsilon}_p) = \varepsilon_0 \cdot [\varepsilon_{p0} + \bar{\varepsilon}_p]^n \quad (10) $$

The hardening parameters are identified based on the quasi-static static response of Ti4Al4V [8]. Where $\bar{\varepsilon}_p$ is the plastic strain, $\varepsilon_0 = 1472.15$ MPa, $\varepsilon_{p0} = 0.039$, $n = 0.127$ are hardening parameters. Standard values for the elastic properties are used, namely Young's modulus $E=117$ GPa, Poisson's ratio $\nu = 0.3$, and density $\rho = 4500$ kg/m$^3$.

## Finite element model

Finite element simulations of tensile experiments are performed using the commercial finite element code Abaqus/Explicit. The CB2001 yield criterion previously presented is implemented as a user-defined material subroutine (VUMAT). In the simulation, the gauge length of the sample and part of the clamping region are modeled (see Fig. 2). The boundary conditions are imposed to the sample through control points. The control points are coupled to the right-hand and left-hand sides of the sample. A kinematic coupling type is used to couple all degrees of freedom of the control points, i.e., displacements (Ux,Uy,Uz), and rotations (URx, URy, URz), to the two sample end surfaces (see also Fig. 2). All degrees of freedom of the left-hand control point are fixed. A horizontal linear displacement of 3mm is imposed to the control point at the right-hand side of the sample. The tensile specimen is meshed with hexahedral elements with reduced integration (C3D8R in Abaqus) of size 0.1mm x 0.1mm x 0.1mm, resulting in a total of 60852 elements. FE simulations are performed for several in-plane orientations, namely $\theta = 0^\circ, 15^\circ, 30^\circ, 45^\circ, 60^\circ, 75^\circ$ and $90^\circ$ to RD.

<figure><img src="image_1038.png" /><figcaption>Fig. 2. Schematic representation FE simulation uniaxial tensile experiment.</figcaption></figure>

## Results and discussion

In this section, FE results of uniaxial tensile simulations for all anisotropic materials listed in Table 1 are presented and discussed. First, the effect of both Lankford coefficient and normalized flow stress on the sign of the necking bands, is investigated. Next, the analytic predictions of the inclination of the necking band calculated using Eq. (2) are compared to the FE results. The angle $\beta$ between the loading direction and the localized necking band is considered to be positive for a counter-clockwise band orientation and negative for a clockwise band inclination (see also Fig. 3).
FE contour plots of the equivalent plastic strain associated with the CB2001 yield criterion for all considered materials are shown in Fig. 3. The necking bands along seven loading directions, namely $\theta = 0°, 15°, 30°, 45°, 60°, 75°$ and $90°$ to the rolling direction (RD), are reported along with the predictions of the r-values (blue solid line) and normalized flow stresses (red dashed line), according to CB2001. FE contour plots are given for a plastic strain at the sample surface of 1. Moreover, in Fig. 3, two different regions are distinguished: the blue region indicates the loading directions for which a negative clockwise band angle ($\beta < 0$) can be detected, while the red region includes all directions for which a positive counter-clockwise band angle ($\beta > 0$) is observed. In all contour plots, two bands of equal inclination to the loading axis develop for samples loaded along the principal axis of orthotropy RD and TD. Moreover, there exists an additional off-axis specimen orientation $\theta^*$ for which bands of equal inclination develop.

<figure><img src="image_1034.png" /><figcaption>Fig. 3. Contour plot at 1% of equivalent plastic strain of uniaxial experiments for seven in-plane orientations, from RD to TD, for (a) Material 1 and (b) Material 2.</figcaption></figure>

In Fig. 3a, Fig. 3b the extrema of the uniaxial flow stresses correspond to $\theta^* \sim 60°$ and $\theta^* \sim 45°$ for the Material 1 and 2, respectively. In particular, FE simulations for the Material 1 result in necking bands being negative ($\beta < 0$) for $\theta < \theta^*$ and positive ($\beta > 0$) for $\theta > \theta^*$. On the contrary, for the Material 2, localized bands are found to be positive ($\beta > 0$) for $\theta < \theta^*$ and negative ($\beta < 0$) for $\theta > \theta^*$. It is worth to highlight that both considered materials present very similar anisotropy in r-values but opposite trends in yield stresses. Indeed, for the Material 1 the predicted $\sigma_\theta / \sigma_0$ vs. $\theta$ is a monotonically decreasing function in the range $0° < \theta < \theta^*$, while it is a monotonically increasing function in the range $\theta^* < \theta < 90°$. On the other hand, for the Material 2, the opposite is true. Particularly, negative necking bands ($\beta < 0$) are found for $\sigma_\theta / \sigma_0$ vs. $\theta$ being a monotonically decreasing function, while positive necking bands ($\beta > 0$) are found for $\sigma_\theta / \sigma_0$ vs. $\theta$ being a monotonically increasing function. As such, FE results again suggest that the plastic flow anisotropy does not influence the sign of the necking bands. It is possible to conclude that the anisotropy in flow stress, i.e., $\sigma_\theta / \sigma_0$ vs. $\theta$, is the only responsible for the sign of the necking bands.
FE contour plots of the equivalent plastic strain associated with the CB2001 yield criterion for all considered materials are shown in Fig. 3. The necking bands along seven loading directions, namely $\theta = 0°, 15°, 30°, 45°, 60°, 75°$ and $90°$ to the rolling direction (RD), are reported along with the predictions of the r-values (blue solid line) and normalized flow stresses (red dashed line), according to CB2001. FE contour plots are given for a plastic strain at the sample surface of 1. Moreover, in Fig. 3, two different regions are distinguished: the blue region indicates the loading directions for which a negative clockwise band angle ($\beta < 0$) can be detected, while the red region includes all directions for which a positive counter-clockwise band angle ($\beta > 0$) is observed. In all contour plots, two bands of equal inclination to the loading axis develop for samples loaded along the principal axis of orthotropy RD and TD. Moreover, there exists an additional off-axis specimen orientation $\theta^*$ for which bands of equal inclination develop.

<figure><img src="image_1039.png" /><figcaption>Fig. 3. Contour plot at 1% of equivalent plastic strain of uniaxial experiments for seven in-plane orientations, from RD to TD, for (a) Material 1 and (b) Material 2.</figcaption></figure>

In Fig. 3a, Fig. 3b the extrema of the uniaxial flow stresses correspond to $\theta^* \sim 60°$ and $\theta^* \sim 45°$ for the Material 1 and 2, respectively. In particular, FE simulations for the Material 1 result in necking bands being negative ($\beta < 0$) for $\theta < \theta^*$ and positive ($\beta > 0$) for $\theta > \theta^*$. On the contrary, for the Material 2, localized bands are found to be positive ($\beta > 0$) for $\theta < \theta^*$ and negative ($\beta < 0$) for $\theta > \theta^*$. It is worth to highlight that both considered materials present very similar anisotropy in r-values but opposite trends in yield stresses. Indeed, for the Material 1 the predicted $\sigma_\theta / \sigma_0$ vs. $\theta$ is a monotonically decreasing function in the range $0° < \theta < \theta^*$, while it is a monotonically increasing function in the range $\theta^* < \theta < 90°$. On the other hand, for the Material 2, the opposite is true. Particularly, negative necking bands ($\beta < 0$) are found for $\sigma_\theta / \sigma_0$ vs. $\theta$ being a monotonically decreasing function, while positive necking bands ($\beta > 0$) are found for $\sigma_\theta / \sigma_0$ vs. $\theta$ being a monotonically increasing function. As such, FE results again suggest that the plastic flow anisotropy does not influence the sign of the necking bands. It is possible to conclude that the anisotropy in flow stress, i.e., $\sigma_\theta / \sigma_0$ vs. $\theta$, is the only responsible for the sign of the necking bands.
It is worth to mention that the analytical prediction of the necking bands angle obtained solving Eq. (2) for a generic loading direction θ results in two solutions β₁ and β₂ being the inclinations of the positive and negative necking bands. Therefore, if the trend of the yield stress anisotropy is known, also the sign of the band for a generic orientation θ in the plane of the sheet is known beforehand. Fig. 4a and Fig. 4b, show the analytical and FE predictions of the band inclinations according to CB2001 for the Material 1 and Material 2, respectively. In FE simulations, the band inclinations are determined at a plastic strain of 0.2. A vertical dashed line indicates the loading direction θ* for which two equally inclined bands are expected. The trends of the necking angle |β| vs. θ obtained by FE and the analytical prediction are in good agreement for both the considered materials. In the case of uniaxial loadings along RD and TD (i.e., θ = 0° or θ = 90°) B = 0 (see Eq. (4)), and two bands of equal inclinations are obtained. Moreover, for the Material 1, since the flow curve admits a unique minimum in the range 0° < θ < 90° (see also Fig. 4a), then it exists another orientation θ*=57° for which two equally inclined bands appear. For the Material 1 the following is true:

* for specimens with orientations θ < θ* = 57°, the inclination of the band which makes an acute counter-clockwise angle with the loading direction is greater than the inclination of the band which makes an acute clockwise angle with the loading direction, i.e., β₁ > |β₂|.

* for specimen orientations θ > θ* = 57° the reverse occurs, i.e., β₁ < |β₂|.

On the other hand, in the case of the Material 2, since the flow curve admits a maximum at ~45°, the predominant band making an acute angle for off-axis specimen of orientation is greater than the inclination of the band which makes an acute counter-clockwise angle.

<figure><img src="image_1035.png" /><figcaption>Fig. 4. Localization band inclination as a function of the specimen orientation θ with respect to RD for the a) Material 1 and b) Material 2, using Eq. (2) with the CB2001 effective stress given by Eq. (8). The vertical dashed line corresponds to the orientation θ*. FE results are represented by hollow circles.</figcaption></figure>

## Summary and conclusions

In present work, FE simulations were carried out to investigate the role of the material anisotropy on the inception of necking bands in uniaxial tensile dog bone samples for two virtual materials having mild anisotropy in yield stresses and pronounced anisotropy in r-values. The two considered materials have same anisotropy in r-values but different trends in yield stress anisotropy. For the Material 1, the function σθ / σ₀ vs. θ presents a minimum for θ = 45°, while for the same angle θ, the function σθ / σ₀ vs. θ associated to the Material 2 presents a maximum. The CB2001 yield criterion was used to model the anisotropic response of the two materials, and
It is worth to mention that the analytical prediction of the necking bands angle obtained solving Eq. (2) for a generic loading direction θ results in two solutions β₁ and β₂ being the inclinations of the positive and negative necking bands. Therefore, if the trend of the yield stress anisotropy is known, also the sign of the band for a generic orientation θ in the plane of the sheet is known beforehand. Fig. 4a and Fig. 4b, show the analytical and FE predictions of the band inclinations according to CB2001 for the Material 1 and Material 2, respectively. In FE simulations, the band inclinations are determined at a plastic strain of 0.2. A vertical dashed line indicates the loading direction θ* for which two equally inclined bands are expected. The trends of the necking angle |β| vs. θ obtained by FE and the analytical prediction are in good agreement for both the considered materials. In the case of uniaxial loadings along RD and TD (i.e., θ = 0° or θ = 90°) B = 0 (see Eq. (4)), and two bands of equal inclinations are obtained. Moreover, for the Material 1, since the flow curve admits a unique minimum in the range 0° < θ < 90° (see also Fig. 4a), then it exists another orientation θ*=57° for which two equally inclined bands appear. For the Material 1 the following is true:

* for specimens with orientations θ < θ* = 57°, the inclination of the band which makes an acute counter-clockwise angle with the loading direction is greater than the inclination of the band which makes an acute clockwise angle with the loading direction, i.e., β₁ > |β₂|.

* for specimen orientations θ > θ* = 57° the reverse occurs, i.e., β₁ < |β₂|.

On the other hand, in the case of the Material 2, since the flow curve admits a maximum at ~45°, the predominant band making an acute angle for off-axis specimen of orientation is greater than the inclination of the band which makes an acute counter-clockwise angle.

<figure><img src="image_1040.png" /><figcaption>Fig. 4. Localization band inclination as a function of the specimen orientation θ with respect to RD for the a) Material 1 and b) Material 2, using Eq. (2) with the CB2001 effective stress given by Eq. (8). The vertical dashed line corresponds to the orientation θ*. FE results are represented by hollow circles.</figcaption></figure>

## Summary and conclusions

In present work, FE simulations were carried out to investigate the role of the material anisotropy on the inception of necking bands in uniaxial tensile dog bone samples for two virtual materials having mild anisotropy in yield stresses and pronounced anisotropy in r-values. The two considered materials have same anisotropy in r-values but different trends in yield stress anisotropy. For the Material 1, the function σθ / σ₀ vs. θ presents a minimum for θ = 45°, while for the same angle θ, the function σθ / σ₀ vs. θ associated to the Material 2 presents a maximum. The CB2001 yield criterion was used to model the anisotropic response of the two materials, and
FE simulations were performed in Abaqus/Explicit. Simulations of tensile experiments were carried out along seven loading directions, namely $θ = 0°, 15°, 30°, 45°, 60°, 75°$ and $90°$ to the rolling direction (RD). Necking bands of equal inclination to the loading axis developed for samples loaded along the principal axes of orthotropy RD and TD. Moreover, it was found that an additional off-axis orientation $θ*$ in the plane of the sheet exists for which bands of equal inclination developed. However, this angle changes with the material anisotropy. This result is in agreement with the work of Cazacu and Rodríguez-Martínez [6], which proved that only for specimens with orientations $θ*$ corresponding to the extrema, i.e., minimum or maximum, of the uniaxial tensile flow stresses in the plane of the sheet, the two localization bands are equally inclined. Concerning the other off-axis loading directions, numerical results put into evidence the existence of a relation between the evolution of the anisotropy in flow stress, i.e., $σ_θ / σ_0$ vs. $θ$, and the sign of the necking bands. In particular, negative necking bands ($β < 0$) were found for $σ_θ / σ_0$ vs. $θ$ being a monotonically decreasing function, while positive necking bands ($β > 0$) are found for $σ_θ / σ_0$ vs. $θ$ being a monotonically increasing function. Next, the theoretical analytical model proposed by Cazacu and Rodríguez-Martínez [1] was used to calculate the inclination of the necking bands for the two considered materials. The good agreement between FE results and analytical calculations put into evidence the good predictive capability of the analytical calculations. As a general remark, FE simulations put into evidence the key role of material anisotropy, and the importance of its correct modelling, to effectively predict the localization phenomenon arising during uniaxial tensile experiments.

## References

[1] Cazacu O, Rodríguez-Martínez JA. Effects of plastic anisotropy on localization in orthotropic materials: new explicit expressions for the orientation of localization bands in flat specimens subjected to uniaxial tension. Journal of the Mechanics and Physics of Solids, 2019;126:272–84. https://doi.org/10.1016/j.jmps.2019.03.002

[2] Nadai, A. Plasticity 1931;Vol.1:Mc Graw-Hill, New York.

[3] Bell JF. The Physics of Large Deformation of Crystalline Solids. vol. 14. Berlin, Heidelberg: Springer; 1968. https://doi.org/10.1007/978-3-642-88440-5

[4] Neilsen MK, Schreyer HL. Bifurcations in elastic-plastic materials. International Journal of Solids and Structures 1993;30:521–44. https://doi.org/10.1016/0020-7683(93)90185-A

[5] Korber, F., Hoff, H. Ueber die Festigkeitseigenschaften und den reißwinkel kaltgewalzter metalle. Mitteilungen Aus Dem Kaiser-Wilhelm-Institut Fuer Eisenforschung Zu Duesseldorf 10 1928;175–187.

[6] Hill R. A theory of the yielding and plastic flow of anisotropic metals. Proc R Soc Lond A 1948;193:281–97. https://doi.org/10.1098/rspa.1948.0045

[7] N’souglo KE, Rodríguez-Martínez JA, Vaz-Romero A, Cazacu O. The combined effect of plastic orthotropy and tension-compression asymmetry on the development of necking instabilities in flat tensile specimens subjected to dynamic loading. International Journal of Solids and Structures 2019;159:272–88. https://doi.org/10.1016/j.ijsolstr.2018.10.006

[8] Corallo L, Revil-Baudard B, Cazacu O, Verleysen P. Role of anisotropy on strain localization at high-strain rates in Ti6Al4V. European Journal of Mechanics - A/Solids 2023;97:104856. https://doi.org/10.1016/j.euromechsol.2022.104856

[9] Tvergaard V. Necking in tensile bars with rectangular cross-section. Computer Methods in Applied Mechanics and Engineering 1993;103:273–90. https://doi.org/10.1016/0045-7825(93)90049-4
FE simulations were performed in Abaqus/Explicit. Simulations of tensile experiments were carried out along seven loading directions, namely $θ = 0°, 15°, 30°, 45°, 60°, 75°$ and $90°$ to the rolling direction (RD). Necking bands of equal inclination to the loading axis developed for samples loaded along the principal axes of orthotropy RD and TD. Moreover, it was found that an additional off-axis orientation $θ*$ in the plane of the sheet exists for which bands of equal inclination developed. However, this angle changes with the material anisotropy. This result is in agreement with the work of Cazacu and Rodríguez-Martínez [6], which proved that only for specimens with orientations $θ*$ corresponding to the extrema, i.e., minimum or maximum, of the uniaxial tensile flow stresses in the plane of the sheet, the two localization bands are equally inclined. Concerning the other off-axis loading directions, numerical results put into evidence the existence of a relation between the evolution of the anisotropy in flow stress, i.e., $σ_θ / σ_0$ vs. $θ$, and the sign of the necking bands. In particular, negative necking bands ($β < 0$) were found for $σ_θ / σ_0$ vs. $θ$ being a monotonically decreasing function, while positive necking bands ($β > 0$) are found for $σ_θ / σ_0$ vs. $θ$ being a monotonically increasing function. Next, the theoretical analytical model proposed by Cazacu and Rodríguez-Martínez [1] was used to calculate the inclination of the necking bands for the two considered materials. The good agreement between FE results and analytical calculations put into evidence the good predictive capability of the analytical calculations. As a general remark, FE simulations put into evidence the key role of material anisotropy, and the importance of its correct modelling, to effectively predict the localization phenomenon arising during uniaxial tensile experiments.

## References

[1] Cazacu O, Rodríguez-Martínez JA. Effects of plastic anisotropy on localization in orthotropic materials: new explicit expressions for the orientation of localization bands in flat specimens subjected to uniaxial tension. Journal of the Mechanics and Physics of Solids, 2019;126:272-84. https://doi.org/10.1016/j.jmps.2019.03.002

[2] Nadai, A. Plasticity 1931;Vol.1:Mc Graw-Hill, New York.

[3] Bell JF. The Physics of Large Deformation of Crystalline Solids. vol. 14. Berlin, Heidelberg: Springer; 1968. https://doi.org/10.1007/978-3-642-88440-5

[4] Neilsen MK, Schreyer HL. Bifurcations in elastic-plastic materials. International Journal of Solids and Structures 1993;30:521-44. https://doi.org/10.1016/0020-7683(93)90185-A

[5] Korber, F., Hoff, H. Ueber die Festigkeitseigenschaften und den reißwinkel kaltgewalzter metalle. Mitteilungen Aus Dem Kaiser-Wilhelm-Institut Fuer Eisenforschung Zu Duesseldorf 10 1928;175-187.

[6] Hill R. A theory of the yielding and plastic flow of anisotropic metals. Proc R Soc Lond A 1948;193:281-97. https://doi.org/10.1098/rspa.1948.0045

[7] N'souglo KE, Rodríguez-Martínez JA, Vaz-Romero A, Cazacu O. The combined effect of plastic orthotropy and tension-compression asymmetry on the development of necking instabilities in flat tensile specimens subjected to dynamic loading. International Journal of Solids and Structures 2019;159:272-88. https://doi.org/10.1016/j.ijsolstr.2018.10.006

[8] Corallo L, Revil-Baudard B, Cazacu O, Verleysen P. Role of anisotropy on strain localization at high-strain rates in Ti6Al4V. European Journal of Mechanics - A/Solids 2023;97:104856. https://doi.org/10.1016/j.euromechsol.2022.104856

[9] Tvergaard V. Necking in tensile bars with rectangular cross-section. Computer Methods in Applied Mechanics and Engineering 1993;103:273-90. https://doi.org/10.1016/0045-7825(93)90049-4
[10] Belytschko T, Chiang H-Y, Plaskacz E. High resolution two-dimensional shear band computations: imperfections and mesh dependence. Computer Methods in Applied Mechanics and Engineering 1994;119:1–15. https://doi.org/10.1016/0045-7825(94)00073-5

[11] Petryk H. Plastic instability: Criteria and computational approaches. ARCO 1997;4:111–51. https://doi.org/10.1007/BF03020127

[12] Needleman A, Rice J. Limits to Ductility Set by Plastic Flow Localization 1978. https://doi.org/10.1007/978-1-4613-2880-3_10

[13] Labergere C, Guelorget B, Francois M. Strain rate distribution and localization band width evolution during tensile test. International Journal of Solids and Structures 2014;51:3944–61. https://doi.org/10.1016/j.ijsolstr.2014.07.017

[14] Verleysen P, Peirs J. Quasi-static and high strain rate fracture behaviour of Ti6Al4V. International Journal of Impact Engineering 2017;108:370–88. https://doi.org/10.1016/j.ijimpeng.2017.03.001

[15] Cazacu O, Barlat F. Generalization of Drucker's Yield Criterion to Orthotropy. Mathematics and Mechanics of Solids 2001;6:613–30. https://doi.org/10.1177/108128650100600603

[16] Cazacu O, Revil-Baudard B. Plasticity of Metallic Materials: Modelling and Applications to Metal Forming. Elsevier. ISBN: 9780128179840.; 2021.

[17] Corallo L, Verleysen P. On the Anisotropic Mechanical Response of Ti6Al4V Sheet at High Strain Rates. KEM 2022;926:970–9. https://doi.org/10.4028/p-xq42f4

[18] Gilles G, Hammami W, Libertiaux V, Cazacu O, Yoon JH, Kuwabara T, et al. Experimental characterization and elasto-plastic modeling of the quasi-static mechanical response of TA-6V at room temperature. International Journal of Solids and Structures 2011;48:1277–89. https://doi.org/10.1016/j.ijsolstr.2011.01.011

[19] Cazacu O, Revil-Baudard B, Chandola N. Plasticity-Damage Couplings: From Single Crystal to Polycrystalline Materials. vol. 253. Cham: Springer International Publishing; 2019. https://doi.org/10.1007/978-3-319-92922-4
[10] Belytschko T, Chiang H-Y, Plaskacz E. High resolution two-dimensional shear band computations: imperfections and mesh dependence. Computer Methods in Applied Mechanics and Engineering 1994;119:1–15. https://doi.org/10.1016/0045-7825(94)00073-5

[11] Petryk H. Plastic instability: Criteria and computational approaches. ARCO 1997;4:111–51. https://doi.org/10.1007/BF03020127

[12] Needleman A, Rice J. Limits to Ductility Set by Plastic Flow Localization 1978. https://doi.org/10.1007/978-1-4613-2880-3_10

[13] Labergere C, Guelorget B, Francois M. Strain rate distribution and localization band width evolution during tensile test. International Journal of Solids and Structures 2014;51:3944–61. https://doi.org/10.1016/j.ijsolstr.2014.07.017

[14] Verleysen P, Peirs J. Quasi-static and high strain rate fracture behaviour of Ti6Al4V. International Journal of Impact Engineering 2017;108:370–88. https://doi.org/10.1016/j.ijimpeng.2017.03.001

[15] Cazacu O, Barlat F. Generalization of Drucker's Yield Criterion to Orthotropy. Mathematics and Mechanics of Solids 2001;6:613–30. https://doi.org/10.1177/108128650100600603

[16] Cazacu O, Revil-Baudard B. Plasticity of Metallic Materials: Modelling and Applications to Metal Forming. Elsevier. ISBN: 9780128179840.; 2021.

[17] Corallo L, Verleysen P. On the Anisotropic Mechanical Response of Ti6Al4V Sheet at High Strain Rates. KEM 2022;926:970–9. https://doi.org/10.4028/p-xq42f4

[18] Gilles G, Hammami W, Libertiaux V, Cazacu O, Yoon JH, Kuwabara T, et al. Experimental characterization and elasto-plastic modeling of the quasi-static mechanical response of TA-6V at room temperature. International Journal of Solids and Structures 2011;48:1277–89. https://doi.org/10.1016/j.ijsolstr.2011.01.011

[19] Cazacu O, Revil-Baudard B, Chandola N. Plasticity-Damage Couplings: From Single Crystal to Polycrystalline Materials. vol. 253. Cham: Springer International Publishing; 2019. https://doi.org/10.1007/978-3-319-92922-4
# From strain to stress using full-field data:
## Computationally efficient stress reconstruction

HALILOVIČ Miroslav<fnref n="1" />,ᵃ, STARMAN Bojan<fnref n="1" />,<fnref n="2" />,ᵇ and COPPIETERS Sam<fnref n="2" />,ᶜ*

**Keywords:** Stress Reconstruction, Plasticity, Digital Image Correlation, Virtual Field Method

**Abstract.** Conventional stress reconstruction based on full-field strain measurements presents a major computational burden, especially when using standard implicit stress integration methods. This presents a notable challenge for inverse identification methods used to characterize the plasticity of metallic materials, particularly those reliant on stress reconstruction, such as the nonlinear sensitivity-based Virtual Fields Method (VFM). To reduce the computational effort, the full-field strain data are usually spatially and temporally down-sampled. However, for metals subject to nonlinear strain paths, this practice can lead to errors in the resulting stress states and compromise the accuracy of the nonlinear VFM. In this work, we introduce a highly efficient explicit stress reconstruction algorithm to reduce the computational challenges of repeated stress reconstruction which can be utilized in inverse identification methods such as nonlinear VFM.

## Introduction

The development of optical technologies has opened up new possibilities to characterize the behavior of materials in solid mechanics. The use of 3D scanners and cameras has enabled precise measurements of body displacements and deformations under load. The ability to accurately track and quantify displacements and strains allows engineers to assess material behavior and make informed design decisions. Digital image correlation (DIC) is one of the most commonly used techniques for calculating strain. In this method, sequentially captured images are processed to extract valuable information that facilitates the study of local behavior and global responses in structures. In addition, DIC can be used in real time and offers the possibility to monitor structural responses under dynamic loading conditions. This feature makes DIC a powerful tool for various applications, including material testing, structural analyses and product development.

Despite significant advances in metrology that enable precise measurements of the full evolution of the strain field during loading, the challenge of efficiently reconstructing the corresponding evolution of the stress field in real time remains. In the case of linear elastic materials, the reconstruction of the stress field from the strain field can be performed seamlessly. However, in the case of nonlinear material behavior, such as plasticity, the mechanical response depends on the loading path and is described by algebraic differential equations. The effective solution of these equations to obtain the stress state for a given time frame remains a challenge.

## State of the art

The effectiveness of stress reconstruction algorithms is of crucial importance in today's applications. The need arises, for example, in materials testing, where inverse identification methods are used to determine material behavior from measured heterogeneous strain data. An effective stress integration algorithm is essential in this context to ensure accurate and timely results.

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

<footnote n="1">¹University of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, 1000 Ljubljana, Slovenia</footnote>

<footnote n="2">²KU Leuven, Department of Materials Engineering, Ghent Campus, 9000 Ghent, Belgium</footnote>

<footnote>ᵃmiroslav.halilovic@fs.uni-lj.si, ᵇbojan.starman@fs.uni-lj.si, ᶜsam.coppieters@kuleuven.be</footnote>
# From strain to stress using full-field data:
## Computationally efficient stress reconstruction

HALILOVIČ Miroslav<fnref n="1" />,ª, STARMAN Bojan<fnref n="1" />,<fnref n="2" />,ᵇ and COPPIETERS Sam<fnref n="2" />,ᶜ*

**Keywords:** Stress Reconstruction, Plasticity, Digital Image Correlation, Virtual Field Method

**Abstract.** Conventional stress reconstruction based on full-field strain measurements presents a major computational burden, especially when using standard implicit stress integration methods. This presents a notable challenge for inverse identification methods used to characterize the plasticity of metallic materials, particularly those reliant on stress reconstruction, such as the nonlinear sensitivity-based Virtual Fields Method (VFM). To reduce the computational effort, the full-field strain data are usually spatially and temporally down-sampled. However, for metals subject to nonlinear strain paths, this practice can lead to errors in the resulting stress states and compromise the accuracy of the nonlinear VFM. In this work, we introduce a highly efficient explicit stress reconstruction algorithm to reduce the computational challenges of repeated stress reconstruction which can be utilized in inverse identification methods such as nonlinear VFM.

## Introduction

The development of optical technologies has opened up new possibilities to characterize the behavior of materials in solid mechanics. The use of 3D scanners and cameras has enabled precise measurements of body displacements and deformations under load. The ability to accurately track and quantify displacements and strains allows engineers to assess material behavior and make informed design decisions. Digital image correlation (DIC) is one of the most commonly used techniques for calculating strain. In this method, sequentially captured images are processed to extract valuable information that facilitates the study of local behavior and global responses in structures. In addition, DIC can be used in real time and offers the possibility to monitor structural responses under dynamic loading conditions. This feature makes DIC a powerful tool for various applications, including material testing, structural analyses and product development.

Despite significant advances in metrology that enable precise measurements of the full evolution of the strain field during loading, the challenge of efficiently reconstructing the corresponding evolution of the stress field in real time remains. In the case of linear elastic materials, the reconstruction of the stress field from the strain field can be performed seamlessly. However, in the case of nonlinear material behavior, such as plasticity, the mechanical response depends on the loading path and is described by algebraic differential equations. The effective solution of these equations to obtain the stress state for a given time frame remains a challenge.

## State of the art

The effectiveness of stress reconstruction algorithms is of crucial importance in today's applications. The need arises, for example, in materials testing, where inverse identification methods are used to determine material behavior from measured heterogeneous strain data. An effective stress integration algorithm is essential in this context to ensure accurate and timely results.

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

<footnote n="1">¹University of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, 1000 Ljubljana, Slovenia</footnote>

<footnote n="2">²KU Leuven, Department of Materials Engineering, Ghent Campus, 9000 Ghent, Belgium</footnote>

<footnote>ªmiroslav.halilovic@fs.uni-lj.si, ᵇbojan.starman@fs.uni-lj.si, ᶜsam.coppieters@kuleuven.be</footnote>
In other words, the development of full-field measurement techniques which provide a wealth of information has played a crucial role in the development of novel identification strategies. These strategies aim to reduce the calibration effort compared to conventional material testing. Information-rich experiments take advantage of the heterogeneity of deformation to determine multiple material model parameters simultaneously. Most of these experiments use DIC to capture heterogeneous strain fields and subject them to an inverse identification process. This approach has unlocked the potential for inverse material characterization and led to the development of new identification methods, such as the Virtual Field Method (VFM) [1], the Equilibrium Gap Method (EGM) [2] and the Constitutive Equation Gap Method (CEGM) [3], Reciprocity Gap Method (RGM) [4], Finite Element Model Updating (FEMU) [5], Integrated Mechanical Image Correlation (I-MIC) [6], Integrated Digital Image Correlation (I-DIC) [7] and various combined theoretical/experimental approaches such as [8], [9].

Initially, these methods were applied to elasticity problems, but today, the most commonly used methods for plasticity problems are FEMU and VFM. FEMU is a user-friendly method that uses an optimization algorithm to minimize the difference between measured and computed Finite Element Analysis (FEA) results. It has been used to identify parameters for various plasticity models, such as Hill48 [10] and YLD2000-2d [11]. VFM, on the other hand, is based on the principle of virtual work and force equilibrium. It was initially used to determine elastic properties [12], but later, the method was also successfully applied to plasticity problems [13], [14].

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

As an example, Lattanzi et al. [15] reported that the computational time for the VFM ranges from 11 to 83 minutes on an Intel Xeon CPU E5-1650 v2 @ 3.50 GHz and 128 GB of RAM. However, it's important to note that the computational time can vary greatly depending on various factors such as the size of the model, the complexity of material behavior, and the computational resources available. Despite these variations, the superiority of VFM over FEMU has been demonstrated by Rahmani et al. [16] and Zhang et al. [17], who reported that VFM is 125 times faster than conventional FEMU for their specific application. In a study by Marek et al. [18], the identification time was found to range from 4 to 278 hours, depending on the material model and identification algorithm used. They concluded that the stress reconstruction algorithm was the main limiting factor in efficiency and suggested replacing the implicit stress reconstruction algorithm with a more efficient explicit one, as proposed by Rossi et al. [19]. A similar conclusion was found in [20], where the computational time for VFM was reported to range from 25 to 35 minutes.

As presented, both methods, FEMU and VFM, utilize the inhomogeneous strain field processed by DIC, but they face efficiency issues when it comes to reconstructing the stress state from a given strain field. This limitation hampers the efficiency of material identification when a large number of strain increments are used to describe the loading path.

In this paper, our primary objective is to improve the efficiency of conventional stress reconstruction algorithms commonly employed in inverse material identification methods, such as the Virtual Fields Method (VFM), while addressing and overcoming the limitations inherent in traditional integration techniques. Dense temporal sampling drastically increases the
In other words, the development of full-field measurement techniques which provide a wealth of information has played a crucial role in the development of novel identification strategies. These strategies aim to reduce the calibration effort compared to conventional material testing. Information-rich experiments take advantage of the heterogeneity of deformation to determine multiple material model parameters simultaneously. Most of these experiments use DIC to capture heterogeneous strain fields and subject them to an inverse identification process. This approach has unlocked the potential for inverse material characterization and led to the development of new identification methods, such as the Virtual Field Method (VFM) [1], the Equilibrium Gap Method (EGM) [2] and the Constitutive Equation Gap Method (CEGM) [3], Reciprocity Gap Method (RGM) [4], Finite Element Model Updating (FEMU) [5], Integrated Mechanical Image Correlation (I-MIC) [6], Integrated Digital Image Correlation (I-DIC) [7] and various combined theoretical/experimental approaches such as [8], [9].

Initially, these methods were applied to elasticity problems, but today, the most commonly used methods for plasticity problems are FEMU and VFM. FEMU is a user-friendly method that uses an optimization algorithm to minimize the difference between measured and computed Finite Element Analysis (FEA) results. It has been used to identify parameters for various plasticity models, such as Hill48 [10] and YLD2000-2d [11]. VFM, on the other hand, is based on the principle of virtual work and force equilibrium. It was initially used to determine elastic properties [12], but later, the method was also successfully applied to plasticity problems [13], [14].

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

As an example, Lattanzi et al. [15] reported that the computational time for the VFM ranges from 11 to 83 minutes on an Intel Xeon CPU E5-1650 v2 @ 3.50 GHz and 128 GB of RAM. However, it's important to note that the computational time can vary greatly depending on various factors such as the size of the model, the complexity of material behavior, and the computational resources available. Despite these variations, the superiority of VFM over FEMU has been demonstrated by Rahmani et al. [16] and Zhang et al. [17], who reported that VFM is 125 times faster than conventional FEMU for their specific application. In a study by Marek et al. [18], the identification time was found to range from 4 to 278 hours, depending on the material model and identification algorithm used. They concluded that the stress reconstruction algorithm was the main limiting factor in efficiency and suggested replacing the implicit stress reconstruction algorithm with a more efficient explicit one, as proposed by Rossi et al. [19]. A similar conclusion was found in [20], where the computational time for VFM was reported to range from 25 to 35 minutes.

As presented, both methods, FEMU and VFM, utilize the inhomogeneous strain field processed by DIC, but they face efficiency issues when it comes to reconstructing the stress state from a given strain field. This limitation hampers the efficiency of material identification when a large number of strain increments are used to describe the loading path.

In this paper, our primary objective is to improve the efficiency of conventional stress reconstruction algorithms commonly employed in inverse material identification methods, such as the Virtual Fields Method (VFM), while addressing and overcoming the limitations inherent in traditional integration techniques. Dense temporal sampling drastically increases the
computational cost. Such a situation is forcing to reconsider the stress reconstruction algorithms: are the algorithms known for stress updating in FEM also optimal for stress reconstruction associated with full field strain measurements? Not necessarily; a study by Rossi et al. [21] found that there is an approximated method for stress reconstruction from measured strain data that does not require iteration and is capable of handling large increments.

The stress update algorithms used in FEM are typically based on unconditionally stable implicit integration methods, with the implicit backward Euler scheme being the most widely used approach. This method is widely used in commercial software and among researchers for various models (referenced in [22]). However, when using implicit integration for stress reconstruction in each time frame, the computation time depends directly on the number of time frames, creating a dilemma of how many time frames should be considered to achieve satisfactory accuracy. There is currently no clear answer to this question in the field of plasticity.

The main message we want to convey is that the dilemma of determining the number of time frames for stress path integration can be resolved by using explicit schemes. These explicit schemes solve the problem by using numerous, computationally effective increments, and if the increments are small due to the algorithm itself, this integration approach can cover all recorded time frames, eliminating the need to decide on the amount of time frames to consider.

In the plasticity domain, explicit schemes have a low reputation due to their conditional stability and permanent drift. However, in [23], [24], we presented an explicit integration scheme Next Increment Corrects Error (NICE), which is suitable for explicit FEM and does not produce permanent drift. The scheme has been further developed for higher order accuracy [25] and extended to be suitable for implicit FEM, integrating with the maximum possible stable subincrements and preserving accuracy within each increment [26]. Guaranteed stability, accuracy, and fast computation are crucial for efficient stress reconstruction from strain data.

In this paper, we present an adaptation of the NICE algorithm to be used with strain field data obtained through image processing. This adaptation is independent of the number of frames, allowing for full use of full-field data. The proposed algorithm differs from conventional stress reconstruction algorithms as it possesses the necessary properties for inverse identification.

## Methodology

In this paper we propose a methodology to meet the main objectives of examining the influence of strain temporal down-sampling on the stress reconstruction and avoiding the efficiency bottleneck of conventional stress reconstruction algorithms. We conducted a simple synthetic test, an open-hole test experiment, which is commonly used in heterogeneous specimen testing. The experiment is mimicked in ABAQUS/Standard, from which the displacement field is extracted in each increment for each nodal position in FE mesh. In the next step, the extracted displacement field is used to numerically deform a reference speckle pattern, which is used to generate synthetically deformed images. These numerically deformed images are afterwards passed into DIC software to process the corresponding strain fields. In the final step, the extracted strain data is used for stress computation using both conventional Implicit Radial Mapping Algorithm (IRMA) and NICE algorithms for different sets of strain increments. For instance, the whole strain path can be integrated at each image frame or using only a selected number of frames. The difference in stress error due to the strain path approximation is monitored. This scenario is commonly encountered in real DIC application where numerous images representing incremental strain evolution are available, but stress reconstruction algorithms like VFM do not consider all increments due to efficiency limitations.

Finally, we analyze the efficiency of the IRMA and NICE algorithms when dealing with the entire strain path divided into different numbers of strain increments used in the integration process. The NICE algorithm, being explicit in nature, is expected to offer higher computational efficiency.
computational cost. Such a situation is forcing to reconsider the stress reconstruction algorithms: are the algorithms known for stress updating in FEM also optimal for stress reconstruction associated with full field strain measurements? Not necessarily; a study by Rossi et al. [21] found that there is an approximated method for stress reconstruction from measured strain data that does not require iteration and is capable of handling large increments.

The stress update algorithms used in FEM are typically based on unconditionally stable implicit integration methods, with the implicit backward Euler scheme being the most widely used approach. This method is widely used in commercial software and among researchers for various models (referenced in [22]). However, when using implicit integration for stress reconstruction in each time frame, the computation time depends directly on the number of time frames, creating a dilemma of how many time frames should be considered to achieve satisfactory accuracy. There is currently no clear answer to this question in the field of plasticity.

The main message we want to convey is that the dilemma of determining the number of time frames for stress path integration can be resolved by using explicit schemes. These explicit schemes solve the problem by using numerous, computationally effective increments, and if the increments are small due to the algorithm itself, this integration approach can cover all recorded time frames, eliminating the need to decide on the amount of time frames to consider.

In the plasticity domain, explicit schemes have a low reputation due to their conditional stability and permanent drift. However, in [23], [24], we presented an explicit integration scheme Next Increment Corrects Error (NICE), which is suitable for explicit FEM and does not produce permanent drift. The scheme has been further developed for higher order accuracy [25] and extended to be suitable for implicit FEM, integrating with the maximum possible stable subincrements and preserving accuracy within each increment [26]. Guaranteed stability, accuracy, and fast computation are crucial for efficient stress reconstruction from strain data.

In this paper, we present an adaptation of the NICE algorithm to be used with strain field data obtained through image processing. This adaptation is independent of the number of frames, allowing for full use of full-field data. The proposed algorithm differs from conventional stress reconstruction algorithms as it possesses the necessary properties for inverse identification.

## Methodology

In this paper we propose a methodology to meet the main objectives of examining the influence of strain temporal down-sampling on the stress reconstruction and avoiding the efficiency bottleneck of conventional stress reconstruction algorithms. We conducted a simple synthetic test, an open-hole test experiment, which is commonly used in heterogeneous specimen testing. The experiment is mimicked in ABAQUS/Standard, from which the displacement field is extracted in each increment for each nodal position in FE mesh. In the next step, the extracted displacement field is used to numerically deform a reference speckle pattern, which is used to generate synthetically deformed images. These numerically deformed images are afterwards passed into DIC software to process the corresponding strain fields. In the final step, the extracted strain data is used for stress computation using both conventional Implicit Radial Mapping Algorithm (IRMA) and NICE algorithms for different sets of strain increments. For instance, the whole strain path can be integrated at each image frame or using only a selected number of frames. The difference in stress error due to the strain path approximation is monitored. This scenario is commonly encountered in real DIC application where numerous images representing incremental strain evolution are available, but stress reconstruction algorithms like VFM do not consider all increments due to efficiency limitations.

Finally, we analyze the efficiency of the IRMA and NICE algorithms when dealing with the entire strain path divided into different numbers of strain increments used in the integration process. The NICE algorithm, being explicit in nature, is expected to offer higher computational efficiency.
## Theoretical background

*Kinematics of incremental finite deformation.* In finite element computations, the equations of motion are solved numerically to obtain the nodal displacements at the end of each increment. To achieve full agreement with the kinematics of commercial FEM codes (e.g. ABAQUS) in our external code, we follow the same principles here. This step is crucial to accurately analyze the effects of the strain path when only a selected number of increments are considered in the stress integration algorithm. In other words, the input of the external code, which uses either the IRMA or the NICE algorithm, is based on the known stress state at the beginning of the increment $\sigma^n$ and the incremental deformation gradient $\mathbf{F} = \mathbf{F}^{n+1}(\mathbf{F}^n)^{-1}$, where $\mathbf{F}^n$ and $\mathbf{F}^{n+1}$ correspond to deformation gradients at the beginning and at the end of the increment with respect to the initial configuration, respectively. As implemented ABAQUS [27], a constant midpoint velocity gradient $\mathbf{L} = \mathbf{L}(t + \Delta t/2)$ and incremental deformation gradient $\mathbf{F}$ are related by $\mathbf{L} = 2(\mathbf{F} - \mathbf{1})(\mathbf{F} + \mathbf{1})^{-1}$, where **1** presents second order identity tensor. According to Hughes and Winget [27], a constant strain increment tensor can be computed as $\Delta \boldsymbol{\varepsilon} = 1/2(\mathbf{L} + \mathbf{L}^T)$ and the incremental rotation tensor as $\mathbf{R} = (\mathbf{1} - 1/2 \boldsymbol{\omega})^{-1}(\mathbf{1} + 1/2 \boldsymbol{\omega})$, with a spin tensor $\boldsymbol{\omega} = 1/2 (\mathbf{L} - \mathbf{L}^T)$. According to this procedure, the total strain is computed by approximatively integrating the strain rate according to the midpoint scheme. Since the components of the total strain and stress tensors are related to the fixed coordinate basis, the Hughes-Winget method is used to account for rigid body rotations by computing $\boldsymbol{\varepsilon}^n \leftarrow \mathbf{R} \boldsymbol{\varepsilon}^n \mathbf{R}^T$ and $\boldsymbol{\sigma}^n \leftarrow \mathbf{R} \boldsymbol{\sigma}^n \mathbf{R}^T$ before passing the states to the stress integration algorithm.

*Stress integration.* In plasticity, the evolution of stress tensor is conventionally calculated from the strain tensor evolution using radial return integration algorithm [28], which is implemented in FEM using iterative method [29]. Stress reconstruction from a measured strain field is commonly performed using the IRMA algorithm [22], based on the radial return algorithm. However, recent advances in the field have led to the development of non-iterative stress integration methods. Rossi et al. [21] proposed an approximation of the stresses without time-consuming iterations, while Yoon and Barlat [30] introduced a method that redefines the effective plastic strain increment. These approaches focus on approximating the stress state with large increments and do not address the issue of inaccurate integration caused by the coarseness of time integration and the acquired full-field data. In contrast, the algorithm described here integrates the entire strain path in small but fast increments.

The NICE stress integration algorithm is based on the work presented in [26], where we introduced an explicit scheme for integrating plasticity models that is also suitable for implicit FEM with large increments. The algorithm's stability is conditional and it marches with the maximum possible stable subincrements. Unlike other explicit schemes, the NICE scheme does not generate permanent drift from the consistency condition, making it a promising choice for stress reconstruction from measured strain fields. The NICE scheme has been successfully applied to various rate-independent plasticity models, including models with inhomogeneous yield functions, such as the Gurson model. Its accuracy, simple expressions, and fast calculations make it a highly attractive option for stress reconstruction from a measured strain field.

A comprehensive explanation of the scheme is provided in [26], which includes the formulation for the plane stress state. Although the algorithm is versatile and can be applied to different models, the formulation is simplified to von Mises plasticity with isotropic hardening for brevity.

The variables in the plane stress state von Mises plasticity model are updated using the following expressions: $\Phi = \sigma_{eq} - \sigma_Y = 0$ represents the yield condition with $\sigma_{eq}^2 = \sigma_{11}^2 + \sigma_{22}^2 + \sigma_{11}\sigma_{22} + 3\sigma_{12}^2$ being the equivalent stress and $\sigma_Y$ the current yield stress. The components of the stress deviator are defined as $s_{11} = \sigma_{11} - \sigma_H$, $s_{22} = \sigma_{22} - \sigma_H$, $s_{33} = \sigma_{33} - \sigma_H$, $s_{12} = \sigma_{12}$, where
## Theoretical background

*Kinematics of incremental finite deformation.* In finite element computations, the equations of motion are solved numerically to obtain the nodal displacements at the end of each increment. To achieve full agreement with the kinematics of commercial FEM codes (e.g. ABAQUS) in our external code, we follow the same principles here. This step is crucial to accurately analyze the effects of the strain path when only a selected number of increments are considered in the stress integration algorithm. In other words, the input of the external code, which uses either the IRMA or the NICE algorithm, is based on the known stress state at the beginning of the increment $\sigma^n$ and the incremental deformation gradient $\mathbf{F} = \mathbf{F}^{n+1}(\mathbf{F}^n)^{-1}$, where $\mathbf{F}^n$ and $\mathbf{F}^{n+1}$ correspond to deformation gradients at the beginning and at the end of the increment with respect to the initial configuration, respectively. As implemented ABAQUS [27], a constant midpoint velocity gradient $\mathbf{L} = \mathbf{L}(t + \Delta t/2)$ and incremental deformation gradient $\mathbf{F}$ are related by $\mathbf{L} = 2(\mathbf{F} - \mathbf{1})(\mathbf{F} + \mathbf{1})^{-1}$, where **1** presents second order identity tensor. According to Hughes and Winget [27], a constant strain increment tensor can be computed as $\Delta\mathbf{\epsilon} = 1/2(\mathbf{L} + \mathbf{L}^T)$ and the incremental rotation tensor as $\mathbf{R} = (\mathbf{1} - 1/2 \mathbf{\omega})^{-1}(\mathbf{1} + 1/2 \mathbf{\omega})$, with a spin tensor $\mathbf{\omega} = 1/2 (\mathbf{L} - \mathbf{L}^T)$. According to this procedure, the total strain is computed by approximatively integrating the strain rate according to the midpoint scheme. Since the components of the total strain and stress tensors are related to the fixed coordinate basis, the Hughes-Winget method is used to account for rigid body rotations by computing $\mathbf{\epsilon}^n \leftarrow \mathbf{R} \mathbf{\epsilon}^n \mathbf{R}^T$ and $\mathbf{\sigma}^n \leftarrow \mathbf{R} \mathbf{\sigma}^n \mathbf{R}^T$ before passing the states to the stress integration algorithm.

*Stress integration.* In plasticity, the evolution of stress tensor is conventionally calculated from the strain tensor evolution using radial return integration algorithm [28], which is implemented in FEM using iterative method [29]. Stress reconstruction from a measured strain field is commonly performed using the IRMA algorithm [22], based on the radial return algorithm. However, recent advances in the field have led to the development of non-iterative stress integration methods. Rossi et al. [21] proposed an approximation of the stresses without time-consuming iterations, while Yoon and Barlat [30] introduced a method that redefines the effective plastic strain increment. These approaches focus on approximating the stress state with large increments and do not address the issue of inaccurate integration caused by the coarseness of time integration and the acquired full-field data. In contrast, the algorithm described here integrates the entire strain path in small but fast increments.

The NICE stress integration algorithm is based on the work presented in [26], where we introduced an explicit scheme for integrating plasticity models that is also suitable for implicit FEM with large increments. The algorithm's stability is conditional and it marches with the maximum possible stable subincrements. Unlike other explicit schemes, the NICE scheme does not generate permanent drift from the consistency condition, making it a promising choice for stress reconstruction from measured strain fields. The NICE scheme has been successfully applied to various rate-independent plasticity models, including models with inhomogeneous yield functions, such as the Gurson model. Its accuracy, simple expressions, and fast calculations make it a highly attractive option for stress reconstruction from a measured strain field.

A comprehensive explanation of the scheme is provided in [26], which includes the formulation for the plane stress state. Although the algorithm is versatile and can be applied to different models, the formulation is simplified to von Mises plasticity with isotropic hardening for brevity.

The variables in the plane stress state von Mises plasticity model are updated using the following expressions: $\Phi = \sigma_{eq} - \sigma_Y = 0$ represents the yield condition with $\sigma_{eq}^2 = \sigma_{11}^2 + \sigma_{22}^2 + \sigma_{11}\sigma_{22} + 3\sigma_{12}^2$ being the equivalent stress and $\sigma_Y$ the current yield stress. The components of the stress deviator are defined as $s_{11} = \sigma_{11} - \sigma_H$, $s_{22} = \sigma_{22} - \sigma_H$, $s_{33} = \sigma_{33} - \sigma_H$, $s_{12} = \sigma_{12}$, where
$\sigma_H = (\sigma_{11}+\sigma_{22})/3$ is the hydrostatic stress. The elastic properties are defined as $C_{dia} = E/(1 - \nu^2)$, $C_{offdia} = C_{dia}\nu$, $\mu = E/2/(1 + \nu)$, $C_{3311} = C_{3322} = C_{offdia}(1 - \nu)/(1 - 2\nu)$, $C_{3333} = C_{3311}(\nu^{-1} - 1)$, where $E, \nu$ and $\mu$ are the elastic modulus, the Poisson's coefficient and the shear modulus, respectively. The hardening is described by the hardening modulus $H = d\sigma_Y/d\epsilon_{eq}^p$ based on the given hardening curve.

The algorithm can be outlined as follows:

1. Initialize the elastic constants $E, \nu, \mu, C_{dia}, C_{offdia}, C_{3311}, C_{3322}, C_{3333}$ at the beginning of the step and set $n = 0$.

2. Initialize variables at the beginning of the $n^{th}$ increment $\sigma^n = \{\sigma_{11}^n, \sigma_{22}^n, \sigma_{12}^n\}$, $\sigma_{eq}^n$, $\sigma_Y^n$, $(\epsilon_{eq}^p)^n$, $H^n$, $\Phi^n$ and compute $\Delta\epsilon = \{\Delta\epsilon_{11}, \Delta\epsilon_{22}, \Delta\epsilon_{12}\}$ from deformation gradients $F^n$ and $F^{n+1}$ (Eq. 3)

3. Compute trial stress state $(\boldsymbol{\sigma}^{\text{tr}} \leftarrow \boldsymbol{\sigma}^n + \Delta\boldsymbol{\sigma})$ via $\Delta\sigma_{11} = C_{\text{dia}}\Delta\epsilon_{11} + C_{\text{offdia}}\Delta\epsilon_{22}$, $\Delta\sigma_{22} = C_{\text{offdia}}\Delta\epsilon_{11} + C_{\text{dia}}\Delta\epsilon_{22}$, $\Delta\sigma_{12} = 2\mu \Delta\epsilon_{12}$, compute $\sigma_{eq}^{\text{tr}}$, and $\Phi^{\text{tr}}$.

4. If $\Phi^{\text{tr}} < 0$:
   elastic update: $\sigma^{n+1} \leftarrow \sigma^{\text{tr}}, \Phi^{n+1} \leftarrow \Phi^{\text{tr}} \Rightarrow \text{go to Step 2, next increment } n \leftarrow n + 1$
   else

elastoplastic update $\Rightarrow$ proceed with Step 5

5. Project the stress state onto the yield surface by $\beta = -\Phi^n/(\Phi^{\text{tr}} - \Phi^n)$, update the elastic stress state: $\boldsymbol{\sigma} \leftarrow \boldsymbol{\sigma}^n + \beta \Delta\boldsymbol{\sigma}$, subtract the elastic update: $\Delta\boldsymbol{\epsilon} \leftarrow (1 - \beta)\Delta\boldsymbol{\epsilon}$, $\Delta\boldsymbol{\sigma} \leftarrow (1 - \beta)\Delta\boldsymbol{\sigma}$ and compute $\sigma_{eq}, \Phi, \mathbf{s} = \{s_{11}, s_{22}, s_{33}, s_{12}\}$.

6. Compute the number of subincrements by $N_{\text{sub}} = \text{floor}((1 - \beta)V_1/(\frac{\sigma_Y}{27\mu}V_2 - \Phi) + 1)$, where $V_1 = 3(s_{11}\Delta\sigma_{11} + s_{22}\Delta\sigma_{22} + 2s_{12}\Delta\sigma_{12})/(2\sigma_{eq})$ and $V_2 = 3\mu + H - 9\mu s_{33}^2(1 - 2\nu)/(2\sigma_{eq}^2(1 - \nu))$.

7. Subincrementation and compute $\sigma_{eq}, \Phi, \mathbf{s}, H, V_3, V_4 \Rightarrow \text{Repeat this step } N_{\text{sub}}$-times

$$ \Delta\lambda = (\Phi + V_1/N_{\text{sub}})/V_2, V_3 = 3\mu/\sigma_{eq}, V_4 = s_{33}\nu/(1-\nu) $$

$$ \sigma_{11} \leftarrow \sigma_{11} + \Delta\sigma_{11}/N_{\text{sub}} - V_3(s_{11} - V_4)\Delta\lambda $$

$$ \sigma_{22} \leftarrow \sigma_{22} + \Delta\sigma_{22}/N_{\text{sub}} - V_3(s_{22} - V_4)\Delta\lambda $$

$$ \sigma_{12} \leftarrow \sigma_{12} + \Delta\sigma_{12}/N_{\text{sub}} - V_3 s_{12} \Delta\lambda $$

$$ \epsilon_{eq}^p \leftarrow \epsilon_{eq}^p + \Delta\lambda $$

$$ \sigma_Y \leftarrow \sigma_Y + H \Delta\lambda $$

8. Update variables for the next increment $\sigma^{n+1} \leftarrow \sigma$, $\sigma_Y^n \leftarrow \sigma_Y$, $\Phi^{n+1} \leftarrow \Phi$, $(\epsilon_{eq}^p)^{n+1} \leftarrow \epsilon_{eq}^p \Rightarrow \text{go to Step 2, next increment } n \leftarrow n + 1$

The number of arithmetic operations at each increment is small. There are no iterations, the scheme marches through the integration path. The size of the subincrements is usually smaller than the DIC sampling due to the stability constraint. Therefore, it is meaningless for the presented algorithm to reduce the number of time frames considered, as substepping is activated anyway for large increments; the algorithm integrates the stress path in approximately the same total number of subincrements, regardless of the time frames considered. This means that the computation time is also independent of the DIC sampling (i.e., the number of time frames measured), so it does not matter how many time frames are considered. In order to use all measured data and accurately track the load path, all sampled time frames should be considered, and the presented algorithm will always have the same computational cost for any given sampling.

## Demonstration case via synthetic experiment

In an actual experiment, the deformation gradient is sampled using DIC. To assess the performance of both stress integration algorithms when fed with DIC-sampled deformation gradients, we devised a Digital Virtual Twin (DVT) of the open hole tension test using the FEDEF module of
$\sigma_H = (\sigma_{11}+\sigma_{22})/3$ is the hydrostatic stress. The elastic properties are defined as $C_{dia} = E/(1 - \nu^2)$, $C_{offdia} = C_{dia}\nu$, $\mu = E/2/(1 + \nu)$, $C_{3311} = C_{3322} = C_{offdia}(1 - \nu)/(1 - 2\nu)$, $C_{3333} = C_{3311}(\nu^{-1} - 1)$, where $E, \nu$ and $\mu$ are the elastic modulus, the Poisson's coefficient and the shear modulus, respectively. The hardening is described by the hardening modulus $H = d\sigma_Y/d\epsilon_{eq}^p$ based on the given hardening curve.

The algorithm can be outlined as follows:

1. Initialize the elastic constants $E, \nu, \mu, C_{dia}, C_{offdia}, C_{3311}, C_{3322}, C_{3333}$ at the beginning of the step and set $n = 0$.

2. Initialize variables at the beginning of the $n^{th}$ increment $\sigma^n = \{\sigma_{11}^n, \sigma_{22}^n, \sigma_{12}^n\}$, $\sigma_{eq}^n$, $\sigma_Y^n$, $(\epsilon_{eq}^p)^n$, $H^n$, $\Phi^n$ and compute $\Delta\epsilon = \{\Delta\epsilon_{11}, \Delta\epsilon_{22}, \Delta\epsilon_{12}\}$ from deformation gradients $F^n$ and $F^{n+1}$ (Eq. 3)

3. Compute trial stress state $(\boldsymbol{\sigma}^{\text{tr}} \leftarrow \boldsymbol{\sigma}^n + \Delta\boldsymbol{\sigma})$ via $\Delta\sigma_{11} = C_{\text{dia}}\Delta\epsilon_{11} + C_{\text{offdia}}\Delta\epsilon_{22}$, $\Delta\sigma_{22} = C_{\text{offdia}}\Delta\epsilon_{11} + C_{\text{dia}}\Delta\epsilon_{22}$, $\Delta\sigma_{12} = 2\mu \Delta\epsilon_{12}$, compute $\sigma_{\text{eq}}^{\text{tr}}$, and $\Phi^{\text{tr}}$.

4. If $\Phi^{\text{tr}} < 0$:

elastic update: $\sigma^{n+1} \leftarrow \sigma^{\text{tr}}, \Phi^{n+1} \leftarrow \Phi^{\text{tr}} \Rightarrow \text{go to Step 2, next increment } n \leftarrow n + 1$

else elastoplastic update $\Rightarrow$ proceed with Step 5

5. Project the stress state onto the yield surface by $\beta = -\Phi^n/(\Phi^{\text{tr}} - \Phi^n)$, update the elastic stress state: $\boldsymbol{\sigma} \leftarrow \boldsymbol{\sigma}^n + \beta \Delta\boldsymbol{\sigma}$, subtract the elastic update: $\Delta\boldsymbol{\epsilon} \leftarrow (1 - \beta)\Delta\boldsymbol{\epsilon}$, $\Delta\boldsymbol{\sigma} \leftarrow (1 - \beta)\Delta\boldsymbol{\sigma}$ and compute $\sigma_{\text{eq}}$, $\Phi$, $\mathbf{s} = \{s_{11}, s_{22}, s_{33}, s_{12}\}$.

6. Compute the number of subincrements by $N_{\text{sub}} = \text{floor}((1 - \beta)V_1/(\frac{\sigma_Y}{27\mu}V_2 - \Phi) + 1)$, where $V_1 = 3(s_{11}\Delta\sigma_{11} + s_{22}\Delta\sigma_{22} + 2s_{12}\Delta\sigma_{12})/(2\sigma_{\text{eq}})$ and $V_2 = 3\mu + H - 9\mu s_{33}^2(1 - 2\nu)/(2\sigma_{\text{eq}}^2(1 - \nu))$.

7. Subincrementation and compute $\sigma_{\text{eq}}$, $\Phi$, $\mathbf{s}$, $H$, $V_3$, $V_4 \Rightarrow$ Repeat this step $N_{\text{sub}}$-times

$$ \Delta\lambda = (\Phi + V_1/N_{\text{sub}})/V_2, V_3 = 3\mu/\sigma_{\text{eq}}, V_4 = s_{33}\nu/(1-\nu) $$

$$ \sigma_{11} \leftarrow \sigma_{11} + \Delta\sigma_{11}/N_{\text{sub}} - V_3(s_{11} - V_4)\Delta\lambda $$

$$ \sigma_{22} \leftarrow \sigma_{22} + \Delta\sigma_{22}/N_{\text{sub}} - V_3(s_{22} - V_4)\Delta\lambda $$

$$ \sigma_{12} \leftarrow \sigma_{12} + \Delta\sigma_{12}/N_{\text{sub}} - V_3 s_{12} \Delta\lambda $$

$$ \epsilon_{eq}^p \leftarrow \epsilon_{eq}^p + \Delta\lambda $$

$$ \sigma_Y \leftarrow \sigma_Y + H \Delta\lambda $$

8. Update variables for the next increment $\sigma^{n+1} \leftarrow \boldsymbol{\sigma}$, $\sigma_Y^n \leftarrow \sigma_Y$, $\Phi^{n+1} \leftarrow \Phi$, $(\epsilon_{eq}^p)^{n+1} \leftarrow \epsilon_{eq}^p \Rightarrow$ go to Step 2, next increment $n \leftarrow n + 1$

The number of arithmetic operations at each increment is small. There are no iterations, the scheme marches through the integration path. The size of the subincrements is usually smaller than the DIC sampling due to the stability constraint. Therefore, it is meaningless for the presented algorithm to reduce the number of time frames considered, as substepping is activated anyway for large increments; the algorithm integrates the stress path in approximately the same total number of subincrements, regardless of the time frames considered. This means that the computation time is also independent of the DIC sampling (i.e., the number of time frames measured), so it does not matter how many time frames are considered. In order to use all measured data and accurately track the load path, all sampled time frames should be considered, and the presented algorithm will always have the same computational cost for any given sampling.

## Demonstration case via synthetic experiment

In an actual experiment, the deformation gradient is sampled using DIC. To assess the performance of both stress integration algorithms when fed with DIC-sampled deformation gradients, we devised a Digital Virtual Twin (DVT) of the open hole tension test using the FEDEF module of
the MatchID software. In this module, a speckle pattern is numerically deformed based on the nodal displacements generated by ABAQUS/Standard. Given that the ABAQUS simulation was divided into 300 frames, 300 deformed images were created. The optimized speckle pattern presented in [31] was used in this study. The 300 images were then post-processed using local DIC to reconstruct the deformation gradient field in each image. We should mention that in this case the location of the ROI was chosen as a function of plastic deformation occurring in the open hole tensile test. For this purpose, we avoid the regions which do not undergo plastic straining. In addition, we only consider trustworthy DIC data hence discard data affected by edge effects. In Fig. 1 the location of the considered is depicted, where 1 px 10 µm, which corresponds to typical DIC settings.

<figure><img src="image_1037.png" /><figcaption>Figure 1: The Digital Virtual Twin. (a) a numerical reference speckle pattern used to generate synthetically deformed images based on the finite element simulation of the open hole tensile test. (b)-(d) reference stress field solution computed by ABAQUS/Standard.</figcaption></figure>

During the procedure, we extracted the deformation gradient **F** for all frames and points within the ROI from the DIC analysis and repeated the process. The deformation gradient was then used as input for the external stress reconstruction software utilizing both the IRMA algorithm and the NICE algorithm. Fig. 2a through Fig. 2f display the distributions of the stress components, which were computed by IRMA and NICE from the spatiotemporal distribution of the deformation gradient **F**. From the figure it can be observed that both algorithms produce a closely matching solution. By comparing the field distributions with the referential stress state shown in the Fig. 1a to Fig. 1d, it becomes evident that solutions by both numerical algorithms also match with the solution computed by ABAQUS. Hence, it can be stated that the kinematics are adequately captured by the adopted DIC settings. Naturally, the smoothness of the numerical stress fields shown in Fig. 1 cannot be completely reproduced by the reconstructed “experimental” stress fields.
the MatchID software. In this module, a speckle pattern is numerically deformed based on the nodal displacements generated by ABAQUS/Standard. Given that the ABAQUS simulation was divided into 300 frames, 300 deformed images were created. The optimized speckle pattern presented in [31] was used in this study. The 300 images were then post-processed using local DIC to reconstruct the deformation gradient field in each image. We should mention that in this case the location of the ROI was chosen as a function of plastic deformation occurring in the open hole tensile test. For this purpose, we avoid the regions which do not undergo plastic straining. In addition, we only consider trustworthy DIC data hence discard data affected by edge effects. In Fig. 1 the location of the considered is depicted, where 1 px 10 µm, which corresponds to typical DIC settings.

<figure><img src="image_1042.png" /><figcaption>Figure 1: The Digital Virtual Twin. (a) a numerical reference speckle pattern used to generate synthetically deformed images based on the finite element simulation of the open hole tensile test. (b)-(d) reference stress field solution computed by ABAQUS/Standard.</figcaption></figure>

During the procedure, we extracted the deformation gradient **F** for all frames and points within the ROI from the DIC analysis and repeated the process. The deformation gradient was then used as input for the external stress reconstruction software utilizing both the IRMA algorithm and the NICE algorithm. Fig. 2a through Fig. 2f display the distributions of the stress components, which were computed by IRMA and NICE from the spatiotemporal distribution of the deformation gradient **F**. From the figure it can be observed that both algorithms produce a closely matching solution. By comparing the field distributions with the referential stress state shown in the Fig. 1a to Fig. 1d, it becomes evident that solutions by both numerical algorithms also match with the solution computed by ABAQUS. Hence, it can be stated that the kinematics are adequately captured by the adopted DIC settings. Naturally, the smoothness of the numerical stress fields shown in Fig. 1 cannot be completely reproduced by the reconstructed “experimental” stress fields.
<figure><img src="image_1038.png" /><figcaption>Figure 2: Comparison of referential and DIC-based reconstructed stress components: While the 1st row represents the stress fields computed from the reconstructed kinematic field using DIC by IRMA, the second row presents the solution obtained by NICE using the same procedure. Both schemes produce closely matching results.</figcaption></figure>

The blue and red curves in Fig. 3 relate to different number of frames used in stress reconstruction algorithm. It can be observed that time consumption increases exponentially with the number of frames for IRMA algorithm, but remains constant for the NICE algorithm. Moreover, we should mention that for both algorithms the consumed time increases proportionally with the number of data points – the higher the number of points, the higher the time consumption. Finally, as shown, the NICE scheme significantly outperforms the IRMA scheme when a large number of images is used. In other words, with the NICE, the experimentalist can achieve the same stress reconstruction accuracy at a computational cost which is independent of the considered number of frames. The speedup factor of the NICE compared to the IRMA in this particular case is approximately 10. The more images are considered (e.g. to capture highly nonlinear strain paths or material behavior), the higher the speedup factor will be. As such, it must be emphasized that the speedup factor highly depends on the considered problem. In this regard, for MT2.0 applications requiring repetitive stress reconstructions (e.g. nonlinear VFM), we consider the NICE scheme to be a breakthrough in terms of computational efficiency.
<figure><img src="image_1043.png" /><figcaption>Figure 2: Comparison of referential and DIC-based reconstructed stress components: While the 1st row represents the stress fields computed from the reconstructed kinematic field using DIC by IRMA, the second row presents the solution obtained by NICE using the same procedure. Both schemes produce closely matching results.</figcaption></figure>

The blue and red curves in Fig. 3 relate to different number of frames used in stress reconstruction algorithm. It can be observed that time consumption increases exponentially with the number of frames for IRMA algorithm, but remains constant for the NICE algorithm. Moreover, we should mention that for both algorithms the consumed time increases proportionally with the number of data points – the higher the number of points, the higher the time consumption. Finally, as shown, the NICE scheme significantly outperforms the IRMA scheme when a large number of images is used. In other words, with the NICE, the experimentalist can achieve the same stress reconstruction accuracy at a computational cost which is independent of the considered number of frames. The speedup factor of the NICE compared to the IRMA in this particular case is approximately 10. The more images are considered (e.g. to capture highly nonlinear strain paths or material behavior), the higher the speedup factor will be. As such, it must be emphasized that the speedup factor highly depends on the considered problem. In this regard, for MT2.0 applications requiring repetitive stress reconstructions (e.g. nonlinear VFM), we consider the NICE scheme to be a breakthrough in terms of computational efficiency.
<figure><img src="image_1039.png" /><figcaption>Figure 3: Total CPU time for all plastic updates of stress reconstruction using IRMA and NICE algorithms for varying number of frames. The relationship between the number of frames, number of data points and CPU time is presented in a Log-Log scale. As presented the CPU time increases exponentially for with the number of frames for IRMA as for NICE algorithm in remains constant, resulting in a speed-up factor of 10 for 300 frames.</figcaption></figure>

## Summary

To summarize, we have introduced a highly efficient explicit stress reconstruction algorithm based on the NICE scheme to reduce the computational challenges of repeated stress reconstruction within inverse identification methods such as the nonlinear VFM. The main findings of our study are as follows:

- The proposed method for stress reconstruction shows comparable accuracy to the conventional fully implicit return mapping algorithm.

- In contrast to the fully implicit approach, our method maintains consistent computational efficiency regardless of the number of images used to derive the strain fields by DIC.

- Specifically applied to the open-hole tensile test problem, our stress reconstruction method achieves a remarkable speed-up factor of 10 compared to the traditional fully implicit return mapping algorithm. It is important to emphasize that the observed speed-up factor depends on the nonlinearity of the problem. This indicates the potential for even greater speed up in scenarios with pronounced nonlinearity, intricate anisotropic plasticity and nonlinear loading conditions.

## Acknowledgement

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

## Disclaimer

The results reflect only the authors' view, and the European Commission is not responsible for any use that may be made of the information it contains.
<figure><img src="image_1044.png" /><figcaption>Figure 3: Total CPU time for all plastic updates of stress reconstruction using IRMA and NICE algorithms for varying number of frames. The relationship between the number of frames, number of data points and CPU time is presented in a Log-Log scale. As presented the CPU time increases exponentially for with the number of frames for IRMA as for NICE algorithm in remains constant, resulting in a speed-up factor of 10 for 300 frames.</figcaption></figure>

## Summary

To summarize, we have introduced a highly efficient explicit stress reconstruction algorithm based on the NICE scheme to reduce the computational challenges of repeated stress reconstruction within inverse identification methods such as the nonlinear VFM. The main findings of our study are as follows:

- The proposed method for stress reconstruction shows comparable accuracy to the conventional fully implicit return mapping algorithm.

- In contrast to the fully implicit approach, our method maintains consistent computational efficiency regardless of the number of images used to derive the strain fields by DIC.

- Specifically applied to the open-hole tensile test problem, our stress reconstruction method achieves a remarkable speed-up factor of 10 compared to the traditional fully implicit return mapping algorithm. It is important to emphasize that the observed speed-up factor depends on the nonlinearity of the problem. This indicates the potential for even greater speed up in scenarios with pronounced nonlinearity, intricate anisotropic plasticity and nonlinear loading conditions.

## Acknowledgement

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

## Disclaimer

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

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

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

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

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

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

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

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

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

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

[10] D. Lecompte, S. Cooreman, S. Coppieters, J. Vantomme, H. Sol, and D. Debruyne, "Parameter identification for anisotropic plasticity model using digital image correlation," *Eur. J. Comput. Mech.*, vol. 18, no. 3-4, pp. 393-418. https://doi.org/10.13052/EJCM.18.393-418

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

[12] S. Avril and F. Pierron, "General framework for the identification of constitutive parameters from full-field measurements in linear elasticity," *Int. J. Solids Struct.*, vol. 44, no. 14, pp. 4978-5002. https://doi.org/10.1016/j.ijsolstr.2006.12.018

[13] M. Grédiac and F. Pierron, "Applying the Virtual Fields Method to the identification of elasto-plastic constitutive parameters," *Int. J. Plast.*, vol. 22, no. 4, pp. 602-627. https://doi.org/10.1016/j.ijplas.2005.04.007

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

[15] A. Lattanzi, F. Barlat, F. Pierron, A. Marek, and M. Rossi, "Inverse identification strategies for the characterization of transformation-based anisotropic plasticity models with the non-linear VFM," *Int. J. Mech. Sci.*, p. 105422. https://doi.org/10.1016/j.ijmecsci.2020.105422
## References

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

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

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

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

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

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

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

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

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

[10] D. Lecompte, S. Cooreman, S. Coppieters, J. Vantomme, H. Sol, and D. Debruyne, "Parameter identification for anisotropic plasticity model using digital image correlation," *Eur. J. Comput. Mech.*, vol. 18, no. 3-4, pp. 393-418. https://doi.org/10.13052/EJCM.18.393-418

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

[12] S. Avril and F. Pierron, "General framework for the identification of constitutive parameters from full-field measurements in linear elasticity," *Int. J. Solids Struct.*, vol. 44, no. 14, pp. 4978-5002. https://doi.org/10.1016/j.ijsolstr.2006.12.018

[13] M. Grédiac and F. Pierron, "Applying the Virtual Fields Method to the identification of elasto-plastic constitutive parameters," *Int. J. Plast.*, vol. 22, no. 4, pp. 602-627. https://doi.org/10.1016/j.ijplas.2005.04.007

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

[15] A. Lattanzi, F. Barlat, F. Pierron, A. Marek, and M. Rossi, "Inverse identification strategies for the characterization of transformation-based anisotropic plasticity models with the non-linear VFM," *Int. J. Mech. Sci.*, p. 105422. https://doi.org/10.1016/j.ijmecsci.2020.105422
[16] B. Rahmani, I. Villemure, and M. Levesque, "Regularized virtual fields method for
mechanical properties identification of composite materials," *Comput. Methods Appl. Mech. Eng.*,
vol. 278, pp. 543–566. https://doi.org/10.1016/j.cma.2014.05.010

[17] L. Zhang et al., "Verification of a virtual fields method to extract the mechanical properties of human optic nerve head tissues in vivo," *Biomech. Model. Mechanobiol.*, vol. 16, no. 3, pp. 871–887. https://doi.org/10.1007/s10237-016-0858-2

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

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

[20] A. Marek, F. M. Davis, J.-H. Kim, and F. Pierron, "Experimental Validation of the
Sensitivity-Based Virtual Fields for Identification of Anisotropic Plasticity Models," *Exp. Mech.*,
vol. 60, no. 5, pp. 639–664. https://doi.org/10.1007/s11340-019-00575-3

[21] M. Rossi, A. Lattanzi, L. Cortese, and D. Amodio, "An approximated computational
method for fast stress reconstruction in large strain plasticity," *Int. J. Numer. Methods Eng.*, vol.
121, no. 14, pp. 3048–3065. https://doi.org/10.1002/nme.6346

[22] H. Takizawa et al., "Development of the User Subroutine Library 'Unified Material Model
Driver for Plasticity (UMMDp)' for Various Anisotropic Yield Functions," *J. Phys. Conf. Ser.*,
vol. 1063, no. 1, p. 012099. https://doi.org/10.1088/1742-6596/1063/1/012099

[23] M. Halilovič, M. Vrh, and B. Štok, "NICE—An explicit numerical scheme for efficient
integration of nonlinear constitutive equations," *Math. Comput. Simul.*, vol. 80, no. 294–313.
https://doi.org/10.1016/j.matcom.2009.06.030

[24] M. Vrh, M. Halilovič, and B. Štok, "Improved explicit integration in plasticity," *Int. J. Numer. Methods Eng.*, vol. 81, no. 7, pp. 910–938. https://doi.org/10.1002/nme.2737

[25] M. Halilovič, M. Vrh, and B. Štok, "NICE h: a higher-order explicit numerical scheme for
integration of constitutive models in plasticity," *Eng. Comput.*, vol. 29, no. 1, pp. 55–70.
https://doi.org/10.1007/s00366-011-0243-9

[26] M. Halilovic, B. Starman, M. Vrh, and B. Stok, "A robust explicit integration of elasto-
plastic constitutive models, based on simple subincrement size estimation," *Eng. Comput.*

[27] T. J. R. Hughes and J. Winget, "Finite rotation effects in numerical integration of rate
constitutive equations arising in large-deformation analysis," *Int. J. Numer. Methods Eng.*, vol. 15,
no. 12, pp. 1862–1867. https://doi.org/10.1002/nme.1620151210

[28] Y. Yamada, N. Yoshimura, and T. Sakurai, "Plastic stress-strain matrix and its application
for the solution of elastic-plastic problems by the finite element method," *Int. J. Mech. Sci.*, vol.
10, no. 5, pp. 343–354. https://doi.org/10.1016/0020-7403(68)90001-5

[29] M. C. Oliveira, J. L. Alves, B. M. Chaparro, and L. F. Menezes, "Study on the influence of
work-hardening modeling in springback prediction," *Int. J. Plast.*, vol. 23, no. 3, pp. 516–543.
https://doi.org/10.1016/j.ijplas.2006.07.003

[30] S. Yoon and F. Barlat, "Non-iterative stress integration method for anisotropic materials," *Int. J. Mech. Sci.*, vol. 242, p. 108003. https://doi.org/10.1016/j.ijmecsci.2022.108003

[31] S. Bossuyt, "Optimized Patterns for Digital Image Correlation," in *Imaging Methods for Novel Materials and Challenging Applications, Volume 3*, H. Jin, C. Sciammarella, C. Furlong, and S. Yoshida, Eds., in Conference Proceedings of the Society for Experimental Mechanics Series. New York, NY: Springer, 2013, pp. 239–248. https://doi.org/10.1007/978-1-4614-4235-6_34
[16] B. Rahmani, I. Villemure, and M. Levesque, "Regularized virtual fields method for mechanical properties identification of composite materials," *Comput. Methods Appl. Mech. Eng.*, vol. 278, pp. 543-566. https://doi.org/10.1016/j.cma.2014.05.010

[17] L. Zhang et al., "Verification of a virtual fields method to extract the mechanical properties of human optic nerve head tissues in vivo," *Biomech. Model. Mechanobiol.*, vol. 16, no. 3, pp. 871-887. https://doi.org/10.1007/s10237-016-0858-2

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

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

[20] A. Marek, F. M. Davis, J.-H. Kim, and F. Pierron, "Experimental Validation of the Sensitivity-Based Virtual Fields for Identification of Anisotropic Plasticity Models," *Exp. Mech.*, vol. 60, no. 5, pp. 639-664. https://doi.org/10.1007/s11340-019-00575-3

[21] M. Rossi, A. Lattanzi, L. Cortese, and D. Amodio, "An approximated computational method for fast stress reconstruction in large strain plasticity," *Int. J. Numer. Methods Eng.*, vol. 121, no. 14, pp. 3048-3065. https://doi.org/10.1002/nme.6346

[22] H. Takizawa et al., "Development of the User Subroutine Library 'Unified Material Model Driver for Plasticity (UMMDp)' for Various Anisotropic Yield Functions," *J. Phys. Conf. Ser.*, vol. 1063, no. 1, p. 012099. https://doi.org/10.1088/1742-6596/1063/1/012099

[23] M. Halilovič, M. Vrh, and B. Štok, "NICE—An explicit numerical scheme for efficient integration of nonlinear constitutive equations," *Math. Comput. Simul.*, vol. 80, no. 294-313. https://doi.org/10.1016/j.matcom.2009.06.030

[24] M. Vrh, M. Halilovič, and B. Štok, "Improved explicit integration in plasticity," *Int. J. Numer. Methods Eng.*, vol. 81, no. 7, pp. 910-938. https://doi.org/10.1002/nme.2737

[25] M. Halilovič, M. Vrh, and B. Štok, "NICE h: a higher-order explicit numerical scheme for integration of constitutive models in plasticity," *Eng. Comput.*, vol. 29, no. 1, pp. 55-70. https://doi.org/10.1007/s00366-011-0243-9

[26] M. Halilovic, B. Starman, M. Vrh, and B. Stok, "A robust explicit integration of elasto-plastic constitutive models, based on simple subincrement size estimation," *Eng. Comput.*

[27] T. J. R. Hughes and J. Winget, "Finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis," *Int. J. Numer. Methods Eng.*, vol. 15, no. 12, pp. 1862-1867. https://doi.org/10.1002/nme.1620151210

[28] Y. Yamada, N. Yoshimura, and T. Sakurai, "Plastic stress-strain matrix and its application for the solution of elastic-plastic problems by the finite element method," *Int. J. Mech. Sci.*, vol. 10, no. 5, pp. 343-354. https://doi.org/10.1016/0020-7403(68)90001-5

[29] M. C. Oliveira, J. L. Alves, B. M. Chaparro, and L. F. Menezes, "Study on the influence of work-hardening modeling in springback prediction," *Int. J. Plast.*, vol. 23, no. 3, pp. 516-543. https://doi.org/10.1016/j.ijplas.2006.07.003

[30] S. Yoon and F. Barlat, "Non-iterative stress integration method for anisotropic materials," *Int. J. Mech. Sci.*, vol. 242, p. 108003. https://doi.org/10.1016/j.ijmecsci.2022.108003

[31] S. Bossuyt, "Optimized Patterns for Digital Image Correlation," in *Imaging Methods for Novel Materials and Challenging Applications, Volume 3*, H. Jin, C. Sciammarella, C. Furlong, and S. Yoshida, Eds., in Conference Proceedings of the Society for Experimental Mechanics Series. New York, NY: Springer, 2013, pp. 239-248. https://doi.org/10.1007/978-1-4614-4235-6_34
# Insights from crystal-plasticity-based predictions on deformation of an Mg-10Gd extrudate

STEGLICH Dirk<fnref n="1" />,ª*, JEONG Youngung<fnref n="2" />,ᵇ and ANDAR Mohammad Omar<fnref n="3" />,ᶜ

**Keywords:** Anisotropy, Plastic Work Contour, Differential Hardening

**Abstract.** The mechanical behaviour of an extruded magnesium alloy was quantified in tension and compression along the extrusion direction. With the aid of EBSD measurements the material symmetry was examined, and discrete grain orientations were used to predict the multiaxial plastic behaviour along orientations perpendicular to the extrusion direction by means of the ΔEVPSC formalism. It was found that the material is transversely isotropic, shows a small amount of distortional hardening and modest differences in hardening between the extrusion direction and the tangential direction. A phenomenological plasticity law was calibrated based on the predictions, which allowed performing the finite-element analysis of the expansion process of a vascular stent.

## Introduction

Mg alloys are promising candidates for biodegradable medical implants. Examples of such devices include cardiovascular or temporary orthopaedic fixtures that are required in the body only for the time of healing [1]. These devices are machined from extruded rods showing anisotropic mechanical properties. Hence, additional to the commonly performed mechanical characterisation along the extrusion direction (ED), other loading directions should be considered. The specific loading scenarios expected from a stent during operation involve compression and extension along the circumferential direction, i.e., crimping and expansion (see Fig. 1) by a catheter. Material characteristics under these loading scenarios are difficult to access experimentally as the diameter of the extruded rod amounts to only a few millimetres. Hence, conventional mechanical testing along a direction other than the ED is not feasible in this case.

Several characteristics of magnesium alloys make constitutive modelling challenging:

1) Their low-symmetric hexagonal close-packed structure deformation requires twinning in addition to slip modes in order to accommodate arbitrary deformation.

2) Wrought magnesium alloys usually exhibit strong initial crystallographic texture, which may drastically evolve due to twin-induced reorientation.

These two above-mentioned characteristics may lead to a strength differential (SD) effect [2] and to plastic anisotropy [3], which may further evolve as the applied deformation increases depending on the loading condition. Hence, as mentioned earlier, mechanical tests data acquired along the ED alone cannot completely characterize the anisotropic behaviour of magnesium alloys.

In this study, to estimate the mechanical behaviour along other directions than the ED and under biaxial stress states, numerical simulations using the incremental elasto-visco-plastic self-consistent (ΔEVPSC) scheme [4] were performed. A mean-field crystal plasticity model with elasto-visco-plastic self-consistent formulation bears the advantage of accounting for the crystallographic texture and its evolution; hence the plastic anisotropy induced by dislocation slip

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

<footnote n="1">¹Institute of Material Systems Modeling, Helmholtz-Centre Hereon, Geesthacht, Germany</footnote>

<footnote n="2">²School of Materials Science and Engineering, Changwon National University, Changwon, Republic of Korea</footnote>

<footnote n="3">³GIZ, Dag-Hammarskjöld-Weg 1-5, 65760 Eschborn, Germany</footnote>

<footnote>ªdirk.steglich@hereon.de, ᵇyjeong@changwon.ac.kr, ᶜandaromar@gmail.com</footnote>
# Insights from crystal-plasticity-based predictions on deformation of an Mg-10Gd extrudate

STEGLICH Dirk<fnref n="1" />,ª*, JEONG Youngung<fnref n="2" />,ᵇ and ANDAR Mohammad Omar<fnref n="3" />,ᶜ

**Keywords:** Anisotropy, Plastic Work Contour, Differential Hardening

**Abstract.** The mechanical behaviour of an extruded magnesium alloy was quantified in tension and compression along the extrusion direction. With the aid of EBSD measurements the material symmetry was examined, and discrete grain orientations were used to predict the multiaxial plastic behaviour along orientations perpendicular to the extrusion direction by means of the ΔEVPSC formalism. It was found that the material is transversely isotropic, shows a small amount of distortional hardening and modest differences in hardening between the extrusion direction and the tangential direction. A phenomenological plasticity law was calibrated based on the predictions, which allowed performing the finite-element analysis of the expansion process of a vascular stent.

**Introduction**

Mg alloys are promising candidates for biodegradable medical implants. Examples of such devices include cardiovascular or temporary orthopaedic fixtures that are required in the body only for the time of healing [1]. These devices are machined from extruded rods showing anisotropic mechanical properties. Hence, additional to the commonly performed mechanical characterisation along the extrusion direction (ED), other loading directions should be considered. The specific loading scenarios expected from a stent during operation involve compression and extension along the circumferential direction, i.e., crimping and expansion (see Fig. 1) by a catheter. Material characteristics under these loading scenarios are difficult to access experimentally as the diameter of the extruded rod amounts to only a few millimetres. Hence, conventional mechanical testing along a direction other than the ED is not feasible in this case.

Several characteristics of magnesium alloys make constitutive modelling challenging:

1) Their low-symmetric hexagonal close-packed structure deformation requires twinning in addition to slip modes in order to accommodate arbitrary deformation.

2) Wrought magnesium alloys usually exhibit strong initial crystallographic texture, which may drastically evolve due to twin-induced reorientation.

These two above-mentioned characteristics may lead to a strength differential (SD) effect [2] and to plastic anisotropy [3], which may further evolve as the applied deformation increases depending on the loading condition. Hence, as mentioned earlier, mechanical tests data acquired along the ED alone cannot completely characterize the anisotropic behaviour of magnesium alloys.

In this study, to estimate the mechanical behaviour along other directions than the ED and under biaxial stress states, numerical simulations using the incremental elasto-visco-plastic self-consistent (ΔEVPSC) scheme [4] were performed. A mean-field crystal plasticity model with elasto-visco-plastic self-consistent formulation bears the advantage of accounting for the crystallographic texture and its evolution; hence the plastic anisotropy induced by dislocation slip

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

<footnote n="1">¹Institute of Material Systems Modeling, Helmholtz-Centre Hereon, Geesthacht, Germany</footnote>

<footnote n="2">²School of Materials Science and Engineering, Changwon National University, Changwon, Republic of Korea</footnote>

<footnote n="3">³GIZ, Dag-Hammarskjöld-Weg 1-5, 65760 Eschborn, Germany</footnote>

<footnote>ªdirk.steglich@hereon.de, ᵇyjeong@changwon.ac.kr, ᶜandaromar@gmail.com</footnote>
and twin accompanied by grain reorientation is intrinsically considered. The ΔEVPSC model is relatively efficient and capable of capturing various macro-mechanical features reflecting distinctive micro-mechanical mechanisms of magnesium alloys. The polycrystal aggregate was characterised by the initial texture of the as-received sample of magnesium alloy with 10 wt% gadolinium (Mg-10Gd) after extrusion, and the hardening parameters of the slip and twin mechanisms were calibrated by using the mechanical tests along the ED. The anisotropic and multiaxial behaviours of the as-received extruded rod were analysed in terms of the contours of equal plastic work calculated using the ΔEVPSC model. The polycrystal aggregate was subjected to a multiaxial loading in the principal strain space. Material symmetries were identified by analysing the texture at various positions of a sample. These symmetries were exploited with a phenomenological plasticity model. Its parameters were calibrated based on the ΔEVPSC-predicted mechanical behaviour in the space of tangential and axial stress components.

<figure><img src="image_1041.png" /><figcaption>Fig. 1: Vascular stent machined from an extruded rod with its cylindrical coordinates consisting of axes of r, φ, and z.</figcaption></figure>

## Material, texture and symmetry

The Mg-10Gd sample was produced by direct chill casting, followed by solid solution heat treatment and in-direct extrusion into a round bar with a 6 mm diameter. Metallographic analyses on its longitudinal section revealed a fully recrystallized grain structure with slightly elongated grains (aspect ratio of 0.57 with the long axis being 12 µm ± 4 µm).

The as-received material appears axisymmetric about the ED. To determine the presence or the absence of such symmetry in terms of the microstructure the crystallographic texture was quantified. To do so, EBSD data obtained at five locations on a circular cross section (plane normal is the ED) were analysed. The chosen five locations include the centre and its neighbouring positions at four φ angles (i.e., φ=0°, 90°, 180°, and 270°). Fig. 2 shows the (0002) pole figures corresponding to the five EBSD scans. The maximum intensity was found at an angle tilted away from the ED towards the radial direction (RD) in all five scans examined. Fig. 2 shows that the texture is heterogeneously developed in the as-received sample. Nevertheless, the overall distribution of crystal c-axis is axisymmetric in the scanned plane. Applying location-specific material behaviour to account for such heterogeneity in texture would make the analysis complicated and is not practical. Therefore, in this study, it was assumed that the overall mechanical behaviour exhibits the transverse isotropy due to the overall axisymmetric texture.
and twin accompanied by grain reorientation is intrinsically considered. The ΔEVPSC model is relatively efficient and capable of capturing various macro-mechanical features reflecting distinctive micro-mechanical mechanisms of magnesium alloys. The polycrystal aggregate was characterised by the initial texture of the as-received sample of magnesium alloy with 10 wt% gadolinium (Mg-10Gd) after extrusion, and the hardening parameters of the slip and twin mechanisms were calibrated by using the mechanical tests along the ED. The anisotropic and multiaxial behaviours of the as-received extruded rod were analysed in terms of the contours of equal plastic work calculated using the ΔEVPSC model. The polycrystal aggregate was subjected to a multiaxial loading in the principal strain space. Material symmetries were identified by analysing the texture at various positions of a sample. These symmetries were exploited with a phenomenological plasticity model. Its parameters were calibrated based on the ΔEVPSC-predicted mechanical behaviour in the space of tangential and axial stress components.

<figure><img src="image_1046.png" /><figcaption>Fig. 1: Vascular stent machined from an extruded rod with its cylindrical coordinates consisting of axes of r, φ, and z.</figcaption></figure>

## Material, texture and symmetry

The Mg-10Gd sample was produced by direct chill casting, followed by solid solution heat treatment and in-direct extrusion into a round bar with a 6 mm diameter. Metallographic analyses on its longitudinal section revealed a fully recrystallized grain structure with slightly elongated grains (aspect ratio of 0.57 with the long axis being 12 µm ± 4 µm).

The as-received material appears axisymmetric about the ED. To determine the presence or the absence of such symmetry in terms of the microstructure the crystallographic texture was quantified. To do so, EBSD data obtained at five locations on a circular cross section (plane normal is the ED) were analysed. The chosen five locations include the centre and its neighbouring positions at four φ angles (i.e., φ=0°, 90°, 180°, and 270°). Fig. 2 shows the (0002) pole figures corresponding to the five EBSD scans. The maximum intensity was found at an angle tilted away from the ED towards the radial direction (RD) in all five scans examined. Fig. 2 shows that the texture is heterogeneously developed in the as-received sample. Nevertheless, the overall distribution of crystal c-axis is axisymmetric in the scanned plane. Applying location-specific material behaviour to account for such heterogeneity in texture would make the analysis complicated and is not practical. Therefore, in this study, it was assumed that the overall mechanical behaviour exhibits the transverse isotropy due to the overall axisymmetric texture.
<figure><img src="image_1042.png" /><figcaption>Fig. 2: Basal pole figures obtained at different locations of the extruded rod section.</figcaption></figure>

## Prediction of iso-plastic-work contours

The ΔEVPSC framework was used to predict the plastic behaviour in the stress space. For this, the Voce hardening parameters pertaining to basal, prismatic and pyramidal slip modes as well as (extension) twinning mode were calibrated by fitting the uniaxial tensile and the compressive flow responses of the polycrystalline aggregate with the experiments conducted in the ED. The uniaxial tension and compression simulations were performed using a set of 5000 discrete orientations from one of the EBSD scans shown in Fig. 2.

The plastic flow and the hardening behavior obtained from various biaxial loading simulations can be systematically characterized with respect to the specific plastic work $W_p$. Out of three spatial coordinate axes (r, φ, z), the three different biaxial stress spaces of interest were considered: 1) $σ_r$ vs. $σ_φ$, 2) $σ_z$ vs $σ_r$, and 3) $σ_φ$ vs. $σ_z$. In each 2D biaxial stress space, evolution in the stress components corresponding to a certain plastic work during a proportional strain-ratio-controlled loadings were linearly interpolated. The calculations were repeated under various strain ratios to obtain an iso-work contour as was done in [5]. All the stress coordinates were normalized with respect to its uniaxial tension flow stress, and the results are shown in Fig. 3.
<figure><img src="image_1047.png" /><figcaption>Fig. 2: Basal pole figures obtained at different locations of the extruded rod section.</figcaption></figure>

## Prediction of iso-plastic-work contours

The ΔEVPSC framework was used to predict the plastic behaviour in the stress space. For this, the Voce hardening parameters pertaining to basal, prismatic and pyramidal slip modes as well as (extension) twinning mode were calibrated by fitting the uniaxial tensile and the compressive flow responses of the polycrystalline aggregate with the experiments conducted in the ED. The uniaxial tension and compression simulations were performed using a set of 5000 discrete orientations from one of the EBSD scans shown in Fig. 2.

The plastic flow and the hardening behavior obtained from various biaxial loading simulations can be systematically characterized with respect to the specific plastic work $W_p$. Out of three spatial coordinate axes (r, φ, z), the three different biaxial stress spaces of interest were considered: 1) $σ_r$ vs. $σ_φ$, 2) $σ_z$ vs $σ_r$, and 3) $σ_φ$ vs. $σ_z$. In each 2D biaxial stress space, evolution in the stress components corresponding to a certain plastic work during a proportional strain-ratio-controlled loadings were linearly interpolated. The calculations were repeated under various strain ratios to obtain an iso-work contour as was done in [5]. All the stress coordinates were normalized with respect to its uniaxial tension flow stress, and the results are shown in Fig. 3.
<figure><img src="image_1043.png" /><figcaption>Fig. 3: ΔEVPSC-predicted contours of equal plastic work.</figcaption></figure>

## Phenomenological plasticity model

Mean field crystal plasticity models can be used for solving boundary-value problems, from which the plastic work-contours in Fig. 3 were obtained. However, its application to any complex geometry via finite element simulation may not be practical owing to its low computational speed compared to phenomenological plasticity models. In the current investigation, a chosen phenomenological plasticity model was characterized by the results of the ΔEVPSC simulations. Among various phenomenological models, the CPB model developed by Cazacu et al. [6] was selected, which describes the plastic potential as a function of linearly transformed stress tensor via:

$$ \tilde{B}^a(\boldsymbol{\Sigma}) = (|\Sigma_1| - k\Sigma_1)^a + (|\Sigma_2| - k\Sigma_2)^a + (|\Sigma_3| - k\Sigma_3)^a, \quad (1) $$

where $\Sigma_1$, $\Sigma_2$, and $\Sigma_3$ are the principal values of the transformed stress deviator. The transformed stress deviator $\boldsymbol{\Sigma}$ follows from below linear transformation:

$$ \boldsymbol{\Sigma} = \boldsymbol{C} : \boldsymbol{\sigma}. \quad (2) $$

where $\boldsymbol{\sigma}$ and $\boldsymbol{C}$ are the Cauchy stress tensor and the linear transformation matrix. Using the Voigt notation of stress tensor Eq. (2) can be expressed as the following matrix form:
<figure><img src="image_1048.png" /><figcaption>Fig. 3: ΔEVPSC-predicted contours of equal plastic work.</figcaption></figure>

## Phenomenological plasticity model

Mean field crystal plasticity models can be used for solving boundary-value problems, from which the plastic work-contours in Fig. 3 were obtained. However, its application to any complex geometry via finite element simulation may not be practical owing to its low computational speed compared to phenomenological plasticity models. In the current investigation, a chosen phenomenological plasticity model was characterized by the results of the ΔEVPSC simulations. Among various phenomenological models, the CPB model developed by Cazacu et al. [6] was selected, which describes the plastic potential as a function of linearly transformed stress tensor via:

$$ \tilde{B}^a(\boldsymbol{\Sigma}) = (|\Sigma_1| - k\Sigma_1)^a + (|\Sigma_2| - k\Sigma_2)^a + (|\Sigma_3| - k\Sigma_3)^a, \quad (1) $$

where $\Sigma_1$, $\Sigma_2$, and $\Sigma_3$ are the principal values of the transformed stress deviator. The transformed stress deviator $\boldsymbol{\Sigma}$ follows from below linear transformation:

$$ \boldsymbol{\Sigma} = \boldsymbol{C}:\boldsymbol{\sigma}. \quad (2) $$

where $\boldsymbol{\sigma}$ and $\boldsymbol{C}$ are the Cauchy stress tensor and the linear transformation matrix. Using the Voigt notation of stress tensor Eq. (2) can be expressed as the following matrix form:
$$ \Sigma = \begin{pmatrix} c_{11} & c_{12} & c_{13} & 0 & 0 & 0 \\ c_{12} & c_{11} & c_{13} & 0 & 0 & 0 \\ c_{13} & c_{13} & c_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & 1/2(c_{11} - c_{12}) & 0 & 0 \\ 0 & 0 & 0 & 0 & c_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & c_{66} \end{pmatrix} \begin{pmatrix} \sigma_r \\ \sigma_\varphi \\ \sigma_z \\ \sigma_{r\varphi} \\ \sigma_{rz} \\ \sigma_{\varphi z} \end{pmatrix}. \quad (3) $$

Due to the global axisymmetry in the texture, a typical loading condition of stents in service would have zero shear components, $\sigma_{rz} = \sigma_{\varphi z} = 0$ in Eq. 3. The latter implies that the parameters $c_{55}$ and $c_{66}$ in Eq. (3) are arbitrary and were set to unity. The remaining four coefficients, $c_{11}$, $c_{12}$, $c_{13}$ and $c_{33}$, were determined by fitting the stress-strain responses obtained from $\Delta$EVPSC simulations for tension and compression along the $\varphi (=r)$- and $z$ direction. Isotropic hardening based on two exponential functions was used:

$$ \sigma_Y = R_0 + Q_1 \exp(-b_1 p) + Q_2 \exp(-b_2 p), \quad (4) $$

where $p$ is the equivalent plastic strain; and $R_0$, $Q_1$, $Q_2$, $b_1$ and $b_2$ are fitting parameters. The results are depicted in the left part of Fig. 4. Right in Fig. 4, the computed initial yield surfaces in the stress space of $(\sigma_\varphi, \sigma_z)$ are shown. The hardening parameters and the transformation matrix $C$ are given in Table 1.

<table><caption>Table 1: Model parameters of the phenomenological plasticity model.</caption><thead><tr><th>R<sub>0</sub> [MPa]</th><th>Q<sub>1</sub> [MPa]</th><th>b<sub>1</sub></th><th>Q<sub>2</sub> [MPa]</th><th>b<sub>2</sub></th><th>k</th><th>c<sub>11</sub></th><th>c<sub>12</sub></th><th>c<sub>13</sub></th><th>c<sub>33</sub></th></tr></thead><tbody><tr><td>173.9</td><td>188.6</td><td>14.6</td><td>11.2</td><td>226.6</td><td>0.083</td><td>1.42</td><td>-0.26</td><td>1.34</td><td>0.11</td></tr></tbody></table>

<figure><img src="image_1044.png" /><figcaption>Fig. 4: Flow stress vs. strain curves in tension and compression fitted using the phenomenological plasticity model (left); and the associated initial yield surface (right) in comparison with the $\Delta$EVPSC-predicted initial plastic work-contour.</figcaption></figure>
$$ \Sigma = \begin{pmatrix} c_{11} & c_{12} & c_{13} & 0 & 0 & 0 \\ c_{12} & c_{11} & c_{13} & 0 & 0 & 0 \\ c_{13} & c_{13} & c_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & 1/2(c_{11} - c_{12}) & 0 & 0 \\ 0 & 0 & 0 & 0 & c_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & c_{66} \end{pmatrix} \begin{pmatrix} \sigma_r \\ \sigma_\varphi \\ \sigma_z \\ \sigma_{r\varphi} \\ \sigma_{rz} \\ \sigma_{\varphi z} \end{pmatrix}. \quad (3) $$

Due to the global axisymmetry in the texture, a typical loading condition of stents in service would have zero shear components, $\sigma_{rz} = \sigma_{\varphi z} = 0$ in Eq. 3. The latter implies that the parameters $c_{55}$ and $c_{66}$ in Eq. (3) are arbitrary and were set to unity. The remaining four coefficients, $c_{11}$, $c_{12}$, $c_{13}$ and $c_{33}$, were determined by fitting the stress-strain responses obtained from $\Delta$EVPSC simulations for tension and compression along the $\varphi$ (=r)- and z direction. Isotropic hardening based on two exponential functions was used:

$$ \sigma_Y = R_0 + Q_1 \exp(-b_1 p) + Q_2 \exp(-b_2 p), \quad (4) $$

where $p$ is the equivalent plastic strain; and $R_0$, $Q_1$, $Q_2$, $b_1$ and $b_2$ are fitting parameters. The results are depicted in the left part of Fig. 4. Right in Fig. 4, the computed initial yield surfaces in the stress space of $(\sigma_\varphi, \sigma_z)$ are shown. The hardening parameters and the transformation matrix $C$ are given in Table 1.

<table><caption>Table 1: Model parameters of the phenomenological plasticity model.</caption><thead><tr><th>R<sub>0</sub> [MPa]</th><th>Q<sub>1</sub> [MPa]</th><th>b<sub>1</sub></th><th>Q<sub>2</sub> [MPa]</th><th>b<sub>2</sub></th><th>k</th><th>c<sub>11</sub></th><th>c<sub>12</sub></th><th>c<sub>13</sub></th><th>c<sub>33</sub></th></tr></thead><tbody><tr><td>173.9</td><td>188.6</td><td>14.6</td><td>11.2</td><td>226.6</td><td>0.083</td><td>1.42</td><td>-0.26</td><td>1.34</td><td>0.11</td></tr></tbody></table>

<figure><img src="image_1049.png" /><figcaption>Fig. 4: Flow stress vs. strain curves in tension and compression fitted using the phenomenological plasticity model (left); and the associated initial yield surface (right) in comparison with the $\Delta$EVPSC-predicted initial plastic work-contour.</figcaption></figure>
## Results and discussion

The numerical simulations using the ΔEVPSC-model led to the following conclusions:

* The mechanical property of the as-received Mg-10Gd rod exhibits the transverse isotropy on its cross section, which is perpendicular to the ED. The initial iso-work contour in the stress space of $(\sigma_r, \sigma_\varphi)$ is consistent with such symmetry in that the uniaxial tension stresses along $r$ and $\varphi$ axes are virtually the same.

* The strain hardening behaviour shown in Fig. 4. deviates from the isotropic hardening law:

i. In the case of tension, while the initial yield stresses along the $\varphi$ and the $z$ axes are the same, the flow stress $\sigma_\varphi$ becomes noticeably higher than $\sigma_z$.

ii. Likewise, compression yield stresses along the $\varphi$ and the $z$ axes are the same. However, the flow stress $\sigma_z$ becomes noticeably higher than $\sigma_\varphi$.

iii. An increased hardening rate is related to a reduction in prismatic slip, which is compensated by a rise of pyramidal slip.

iv. In equi-biaxial tension, an increased hardening rate is observed, which is related to the decreasing activity of the prismatic slip systems for crystallographic orientations with the c-axis along the ED [7].

* According to the work contour in the $(\sigma_\varphi, \sigma_z)$ plane, the tensile flow stress level in the tangential direction (TD), which aligns with the $\varphi$ axis, is slightly higher than that in the ED. Under uniaxial compression, the situation is reversed. Consequently, the predicted SD effect, quantified by the ratio of yield stresses in tension and compression $(\sigma_t/\sigma_c)$, becomes 1.01 and 0.98 in the TD and the ED, respectively.

* The plasticity law was tailored for this transversely isotropic material while accounting its SD-effect. As a consequence, the yield locus characterized by the phenomenological model in the $(\sigma_r, \sigma_\varphi)$ plane are symmetric about the line $\sigma_r = \sigma_\varphi$, and cannot capture the dissimilarity in the SD effect along the $\sigma_r$ and $\sigma_\varphi$ axes in the stress space calculated by the above-mentioned ΔEVPSC model prediction.

* The respective parameters were determined by fitting the flow curves of the $(\sigma_\varphi, \sigma_z)$ plane in the interval of total strain [0; 0.2]. The anisotropy could be well described. Changes in the mechanisms of deformation due to grain rotation, however, remain a challenge.

The strain hardening features mentioned above in i.-iv. constitute the distortional hardening law. However, the degree of distortion is not pronounced. Therefore, for the application envisaged, distortional hardening was not further accounted for in the phenomenological model considered in this study.

It is frequently reported that adding rare earth elements weakens the texture, which results in reducing anisotropy and tension-compression asymmetry and an increased ductility, e.g. [8]. This could be confirmed for the material under investigation. It remains in the responsibility of the design engineers to decide, when numerically assessing the performance of implants, whether to use a simple isotropic plasticity law or a more complex model that accounts for anisotropy and the SD effect. While the latter might not be essentially required in the present case considering the computational efficiency, other (stronger textured) magnesium alloys might demand for a careful assessment to determine an adequate model to use. In the case of the 3D-simulation of the crimping (reduction of the stent diameter prior to implementation) and expansion (radial extension by an inflating balloon catheter) processes of a typical vascular stent shown in Fig. 5, the computational time increases by a factor of 2.6 if the anisotropic material model is used (ABAQUS software with UMAT interface, 2.2 million degrees of freedom) compared to the case where von Mises plasticity model is considered. The simulation techniques and boundary conditions for the prediction of stents behavior are subject to ongoing research, see e.g. [9].
## Results and discussion

The numerical simulations using the ΔEVPSC-model led to the following conclusions:

* The mechanical property of the as-received Mg-10Gd rod exhibits the transverse isotropy on its cross section, which is perpendicular to the ED. The initial iso-work contour in the stress space of $(\sigma_r, \sigma_\varphi)$ is consistent with such symmetry in that the uniaxial tension stresses along $r$ and $\varphi$ axes are virtually the same.

* The strain hardening behaviour shown in Fig. 4. deviates from the isotropic hardening law:

i. In the case of tension, while the initial yield stresses along the $\varphi$ and the $z$ axes are the same, the flow stress $\sigma_\varphi$ becomes noticeably higher than $\sigma_z$.

ii. Likewise, compression yield stresses along the $\varphi$ and the $z$ axes are the same. However, the flow stress $\sigma_z$ becomes noticeably higher than $\sigma_\varphi$.

iii. An increased hardening rate is related to a reduction in prismatic slip, which is compensated by a rise of pyramidal slip.

iv. In equi-biaxial tension, an increased hardening rate is observed, which is related to the decreasing activity of the prismatic slip systems for crystallographic orientations with the c-axis along the ED [7].

* According to the work contour in the $(\sigma_\varphi, \sigma_z)$ plane, the tensile flow stress level in the tangential direction (TD), which aligns with the $\varphi$ axis, is slightly higher than that in the ED. Under uniaxial compression, the situation is reversed. Consequently, the predicted SD effect, quantified by the ratio of yield stresses in tension and compression $(\sigma_t/\sigma_c)$, becomes 1.01 and 0.98 in the TD and the ED, respectively.

* The plasticity law was tailored for this transversely isotropic material while accounting its SD-effect. As a consequence, the yield locus characterized by the phenomenological model in the $(\sigma_r, \sigma_\varphi)$ plane are symmetric about the line $\sigma_r = \sigma_\varphi$, and cannot capture the dissimilarity in the SD effect along the $\sigma_r$ and $\sigma_\varphi$ axes in the stress space calculated by the above-mentioned ΔEVPSC model prediction.

* The respective parameters were determined by fitting the flow curves of the $(\sigma_\varphi, \sigma_z)$ plane in the interval of total strain [0; 0.2]. The anisotropy could be well described. Changes in the mechanisms of deformation due to grain rotation, however, remain a challenge.

The strain hardening features mentioned above in i.-iv. constitute the distortional hardening law. However, the degree of distortion is not pronounced. Therefore, for the application envisaged, distortional hardening was not further accounted for in the phenomenological model considered in this study.

It is frequently reported that adding rare earth elements weakens the texture, which results in reducing anisotropy and tension-compression asymmetry and an increased ductility, e.g. [8]. This could be confirmed for the material under investigation. It remains in the responsibility of the design engineers to decide, when numerically assessing the performance of implants, whether to use a simple isotropic plasticity law or a more complex model that accounts for anisotropy and the SD effect. While the latter might not be essentially required in the present case considering the computational efficiency, other (stronger textured) magnesium alloys might demand for a careful assessment to determine an adequate model to use. In the case of the 3D-simulation of the crimping (reduction of the stent diameter prior to implementation) and expansion (radial extension by an inflating balloon catheter) processes of a typical vascular stent shown in Fig. 5, the computational time increases by a factor of 2.6 if the anisotropic material model is used (ABAQUS software with UMAT interface, 2.2 million degrees of freedom) compared to the case where von Mises plasticity model is considered. The simulation techniques and boundary conditions for the prediction of stents behavior are subject to ongoing research, see e.g. [9].
<figure><img src="image_1045.png" /><figcaption>Fig. 5: Equivalent plastic strain in stent segments after crimping and expansion (FE simulation).</figcaption></figure>

## Acknowledgements

The authors thank Henry Ovri for conducting the EBSD scans, and Falk Dorn for sample preparation. The fruitful discussion with Holger Aretz regarding symmetry groups is appreciated.

## References

[1] Staiger MP, Pietak AM, Huadmai J, Dias G. Magnesium and its alloys as orthopedic biomaterials: A review. *Biomaterials*. 2006;27(9):1728-34.
https://doi.org/10.1016/j.biomaterials.2005.10.003

[2] Hosford WF, Allen TJ. Twinning and directional slip as a cause for a strength differential effect. *Metallurgical Transactions*. 1973;4(5):1424-5. https://doi.org/10.1007/BF02644545

[3] Kelley EW, Hosford WF. The deformation characteristics of textured magnesium. *Trans Metall Soc AIME*. 1968;242:654-61.

[4] Jeong Y, Tomé CN. An efficient elasto-visco-plastic self-consistent formulation: Application to steel subjected to loading path changes. *International Journal of Plasticity*. 2020;135:102812. https://doi.org/10.1016/j.ijplas.2020.102812

[5] Steglich D, Jeong Y, Andar MO, Kuwabara T. Biaxial deformation behaviour of AZ31 magnesium alloy: Crystal-plasticity-based prediction and experimental validation. *International Journal of Solids and Structures*. 2012;49(25):3551-61.
https://doi.org/10.1016/j.ijsolstr.2012.06.017

[6] Cazacu O, Plunkett B, Barlat F. Orthotropic yield criterion for hexagonal closed packed metals. *Int J Plast*. 2006;22:1171-94. https://doi.org/10.1016/j.ijplas.2005.06.001

[7] Hama T, Takuda H. Crystal plasticity finite-element simulation of work-hardening behavior in a magnesium alloy sheet under biaxial tension. *Computational Materials Science*. 2012;51(1):156-64. https://doi.org/10.1016/j.commatsci.2011.07.026

[8] Ball EA, Prangnell PB. Tensile-compressive yield asymmetries in high strength wrought magnesium alloys. *Scripta Metallurgica et Materiala*. 1994;31(2):111-6.
https://doi.org/10.1016/0956-716X(94)90159-7

[9] Wiesent L, Schultheiß U, Schmid C, Schratzenstaller T, Nonn A. Experimentally validated simulation of coronary stents considering different dogboning ratios and asymmetric stent positioning. *PLOS ONE*. 2019;14(10):e0224026. https://doi.org/10.1371/journal.pone.0224026
<figure><img src="image_1050.png" /><figcaption>Fig. 5: Equivalent plastic strain in stent segments after crimping and expansion (FE simulation).</figcaption></figure>

## Acknowledgements

The authors thank Henry Ovri for conducting the EBSD scans, and Falk Dorn for sample preparation. The fruitful discussion with Holger Aretz regarding symmetry groups is appreciated.

## References

[1] Staiger MP, Pietak AM, Huadmai J, Dias G. Magnesium and its alloys as orthopedic biomaterials: A review. *Biomaterials*. 2006;27(9):1728-34.
https://doi.org/10.1016/j.biomaterials.2005.10.003

[2] Hosford WF, Allen TJ. Twinning and directional slip as a cause for a strength differential effect. *Metallurgical Transactions*. 1973;4(5):1424-5. https://doi.org/10.1007/BF02644545

[3] Kelley EW, Hosford WF. The deformation characteristics of textured magnesium. *Trans Metall Soc AIME*. 1968;242:654-61.

[4] Jeong Y, Tomé CN. An efficient elasto-visco-plastic self-consistent formulation: Application to steel subjected to loading path changes. *International Journal of Plasticity*. 2020;135:102812. https://doi.org/10.1016/j.ijplas.2020.102812

[5] Steglich D, Jeong Y, Andar MO, Kuwabara T. Biaxial deformation behaviour of AZ31 magnesium alloy: Crystal-plasticity-based prediction and experimental validation. *International Journal of Solids and Structures*. 2012;49(25):3551-61.
https://doi.org/10.1016/j.ijsolstr.2012.06.017

[6] Cazacu O, Plunkett B, Barlat F. Orthotropic yield criterion for hexagonal closed packed metals. *Int J Plast*. 2006;22:1171-94. https://doi.org/10.1016/j.ijplas.2005.06.001

[7] Hama T, Takuda H. Crystal plasticity finite-element simulation of work-hardening behavior in a magnesium alloy sheet under biaxial tension. *Computational Materials Science*. 2012;51(1):156-64. https://doi.org/10.1016/j.commatsci.2011.07.026

[8] Ball EA, Prangnell PB. Tensile-compressive yield asymmetries in high strength wrought magnesium alloys. *Scripta Metallurgica et Materiala*. 1994;31(2):111-6.
https://doi.org/10.1016/0956-716X(94)90159-7

[9] Wiesent L, Schultheiß U, Schmid C, Schratzenstaller T, Nonn A. Experimentally validated simulation of coronary stents considering different dogboning ratios and asymmetric stent positioning. *PLOS ONE*. 2019;14(10):e0224026. https://doi.org/10.1371/journal.pone.0224026
# Tension-compression asymmetry of ferritic stainless steel sheet

Sekiya Koki¹,ª, Higashii Hyuga¹,ᵇ, Taoka Akira²,ᶜ and Kuwabara Toshihiko³,ᵈ*

¹Department of Mechanical Systems Engineering, Graduate School of Engineering, Tokyo
University of Agriculture and Technology, 2-24-16 Naka-cho, Koganei-shi, Tokyo 184-8588,
Japan

²Raw Material Production Engineering Development Section, AISIN Co., 2-1 Asahi-cho, Kariya-
shi, Aichi 448-8650, Japan

³Division of Advanced Mechanical Systems Engineering, Institute of Engineering, Tokyo
University of Agriculture and Technology, 2-24-16 Naka-cho, Koganei-shi, Tokyo 184-8588,
Japan

ªsekiya@st.go.tuat.ac.jp, ᵇs209914v@st.go.tuat.ac.jp, ᶜa taoka@koki.aisin.co.jp,
ᵈkuwabara@cc.tuat.ac.jp

**Keywords:** In-Plane Compression, Strength Differential Effect, Bending Moment

**Abstract.** The strength differential effect (SDE) (i.e., the tension-compression asymmetry of stress-strain curves), for a ferritic stainless steel sheet was measured using an in-plane compression test. It was found that the in-plane compressive flow stress was 4-7% higher than the uniaxial tensile flow stress in both the rolling and transverse directions. In addition, bending moment-curvature curves were obtained from a pure bending test. These curves were found to be in good agreement with those calculated with SDE taken into account. Thus, the validity of the measured SDE was confirmed by the bending experiment.

**Introduction**

Stainless steel sheets are used for automotive exhaust system components because of their superior corrosion and heat resistance. Among the various types of stainless steel, ferritic stainless steel, which is less expensive than austenitic stainless steel and has superior thermal fatigue characteristics, has seen increasing use [1]. Automotive parts are often produced by press forming, in which the materials are subjected to tensile and compressive stresses. The strength differential effect (SDE) has been observed in steel sheets [2]-[9]. Kuwahara et al. [10] measured the SDE in austenitic stainless steel sheets used for electronic components and quantitatively clarified its effect on the prediction of springback in bending. Noma et al. [8] improved the prediction accuracy of springback in forming simulations of curved panels by using a material model that takes into account the SDE of a high-strength steel sheet. Yoon et al. [11] implemented the hydrostatic pressure dependency of the distortional plasticity model HAH₂₀ [12] using a finite element code to account for the SDE observed in advanced high-strength steel sheets. The implementation was validated by comparing the finite element prediction of the material behavior during U-draw bending and B-pillar forming. However, to the authors' knowledge, no studies have measured the SDE of ferritic stainless steel sheets.

The objective of this study is to measure the SDE of a 1.0-mm-thick Fe-18Cr-3Al ferritic stainless steel sheet. First, uniaxial tensile tests and in-plane compression tests (IPCTs) using an in-plane stress reversal test apparatus [7] are conducted to quantitatively evaluate the magnitude of the SDE. Next, bending moment-curvature (M - κ) curves are obtained using a pure bending test apparatus. It is found that the experimental M - κ curves are consistent with those calculated using the tensile and compressive stress-strain (ss) curves in both the rolling direction (RD) and

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