# Review of “A non-local model for the description of twinning in polycrystalline materials in the context of infinitesimal strains: application to a magnesium alloy”

Charles Mareau, Hamidreza Abdolvand, Nicolo Grilli, Shailendra Joshi,
Laurence Brassart

► To cite this version:

Charles Mareau, Hamidreza Abdolvand, Nicolo Grilli, Shailendra Joshi, Laurence Brassart. Review of “A non-local model for the description of twinning in polycrystalline materials in the context of infinitesimal strains: application to a magnesium alloy”. 2022. hal-03603524

HAL Id: hal-03603524

https://hal.science/hal-03603524

Submitted on 15 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Identifiers

Open Review

OAI hal-03603524

Reviewed Article

DOI 10.46298/jtcam.7562

History

Review 9 Jun, 2021

Rebuttal 29 Aug, 2021

Review 2 25 Jan, 2022

Rebuttal 2 6 Feb, 2022

Licence

CC BY 4.0

©The Authors

# Review of “A non-local model for the description of twinning in polycrystalline materials in the context of infinitesimal strains: application to a magnesium alloy”

Charles MAREAU¹, Hamidreza ABDOLVAND², Nicolo GRILLI³, R,

Shailendra Joshi ⁴, R, and Laurence BRASSART⁵, E

¹ Arts et Métiers Institute of Technology, LAMPA, HESAM Université, Angers, France

² Department of Mechanical and Materials Engineering, Western University, London, Canada

³ Department of Mechanical Engineering, University of Bristol, United Kingdom

⁴ Department of Mechanical Engineering, University of Houston, USA

⁵ Department of Engineering Science, University of Oxford, United Kingdom

R Reviewer

E Editor

## Review of version 1

Permalink: arXiv.org:2012.13156v3

**Authors** Dear Editor, we would like to thank the reviewers for their thorough reading of the manuscript. We believe we have addressed the comments raised by the reviewers in the revised version. A detailed response is provided below.

## Reviewer 1 (Nicolo GRILLI)

**Reviewer** Dear Editor, the paper is about twin modelling in polycrystalline AZ31 magnesium using a phase field model coupled with crystal plasticity FEM. The paper is well written and the results are interesting. However, many aspects of the model and the results must be clarified. Therefore, I recommend that major revisions of the paper are carried out before publication. Please find my specific comments below.

There is no novelty statement in this paper, the novelty itself may be very limited but you need to state it anyway. Has the twinning model together with the spectral solver been previously used?

**Authors** The main novelty of the proposed model is that, in contrast with previous implementations of twinning in the context of polycrystalline plasticity within a spectral solver (e.g. Mareau and Daymond 2016, Paramatmuni, and Kanjarla 2019), non-locality is introduced to circumvent the difficulties associated with plastic strain localization resulting from twinning-induced softening. This is now stated in the introduction section.

**Reviewer** In the introduction, neutron and X-ray diffraction to study twinning are mentioned, but no mention of in-situ EBSD, which has become a popular technique to study twinning in the recent years. See for example:

* Nicolò Grilli, Edmund Tarleton, Philip Edmondson, Maxim Gussev, and Alan Cocks, *In situ measurement and modelling of the growth and length scale of twins in alpha-uranium*. Phys. Rev. Materials 4, 043605, 2020. 10.1103/PhysRevMaterials.4.043605

* M.N. Gussev, P.D. Edmondson, K.J. Leonard, *Beam current effect as a potential challenge in SEM-EBSD in situ tensile testing*. Materials Characterization, 146:25-34, 2018. 10.1016/j.matchar.2018.09.037.
**Authors** The possibilities offered by *in situ* EBSD techniques are now mentioned in the introduction section. The above references have also been included.

**Reviewer** Model Section 2: my general comment is that lot of equations are presented in an attempt to show all the derivations. However, phase field theory and the free energy formulations are standard and the authors should try to keep only the equations that are essential to understand the evolution of the state variables. Also, the numerical scheme used to implement the model equations should be reported. A figure showing the steps of the numerical scheme is needed. I would also report the free energy first and then the dissipated power and not the opposite.

**Authors** An appendix has been included in the revised version to present the numerical scheme used for the resolution of PDEs with the spectral method (see Algorithm 1). For the integration of constitutive equations, a fourth order Runge Kutta method is used. Also, in the original submission the free energy is presented first (Section 2.3) and the dissipated power (Section 2.4) immediately after. I guess the reviewer refers to equation (27), which shows the time derivative of the free energy. In the revised version, this equation has been moved to Section 2.3.

1. Set the state variables to their initial values

2. Compute the modified Green tensor in frequency domain

3. **for** each increment *j* of a deformation process **do**

4. Compute time at the end of the increment

5. $j \leftarrow 1$

6. **while** $e_{bc} > e_{bc,c}$ **do**

7. Estimate the macroscopic strain and rotation tensors from boundary conditions

8. **while** $e_{eq} > e_{eq,c}$ **do**

9. **if** $i = 1$ and $j = 1$ **then**

10. Use Voigt approximation to obtain initial estimates for strain and rotation tensor fields

11. **else**

12. Use spectral method to obtain new estimates for strain and rotation tensor fields

13. **end if**

14. **for** each voxel **do**

15. **for** each child or parent constituent **do**

16. Use the Sachs approximation to obtain the corresponding strain tensor

17. Use the Taylor approximation to obtain the corresponding rotation tensor

18. Integrate constitutive relations

19. **end for**

20. **for** each child constituent **do**

21. Compute the divergence of the microstress vector

22. Compute the twin volume fraction

23. **end for**

24. Compute the local stiffness tensor

25. Compute the local plastic strain tensor

26. **end for**

27. Compute the convergence indicator $e_{eq}$

28. $j \leftarrow j + 1$

29. **end while**

30. Compute the macroscopic stress tensor

31. Compute the convergence indicator $e_{bc}$

32. **end while**

33. Save state variables

34. **for** Each voxel **do**

35. **for** Each child or parent constituent **do**

36. Update Euler angles

37. **end for**

38. **end for**

39. **end for**

**Algorithm 1** Computation of the effective behavior of a volume element. The tolerances on the boundary conditions and static equilibrium conditions are denoted by $e_{bc,c}$ and $e_{eq,c}$, respectively.

**Reviewer** Specifically: Equation (3) is not useful, just report (4), which is the one you used.

**Authors** Following the above recommendation, equation (4) has been removed from the revised version.
**Reviewer** The last sentence in Section 2.1 is misleading, you need to specify that a grain is constituted of several material points in your model.

**Authors** The last sentence of Section 2.1 has been modified as follows: "Also, while the stress state is uniform within the different constituents of a material point, the stress state is not identical for the parent and twinned domains of a given grain, which is composed of an infinite number of material points."

**Reviewer** Equation (6), is divergence of sigma, you inverted the symbols.

**Authors** Equation (6) (now (5)) is correct. The divergence of the stress tensor is a vector whose components are $(\sigma \cdot \nabla)_i = \sigma_{ij} \nabla_j$.

**Reviewer** Equation (9), there is slip in the twinned crystal, do you have evidences that this happens in your material? This is not always the case, you need to provide evidences.

**Authors** As discussed by Louca et al. (2021), 3D-XRD measurements indicate that the existence of easy slip systems within twins can lead to stress relaxation, which is why crystallographic slip is considered within twins. This aspect is now mentioned in the revised manuscript: "...The contribution of crystallographic slip to the plastic strain tensor of a child constituent is considered here to describe the potential stress relaxation associated with the possible existence of easy slip systems within twins (Louca et al., 2021)".

**Reviewer** Equations (16)-(17), where do the values of A, B, C come from?

**Authors** In Equations (16-17), now (15-16), no value is given for the A, B and C parameters. The values of these material parameters are provided in Table 1 (in the paper). The value of A, which can be interpreted as a barrier to twin nucleation, has been adjusted to reproduce the experimental stress-strain data. The values of the B and C parameters control the surface energy density for a twin/parent interface that is either parallel or perpendicular to the twin plane. In the present work, the value of B has been selected to obtain a coherent twin boundary energy density $g^{tw}$ of 0.13 J/m², which is the value given by Pei et al. (2017). Specifically, the twin boundary energy $g^{tw}$ is given by:

$$g^{tw} = \pi \sqrt{\frac{AB}{32}}$$

Also, the C parameter has been fixed to 10 × B, which favours twin propagation rather than twin expansion. To clarify these aspects, the following paragraph has been incorporated in Section 3.3 (Material parameters):

The values of the B and C parameters control the energy density for a twin/parent interface that is either parallel or perpendicular to the twin plane<fnref n="1" />. In the present work, the value of B has been selected to obtain a coherent twin boundary energy density of 0.13 J/m² (Pei et al., 2017). Also, the ratio C/B has been fixed to 10 to favour twin propagation over twin expansion.

**Reviewer** You need to specify at the beginning of the model Section if this is a small strain or finite strain framework.

**Authors** The adoption of the infinitesimal strain framework is now mentioned at the beginning of Section 2 (model description).

**Reviewer** Provide a citation for equation (26) and give a brief explanation of where it comes from.

**Authors** The reference to the standard textbook of Gurtin et al. has been incorporated together with the following sentence in the revised version of the manuscript: "The second law of thermodynamics allows defining the density of power being dissipated into heat. In a purely mechanical context, this quantity, which is denoted by d, is obtained from the difference between the power developed by internal forces and the evolution rate for free energy (Gurtin et al., 2010):

$$d = p_i - \dot{a}.$$

<footnote n="1">¹ The energy density is equal to $\pi \sqrt{AB/32}$ (respectively $\pi \sqrt{AC/32}$) for a twin/parent interface that is parallel (respectively perpendicular) to the twin plane.</footnote>
**Reviewer** You need to specify in Table 1 that the critical resolved shear stress values that you report are initial values, then you have hardening.

**Authors** This aspect is now mentioned in the caption of Table 1.

**Reviewer** Where do you take the parameters *K* and *N* from?

**Authors** As indicated in the original manuscript, these parameters have been adjusted to reproduce the experimental stress strain curves.

**Reviewer** Equations (32) and following: if you want to be precise, I think the dissipation potential should be added in equation (26) and then you can derive the power laws for slip and twinning. Otherwise it looks like this part is not connected with the previous theory. Since this is standard crystal plasticity, you could just report the power law equations with citations and say that it is the phenomenological way in which the system dynamically reaches equilibrium.

**Authors** In the authors opinion, it is not possible to add the dissipation potential in Equation (26). Indeed, the dissipative forces driving the evolution of internal variables should first be identified before constructing the dissipation potential, and deriving the evolution equations. As a result, following the framework of standard materials, it is more consistent to (i) establish the expression dissipation source *d*, (ii) determine the dissipative forces, (iii) construct a dissipation potential from the set of dissipative forces and (iv) derive the evolution equations.

**Reviewer** Which software package is used for the spectral solver? And for the constitutive model? In-house solutions can never be really trusted.

**Authors** Our in-house code is used for the simulations. As for any code development, different numerical tests (e.g. comparison with analytical and reference solutions) have been performed to ensure that the numerical implementation is correct.

**Reviewer** Section 3.1: provide element size and total size of the representative volume.

**Authors** The voxel size (there is no element in the context of the spectral method), which is about 1 µm and the volume element size, which is 150 µm, are now provided in Section 3.1 (Microstructure generation).

**Reviewer** Figure 3(b): it is confusing, much better to report the relative slip activity, meaning slip activity of the single mechanisms divided by total slip activity. You should plot relative slip activity as a function of the stress, so stress on the horizontal axis. In that way the curves would cover the plot evenly. Your way of plotting is really ugly.

**Authors** Figure 3(b) has been replotted with the axial strain on the horizontal axis and relative slip activity of the vertical axis. Same has been done for Figure 9(b).

**Reviewer** Equations (41)-(42): is that an average over a single grain? Or grains with same orientation in the polycrystal? This is because you mention it can be obtained by X-Ray diffraction, so it's not clear.

**Authors** Equations (41) and (42) refer to the average stress state over a single (parent or child) grain. To clarify this aspect, the corresponding paragraph has been modified as follows: "...In the recent years, 3D XRD techniques have been used to investigate the evolution of internal stresses within twinned grains. For a given twin-parent pair, these techniques provide some estimates of the average stress tensor for each single grain. In the following, the average stress tensor of a parent grain (respectively child grain) is denoted by $\bar{\sigma}^m$ (respectively $\bar{\sigma}^c$)."

**Reviewer** Figure 4: I doubt this model is not mesh independent, you must show a simulation with a smaller, maybe half, mesh size and discuss what is the difference you see in the twin phase field, I think nobody has discussed this mesh dependence for this type of phase field model yet.

**Authors** The model is mesh independent. While a local version of the model would be mesh dependent (because of the stress softening produced by twinning), the non-local contribution allows circumventing the difficulties associated with strain localization. Indeed, when the twin volume fraction excessively increases within a single material point, the non-local contribution allows "diffusion" in the sense that the growth of the twin domain in neighbouring material points
is favoured. As a result, when looking at the spatial distribution of the twin volume fraction (Fig 4a), one can see that there is no sharp transition between twinned and untwinned regions because of the non-local contribution to twin expansion. To clarify this aspect, the following sentence has been incorporated in Section 2.4 (Evolution equations): "...The term $A (1 - 2\phi_c)$ can be interpreted as the resistance to deformation twinning. The material parameter A therefore represents a barrier to the initiation of deformation twinning. Also, when the twin volume fraction increases, the resistance to deformation twinning decreases. The progression of twinning is therefore accompanied by a softening phenomenon that promotes plastic strain localization. While a local version of the model would be mesh dependent (because of the stress softening produced by twinning), the non-local contribution allows circumventing the difficulties associated with strain localization. Indeed, when the twin volume fraction excessively increases within a single material point, the non-local contribution favours the growth of the twin domain toward neighbouring material points."

When preparing numerical simulations, some preliminary simulations with lower resolutions of 64³, 96³, 128³ and 192³ voxels. As shown in Figure 1, the macroscopic behavior under uniaxial compression is correctly evaluated, even for a resolution of 64³. A minimum resolution of 128³ is however needed to estimate the grain averaged stress/strain tensors. For such a resolution, the average axial stress or strain within a parent grain at the end of a uniaxial compression test is similar to that obtained for a resolution of 192³ voxels, see Figure 2.

<figure><img src="image_2.png" /><figcaption>**Figure 1** Macroscopic stress-strain curves obtained for uniaxial compression with different resolutions.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2** Average axial stress versus average axial strain within a parent grain at the end of a uniaxial compression test with different resolutions.</figcaption></figure>

**Reviewer** What determines the spacing between twin bands? Are the twin bands mostly nucleating at triple junctions? Then, are the triple junctions determining the spacing between twins? In that case, is the grain size determining the spacing between twins in a grain? I had this doubt also about the Liu, Shanthraj, 2018 paper, and nobody has yet given a reasonable explanation.
**Authors** At this stage, we do not have any clear answer to this interesting question. From our numerical results, it seems that there is a connection between grain size and spacing between twin bands. However, our numerical data set is not extensive enough to draw sound conclusions. Specifically, one could imagine that twinning induced-stress relaxation also impacts spacing between twin bands. Indeed, both numerical and experimental indicate that twin formation results in a local stress relaxation. As a consequence, one may imagine that in the vicinity of a freshly nucleated twin, the stress field is largely relaxed, hence preventing twin nucleation in this region. Spacing between twins could therefore also depend on how important local stress relaxation is. Such a scenario is purely speculative at this stage, and further investigations need to be carried out.

**Reviewer** Figure 5: the orange data cover the violet data. This is not the proper way of plotting. You need to randomize the order in which circles are plotted, so that violet and orange data are equally represented. Or just split in two figures.

**Authors** Figure 5 has been modified according to the above suggestion to facilitate reading.

**Reviewer** In the conclusions: you claim that your model is able to reproduce the lenticular shape of twins, but your results do not show that and I judge unlikely that this model can reproduce the lenticular shape. Lenticular shape means at the twin tip there should be a cusp, but there is no mechanism in your model that can reproduce that. Based on equations (16)-(17) for the anisotropic surface free energy, I expect your twin tip to look like an oval, as the simulations by Liu, Shanthraj, 2018 suggest. Your mesh is just too coarse to resolve the twin tip. Either you run a simulation with very small mesh, maybe 20 to 30 elements along the thickness of a single twin and you discuss the shape you get, otherwise just don't make this kind of claim.

**Authors** What is meant in the conclusion is that the proposed model allows replicating the morphological texture evolution resulting from the concentration of plastic shear strains within twin domains. To avoid any confusion, the conclusion has been changed to: "...Also, according to the numerical results, the development of texture and internal stresses is consistent with the experimental observations of the literature. Specifically, the proposed model naturally replicates the morphological texture evolution resulting from the concentration of plastic shear strains within twin domains."

**Reviewer 2 (Anonymous)**

**Reviewer** The contribution is original and quite useful. The modeling of twin growth is convincing. The following remarks should be addressed:

Please explain better the implications of Equation (4) as it seems to be contradicted later in the model description.

**Authors** Equation (4) indicates that locally (i.e. for a material point) the stress state within the different parent and child constituents of a material point is the same. However, for a given grain, which is composed of an infinite number of material points, no such assumption is made and the internal stress field is allowed to develop. As a result, the average stress state is not identical for the parent and twinned domains of a given grain. To clarify this aspect, the following modification has been made to the revised manuscript: "...Also, while the stress state is uniform within the different constituents of a material point, the stress state is not identical for the parent and twinned domains of a given grain, which is composed of an infinite number of material points."

Also, to further indicate that the uniform stress state assumption is used consistently when describing the model, the equations defining the resolved shear stresses acting on the different slip systems (Section 2.3) are now presented as follows:

$$ \tau_c^\alpha = -\frac{\partial a_c}{\partial \gamma_c^\alpha} = \mathbf{m}_c^\alpha \cdot \boldsymbol{\sigma}_c \cdot \mathbf{n}_c^\alpha = \mathbf{m}_c^\alpha \cdot \boldsymbol{\sigma} \cdot \mathbf{n}_c^\alpha, $$

$$ \tau_m^\alpha = -\frac{\partial a_m}{\partial \gamma_m^\alpha} = \mathbf{m}_m^\alpha \cdot \boldsymbol{\sigma}_m \cdot \mathbf{n}_m^\alpha = \mathbf{m}_m^\alpha \cdot \boldsymbol{\sigma} \cdot \mathbf{n}_m^\alpha. $$

In a similar fashion, the density of power developed by internal forces (section 2.4) is now presented as follows:

$$ p_i = \sum_c \phi_c \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}}_c + \phi_m \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}}_m + \sum_c \dot{\phi}_c \boldsymbol{\sigma} : (\boldsymbol{\varepsilon}_c - \boldsymbol{\varepsilon}_m) + \sum_c \xi_c \dot{\phi}_c + \sum_c \boldsymbol{\eta}_c \cdot \nabla \dot{\phi}_c. $$
**Reviewer** The authors claim that they capture texture development. However they do not explain how they handle lattice rotations and the only pole figure shown is not very clear, nor convincing.

**Authors** In the present context, texture evolution is controlled by two different phenomena: (i) lattice rotation due to crystallographic slip and (ii) twinning. Though the contribution of crystallographic slip to texture evolution is very limited, it has however been considered in the present work. For this purpose, the initial orientation of a (child or parent) constituent is specified with a set of three Euler angles that can be updated with the Bunge's formula (Bunge, H. J., 1968, Phys. Status Solidi, 26, pp. 167-172) from the lattice spin tensor. The lattice spin tensor of a (child or parent) constituent is obtained from the difference between the total spin tensor $\dot{\omega}$, which is the skew-symmetric part of the displacement gradient tensor, and the plastic spin tensor, which is calculated from the plastic shear strain rates. The lattice spin tensor of a child or parent constituent ($\dot{\omega}_c^e$ or $\dot{\omega}_m^e$) is therefore calculated according to:

$$ \dot{\omega}_c^e = \dot{\omega} - \sum_{\alpha} \text{skw}(\mathbf{m}_c^{\alpha} \otimes \mathbf{n}_c^{\alpha}) \dot{\gamma}_c^{\alpha}, \quad (1) $$

$$ \dot{\omega}_m^e = \dot{\omega} - \sum_{\alpha} \text{skw}(\mathbf{m}_m^{\alpha} \otimes \mathbf{n}_m^{\alpha}) \dot{\gamma}_m^{\alpha}. \quad (2) $$

The contribution of twinning to texture evolution is naturally accounted for in the present framework. Indeed, a child constituent does not contribute to the initial crystallographic texture since the corresponding volume fraction is zero. When a twinning system is activated, the volume fraction increases and the contribution of a child constituent to the global crystallographic texture of the volume element progressively increases. The above method, which is used to handle texture evolution, is now described in Section 3.4.2 of the revised manuscript.

Also, to facilitate interpretation, the pole figure showing the impact of twinning has been modified to differentiate the contributions of parent constituents from that of child constituents. This clearly indicates that, as a result of twinning, basal poles tend toward the rolling direction at the end of the compression stage.

**Reviewer** The direction of the uniaxial loading (RD?) is not so clear and it could.

**Authors** The loading direction is RD, which is mentioned in Section 3.2 (Loading conditions). To clarify this aspect, the fact that the loading direction is RD is now mentioned in Sections 3.4.1 and 3.4.2: "*The experimental and numerical stress-strain curves obtained under uniaxial compression along the rolling direction are plotted in Figure...*" "*The numerical and experimental stress-strain curves obtained for the compression-tension loading path along the rolling direction are shown in Figure...*"

**Reviewer** The non local approach should predict twin thicknesses independent of the the mesh resolution. It would be worth demonstrating this. Also the influence of grain size on twinning is a subject of debate and the proposed model appears to address this issue. This too should be discussed.

**Authors** Since the proposed approach relies on a continuous description of twinning, there is no straight forward method to compute a quantitative indicator of the twin thickness. However, one can observe the spatial distribution of the twin volume fraction obtained at the end of a uniaxial compression test for different resolutions (128³ and 192³, see Figure 3). From a qualitative point of view, the spatial distributions look similar. Specifically, for the twin domains that can clearly be identified, the thickness is independent of the the mesh resolution. The proposed model does not allow considering the influence of grain size on twinning since the interactions of twins with grain boundaries are not considered. This limitation of the proposed approach is now mentioned in the final section of the paper: "...*The interactions of twins with grain boundaries, which have been ignored here, should also be included in constitutive relations to investigate the impact of grain size on the development of twins.*"

**Reviewer** This statement "A more realistic description of the behavior under tension would require considering kinematic hardening." sounds somewhat speculative. The authors should further develop and motivate their idea.

**Authors** When the loading direction is reversed, the isotropic/kinematic nature of the hardening impacts the flow behavior. Specifically, in contrast with the adopted isotropic hardening rule, a kinematic
<figure><img src="image_4.png" /><figcaption>**Figure 3** Spatial distribution of the twin volume fraction at the end of a compression test for different resolutions: 128³ (left) and 192³ (right).</figcaption></figure>

hardening rule would facilitate the development of plasticity upon load reversal. As a result, when switching from uniaxial compression to uniaxial tension, the adoption of a kinematic hardening rule would promote the development of plasticity, hence reducing the axial stress, which would be more consistent with the experimental results obtained after load reversal.

To clarify the aspect, the second paragraph of section 3.4.2 has been modified as follows: "A more realistic description of the behavior under tension would require considering kinematic hardening. Specifically, in contrast with the adopted isotropic hardening rule, a kinematic hardening rule would facilitate the development of plasticity upon load reversal. As a result, when switching from uniaxial compression to uniaxial tension, the consideration of kinematic hardening rule would promote the development of plasticity, hence reducing the axial stress, which would be more consistent with the experimental results obtained after load reversal."

**Reviewer 3 (Shailendra Joshi)**

**Reviewer** In this paper, the authors present a non-local model of twinning that adopts gradients of twin volume fractions in a thermodynamical framework. The model is implemented in a spectral solver and used to describe the uniaxial mechanical response of a polycrystalline Mg alloy. While the work is of interest to the mechanics community, there are several unresolved questions in the current treatment:

The model is non-local insofar as twinning is concerned. It is unclear as to why a similar treatment is not deemed necessary for dislocation slip. The finite-size effect of twins should, in principle, be of consequence in the dislocation plasticity as well for two reasons:

(i) Relation of grain size in relation to the twin size

(ii) Constraint effects on the dislocation plasticity (aka slip gradients) by finite twin lamellas. How do the authors rationalize ignoring the size effect in dislocation plasticity?

**Authors** We agree that a similar treatment should be adopted for crystallographic slip to consider (i) the impact of grain size on twin size and (ii) the contribution of twin boundaries to strain hardening. These aspects have been ignored here, mostly for simplicity reasons. A possible strategy to evaluate the impact of grain size would consist in adopting the strain gradient plasticity framework. However, this would require treating the spatial gradients of plastic shear strains as additional state variables, which would be computationally expensive. Also, to consider interactions between twin boundaries and dislocations, we have explored the possibility of including the impact of twin volume fraction gradients on the hardening rate. While there is no theoretical difficulty, such a theory involves additional material parameters that are difficult to identify with the present experimental data set.

To mention the above limitations, the conclusion of the proposed paper has been re-written as follows: "*Future work should focus on the role of twin boundaries on strain hardening. Indeed, while experimental studies indicate that twin boundaries contribute to strain hardening (Basinski et al., 1997), the impact of twin boundaries on the resistance to dislocation motion has been ignored in the present work. The coupling between slip and twinning deformation modes should be considered for future developments. For this purpose, the possibility of including the role of twin volume fraction gradients within the hardening rule will be explored. Indeed, twin boundaries correspond to the*
regions where such gradients are important. Twin volume fraction gradients should therefore affect the development of plasticity. Also, the interactions of twins with grain boundaries, which have been ignored here, should be included in constitutive relations to investigate the impact of grain size on the development of twins."

**Reviewer** The title of the paper indicates moderate strains. However, the simulations are carried out up to a strain of 0.5 % (i.e., 0.005). This is hardly a moderate strain level, more so a strain level that causes incipient macroscopic yield. Are the authors referring to local strains, e.g., Figure 4(c)? Given the strong texture, it would have been useful to see a result that leads to saturation of twinning thereby causing a rapid (S-shaped) hardening of the stress-strain response.

Notwithstanding, how do the authors justify using an additive decomposition of the strain if the local strains are high? Wouldn't finite strain effects be important locally?

**Authors** We agree that the term "moderate strain" may be confusing. To avoid confusion, the title of the proposed paper has been changed to "A non-local model for the description of twinning in polycrystalline materials in the context of infinitesimal strains: application to a magnesium alloy" to indicate that the present work relies on the infinitesimal strain theory. Also, to reach twinning saturation, the finite strain framework should be adopted since it would require prescribing very large strains to the volume element. Finally, as shown in Figure 4(c), the maximum equivalent strain is about 7 % so the infinitesimal strain assumption remains reasonable.

**Reviewer** In Equation (30), how are $\xi_c$ and $\zeta_c$ additive? One seems to be a work-conjugate to the rate of twin v.f. while the other is a work-conjugate to the twin v.f. Their dimensions are different. Am I missing something?

**Authors** While $\zeta_c$ is the free energy conjugate to the twin volume fraction, $\xi_c$ is the work conjugate to the twin volume fraction. Both quantities are stress-like quantities. Specifically, as indicated by Equation (5), $\xi_c$ is obtained from the divergence of the work conjugate to the twin volume fraction gradient (i.e. $\xi_c = \nabla \cdot \eta_c$).

**Reviewer** In the polycrystal simulations, how many grains are modeled? Is the grain sampling effect considered? I suspect, there should be statistical effects that may give a broad variability in the macroscopic and micromechanical characteristics. An exemplar of this is Figure 5. While the context related to local versus global Schmid factors is reasonable, the figure itself does not seem to particularly reveal the stated context.

**Authors** As indicated in Section 3.1 (Microstructure generation), the volume element is composed of 50 grains. The choice of the number of grains is guided by a compromise between representativity and sufficient spatial resolution for the application of the spectral method. Preliminary investigations have shown that the relative difference between macroscopic axial stress obtained with different volume elements with 50 grains each is about 10%.

Regarding Figure 5, a more quantitative approach of interpreting the local/global Schmid factors consists of considering the 52 twin variants with the most important activity (i.e. with a volume fraction higher than 1 % at the end of the compression test). Table 1 shows the average Schmid factor and the corresponding standard deviation that have been calculated for these most active variants. According to the results, the standard deviation is much higher for the global Schmid factor, which further confirm that the propensity of a given twin variant to be activated is better described with the local Schmid factor than the global one. These quantitative results have been included in Section 3.4.1 (uniaxial compression) of the revised manuscript.

<table><caption>**Table 1** Average Schmid factor and corresponding standard deviation calculated at the end of the compression test for the most active twin variants (i.e. with a volume fraction higher than 1%) with either the global or local definition of the Schmid factor.</caption><thead><tr><th></th><th>Local</th><th>Global</th></tr></thead><tbody><tr><td>Average</td><td>0.36</td><td>0.38</td></tr><tr><td>Standard deviation</td><td>0.037</td><td>0.064</td></tr></tbody></table>
**Reviewer** Overall, the assessment of the model seems a bit narrow. A few concerns regarding this aspect: It is not clear why the authors chose the particular problems. What if the polycrystal were to be compressed in the rolling plane but perpendicular to the rolling direction. It should still activate twinning, but would the features (e.g., variants) be different. What would their impact be on the micromechanical and macroscopic behaviors?

**Authors** The loading conditions have been selected to reproduce the experimental conditions of Louca et al. (2021). While such loading conditions do not prove the strict validity of the proposed model, they have been selected to see whether the numerical results are consistent with the experimental observations of Louca et al. (2021). It is also worth mentioning that the proposed model is not restricted to any particular type of loading conditions. To clarify this aspect, the following sentence has been incorporated in section 3.2 (Loading conditions): "While the selected loading conditions do not allow a full assessment of the proposed model, the corresponding numerical results can be compared to the experimental data set of Louca et al. (2021) to see whether the description of twinning is consistent with experimental observations or not."

**Reviewer** How critical is the non-local effect to the model predictions? If the non-local term were to be suppressed and the model parameters refitted to the experimental data, how different would the corresponding predictions be?

**Authors** The non-local effect is not very critical as twinning is mostly driven by the mechanical contribution to the driving force. In the present case, non-locality allows circumventing the difficulties associated with mesh dependency. Indeed, the local contribution to the driving force for twinning contains a softening contribution (i.e. $A(1 - 2\phi_c)$). If the model was purely local, the twin volume fraction could excessively increase within a single material point, which would lead to strain localization. The non-local contribution allows "diffusion" in the sense that the growth of the twin domain in neighbouring material points is favoured.

To better explain the motivations for non-locality, the following sentence has been incorporated in Section 2.4 (Evolution equations): "...The term $A(1 - 2\phi_c)$ can be interpreted as the resistance to deformation twinning. The material parameter A therefore represents a barrier to the initiation of deformation twinning. Also, when the twin volume fraction increases, the resistance to deformation twinning decreases. The progression of twinning is therefore accompanied by a softening phenomenon that promotes plastic strain localization. While a local version of the model would be mesh dependent (because of the stress softening produced by twinning), the non-local contribution allows circumventing the difficulties associated with strain localization. Indeed, when the twin volume fraction excessively increases within a single material point, the non-local contribution favours the growth of the twin domain toward neighbouring material points."

**Reviewer** p.11: The statement "For the investigated loading...and twinning" is vague. What kind of competition is being referred to here? Does twinning inhibit slip and vice-versa?

**Authors** We agree that the above statement is vague, it should be more specific. What is meant here is that both crystallographic slip and twinning contribute to plastic deformation for the investigated loading conditions. To clarify this aspect, the conclusion has been modified as follows: "The proposed set of constitutive relations has been implemented in a spectral solver to model the behavior of a AZ31 magnesium alloy. For the investigated loading conditions, both crystallographic slip and twinning contribute to macroscopic plastic deformation. Specifically, the initial yielding during uniaxial compression coincides with the activation of basal slip systems while the transition toward the plastic regime is attributed to twinning."

**Reviewer** Figure 3: While it is a matter of choice, the compressive stress-strain curve shown in the negative quadrant is a bit disorienting. I suppose the reason for this stems from the cyclic loading investigated in the paper. Perhaps the authors want to be consistent in the figures. Just an observation. For clarity, there should be a legend denoting which data are from the experiments (dots) and which one is the simulation result (lines).

**Authors** Indeed, Figure 3 (uniaxial compression) was plotted in the negative quadrant to be consistent with Figure 9 (compression + tension). In the revised version of the manuscript, Figure 3 is
plotted in the positive quadrant. The captions of both figures have been modified to differentiate between experimental and numerical data.

**Reviewer** The paper seems to have missed out on referring to some significant contributions (in a broader context), particularly by Joshi and colleagues (e.g., Zhang and Joshi (2012) 60, JMPS, 945-972, Indurkar et al (2020), 132, IJP, 102762, Ravaji et al (2020), 208, Acta Mater, 1167743).

**Authors** Following the above suggestion, the above references, that deal with local crystal plasticity models, have been included in the introduction section.

## Review of version 2

Permalink: arXiv.org:2012.13156v4

**Authors** Dear Editor, we would like to thank the reviewers for their help in improving the manuscript. The minor revisions suggested by the reviewers have been addressed in the revised version. A detailed response is provided below.

### Reviewer 1 (Nicolo GRILLI)

**Reviewer** Dear Authors, I am happy with your modifications, just in the answer to question about the mesh size dependence of the model you report two figures, please introduce them in the manuscript as well with the descriptions. You don't want to convince only me but also future readers.

**Authors** Following the above suggestion, an appendix that includes the results regarding mesh dependence has been incorporated in the revised version of the manuscript.

**Reviewer** After that modification, I am happy for this paper to be published.

### Reviewer 2 (Anonymous)

**Reviewer** The paper is nicely written. Handling the formation of twin lamellae with a non-local crystal plasticity model is useful and the results presented here are impressive.

The infinitesimal strain assumption is awkward and not really justified in view of the proposed model. It could be applied only to the elastic part of the deformation. Why not adopting an updated Lagrangian framework accounting for progressive lattice rotations? Why not using strain rates, instead of total strain, in Section 2.3?

**Authors** The infinitesimal strain assumption seems adapted to the present work, which focuses on the early stages of plasticity. As shown in Figure 4(c), strains remain small, the maximum equivalent strain being inferior to 7%. It is worth mentioning that the infinitesimal strain assumption has been used in many other studies that deal with deformation twinning, for instance:

* Brown, D., S. Agnew, M. Bourke, T. Holden, S. Vogel, and C. Tomé (2005). Internal strain and texture evolution during deformation twinning in magnesium. Materials Science and Engineering: A 399.1. Measurement and Interpretation of Internal/Residual Stresses, pp. 1-12. 10.1016/j.msea.2005.02.016

* Cherkaoui, M. (2003). Constitutive equations for twinning and slip in low-stacking-fault-energy metals: a crystal plasticity-type model for moderate strains. Philosophical Magazine 83.31-34, pp. 3945-3958. 10.1080/14786430310001603355

* Clausen, B., C. Tomé, D. Brown, and S. Agnew (2008). Reorientation and stress relaxation due to twinning: Modeling and experimental characterization for Mg. Acta Materialia 56.11, pp. 2456-2468. 10.1016/j.actamat.2008.01.057

* Paramatmuni, C. and A. K. Kanjarla (2019). "A crystal plasticity FFT based study of deformation twinning, anisotropy and micromechanics in HCP materials: Application to AZ31 alloy". International Journal of Plasticity 113, pp. 269-290. 10.1016/j.ijplas.2018.10.007

* Xu, F., R. Holt, M. R. Daymond, R. Rogge, and E. Oliver (2008). "Development of internal strains in textured Zircaloy-2 during uni-axial deformation". Materials Science and Engineering: A 488.1, pp. 172-185. 10.1016/j.msea.2007.11.018

Formulations using either the updated Lagrangian framework or strain rates are hypoelastic, i.e. constitutive equations are presented in a rate form only. Such formulations are generally
not obtained within a consistent thermodynamic framework, i.e. from both a state potential and a dissipation potential. Also, the development of constitutive relations relies on the Reuss assumption to partition the strain tensor between child and parent constituents. While the adoption of this assumption is straightforward in the context of infinitesimal transformations, additional difficulties are involved for finite strains. Specifically, for a material point, the problem of determining the deformation gradient tensor associated with each constituent from the Piola stress tensor does not necessarily has a unique solution.

Though it is possible to extend the proposed set of constitutive relations to the context of finite strains, additional efforts are needed to include such equations in a consistent framework.

**Reviewer** Consideration of lattice spin comes very late, in Equations (42)-(43). Also there is a mistake: the spin is the skew symmetric part of the velocity gradient, and not of the deformation gradient tensor. With twinning, comes of course the question of the twin-parent orientation relationship as plastic strain proceeds, as discussed in several publications about texture development influenced by twinning.

**Authors** We agree, the consideration of lattice rotation comes late. In the revised version of the manuscript, the strategy used to handle lattice rotations is described in Section 2, which is dedicated to model description. As a result, all the equations used for the present approach are presented in the same section, which is more consistent. Also, the spin tensor is indeed the skew symmetric part of the velocity gradient, as indicated by Equations (40)-(41), not the displacement gradient. This mistake has been corrected in the revised version. Finally, as mentioned by the Reviewer, the strategy used to handle lattice rotations considers that the child-parent orientation relationship is not impacted by plastic strains, which would be inappropriate in the context of finite strains. This is now mentioned in the revised version of the manuscript (Section 2.5 Texture evolution).

**Reviewer** The question of kinematic hardening is an important one too. I think that it would be wrong to consider that the present model does not include any kinematic hardening, as there is stress heterogeneity between twin lamellae and their surrounding. It would be interesting to pay more attention to the detwinning process and also to evaluate the elastic-plastic transient upon load reversal when detwinning is not allowed.

**Authors** As mentioned in the original submission, the macroscopic behavior is impacted by kinematic hardening. Indeed, the internal stresses resulting from plastic strain incompatibilities provide some contribution to kinematic hardening. However, according to the numerical results, kinematic hardening remains underestimated since the behavior upon load reversal is not fully consistent with experimental results. Specifically, internal stresses allow reactivating basal slip systems when switching from compression to tension, but are not high enough to reactivate prismatic slip systems. To improve the description of the plastic behavior, one could therefore include kinematic hardening at the microscale, i.e. *via* constitutive relations by incorporating an additional contribution to the yield function.

To clarify this aspect, the following comment has been incorporated in Section 3.4.2 (Uniaxial compression and tension) of the revised manuscript: "*While the internal stresses resulting from plastic strain incompatibilities contribute to kinematic hardening at the macroscopic scale, kinematic hardening is not considered at the microscale via constitutive relations. A more realistic description of the behavior under tension would therefore require considering this additional contribution.*"

**Reviewer** The literature review is extended, but somewhat over-focused on contributions from a few (well-known) teams.

**Authors** While there is a profuse literature dedicated to the modelling of deformation twinning, the proposed manuscript is not a review paper. In this perspective, we did our best to include the references (over fifty) that allow (i) underlying the original aspects of the proposed work, (ii) constructing constitutive relations and (iii) indicating the origin of experimental data.

Reviewer 3 (Shailendra Joshi)

**Reviewer** The revised paper does a reasonable job of addressing the concerns of this reviewer.
**Editor's assessment (Laurence BRASSART)**

The paper proposes a novel crystal plasticity-based constitutive model to describe twin growth in polycrystalline materials. The main novelty of the paper is the introduction of non-locality by considering twin volume fractions and their gradients as external variables. The model is implemented in a spectral solver and used to simulate twinning in Magnesium alloys. All three reviewers acknowledged the scientific quality and originality of the work. Major revisions were however requested to clarify a number of technical points. After the first revision, Reviewers 1 and 3 recommended the paper for publication. Reviewer 2 asked for further improvements, which then led the paper to being accepted.

**Open Access** This review is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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