An effective anisotropic visco-plastic model dedicated to high contrast ductile laminated microstructures

In particular types of layer-or lamellar-like microstructures such as pearlite and lath martensite


Introduction
Lamellar structures are known to be an essential building block of various materials.At different scales, they are either formed to minimize the total elastic energy, e.g., in systems with structural variants such as pearlite (Queiros et al., 2019;Vermeij and Hoefnagels, 2022), intermetallic Ti-Al alloys (Yamaguchi et al., 2000;Huang et al., 2018), nano-bainitic structures (Bhadeshia, 2013;Caballero et al., 2019), and lath martensite (Morito et al., 2003), or manufactured to improve mechanical properties, e.g., in thermally-sprayed coating materials (Santana et al., 2008;Klusemann and Svendsen, 2012), in ultra-strong CuNb laminates (Wang et al., 2012), and semi-crystalline polymers (Flory and Yoon, 1978;Bartczak et al., 1992).Highly anisotropic behavior, originating from the different types of crystallographic or morphological heterogeneities of the phases within the laminate, is a typical characteristic of these structures.It has been observed in Ti-Al alloys that the lamellar interfaces can act as obstacles to the homogeneous plastic deformation, resulting in channeling of plasticity and hence an easy deformation mode parallel to the interfaces (Paidar and Yamaguchi, 2007).The molecular chain axes (chain direction) are constrained slip directions in semi-crystalline polymers, forcing extra kinematical constraints on the crystalline phase, resulting in extreme anisotropic behavior (Bartczak et al., 1992).
A prime example of inhomogeneous plastic activity, and the main motivation for the model developed in this paper, is the lamellar https://doi.org/10.1016/j.ijsolstr.2024.112757Received 20 April 2023; Received in revised form 30 December 2023; Accepted 5 March 2024 microstructure of lath martensite (Du et al., 2018).There is increasing experimental evidence that large plastic strains can occur in lath martensite (Du et al., 2016;Kwak et al., 2016), even though it is conventionally regarded as a brittle phase.It has been hypothesized that this plastic deformation, which occurs more favorably parallel to the lath boundaries (Du et al., 2019), is facilitated by thin retained austenite films located between the martensite laths (Maresca et al., 2014a).These films, by virtue of the orientation relationships between the martensite laths and their parent austenite, have slip systems that are approximately parallel to the lath, that can glide easily (Vermeij et al., 2023).In the limit of nanoscale films, plasticity can also be achieved in absence of dislocations within the austenite films, by means of phase transformation carried by interface dislocations (Maresca et al., 2018).Alternatively, the remarkable ductility in particular directions has been attributed to the morphology of the laths, which typically have an aspect ratio of 30 ∶ 1 in the length direction (habit plane direction) compared to the thickness direction, and hence the dislocations on the slip systems parallel to the habit plane can glide comparatively freely (Morsdorf et al., 2016;Ryou et al., 2020).Regardless of the exact physical mechanism underlying this phenomenon, it has been shown experimentally that the plastic activity/slip is mediated favorably in directions parallel to the lath boundaries (Du et al., 2019;Rezazadeh et al., 2022), which results in a significant degree of anisotropy in the plastic response.
It is often unclear how the local anisotropic of these laminate microstructures affects the mesoscale behavior, where multiple grains with different laminate orientations interact with each other.However, mesoscale models, containing multiple grains, are at a much larger length scale compared with the characteristic scale of the lamellar microstructure, and as such a full model that comprises multiple grains and resolves, at the same time, the lamellar microstructural features would be computationally prohibitive.One way to obtain the effective laminate response without resolving the laminate in its full detail is (computational) homogenization (Kouznetsova et al., 2004;Geers et al., 2010).Homogenization of laminates generally involves incorporating the interaction between two (or more) phases (Livingston and Chalmers, 1957;Kocks and Canova, 1981).Neglecting the interactions can result in too stiff (Taylor-type (Taylor, 1938;Hutchinson, 1976)), or too compliant (Sachs-type (Sachs, 1928)) behavior.To satisfy compatibility and equilibrium across the interfaces of phases/grains, intermediate models such as relaxed constrained type models or selfconsistent models have been developed (Lebensohn and Tomé, 1993;Molinari et al., 1987Molinari et al., , 1997)).Similar concepts have also been used to model the interaction between pairs of grains in non-lamellar microstructures such as polycrystals, e.g. in the LAMEL and ALAMEL models (Van Houtte et al., 2005).In these models, the representative volume is assumed to be a single crystal/grain/inclusion in a homogeneous equivalent medium.Lee et al. (1993) proposed a visco-plastic two-phase composite inclusion model for deformation and texture evolution in semi-crystalline polymers (with a lamellar structure of the crystalline and amorphous phases).Later they adopted the idea of using a two-phase composite inclusion for FCC polycrystals, where the composite inclusion is represented by a bicrystal with arbitrary misorientation and with an assumed extended planar interface (Lee et al., 2001).van Dommelen et al. (2003) extended this model to an elasto-viscoplastic formulation.The idea of two-phase laminates is used by Ortiz et al. (2000) for describing the lamellar dislocation structures which develop at large strains.Microstructures of lamellar type have long been treated within the context of the crystallographic theory of martensitic transformation (Pedregal, 1993;Bhattacharya, 1991;Stupkiewicz and Petryk, 2006).Based on the concept of a laminate composed of a martensite plate and an austenite layer, Kouznetsova and Geers (2008) modeled the transformation plasticity of variant formation in the martensitic transformation.Klusemann and Svendsen (2012) proposed a homogenization method in a small strain setting to capture the material behavior of two-phase laminates characterized by a thin-layer-type microstructure found in thermally-sprayed coating materials like WC/Fe.Finally, Maresca et al. (2016b) proposed a laminate-based reduced crystal plasticity model to account for the limited number of active slip systems in the austenite layers trapped between the laths of martensite, and validated the model using lath martensite microstructures.
Here, we focus on the plasticity of two-phase laminate-type microstructures, and particularly on the cases where one of the phases is sufficiently thin for it to be considered as a discrete plane, cf. the boundary sliding mechanism in lath martensite discussed above.Taking benefit of this feature, we model the system as a comparatively hard, isotropic elasto-viscoplastic matrix in which a family of parallel soft (implicitly) discrete plastic planes are embedded.The planes exhibit viscoplastic sliding which is assumed to be isotropic along the plane and which is driven by the resolved shear stress vector.Within their planes, these films experience the same stress as the matrix.However, they respond to it differently, depending on their properties and their orientation with respect to the loading direction.This additional viscoplastic deformation mode complements the total plastic deformation rate of the matrix, resulting in an effective model which has a much lower complexity and computational cost than the homogenized approaches discussed above, e.g., that by Maresca et al. (2016b) for lath martensite.
The paper is organized as follows.In Section 2, we formulate the effective elasto-viscoplastic laminate model within a finite deformation setting.In Section 3, the response of the model is assessed via a comparison to that of an infinite periodic two-phase laminate under shear and tensile loading.The anisotropic yielding behavior of the model is further examined in Section 4 for arbitrary loading directions, by comparing the yield surface predicted by the reduced laminate model to that of the infinite periodic reference model.The influence of the film spacing and of hardening are also investigated.Section 5 presents the application to lath martensite containing inter-lath retained austenite films.Predictions made by the effective laminate model for an infinite periodic microstructure are compared with those of a reference model which takes into account the crystallographic orientations of the martensite laths and the austenite films.Finally, a two-dimensional representative volume element of a dual-phase (DP) steel microstructure with a 40% volume fraction of martensite is simulated, and the results are compared with those of a standard isotropic plasticity model and a standard crystal plasticity model which do not account for the presence of the austenite films.For this problem, the gain in terms of computational effort of the model is discussed by comparing with the isotropic and crystal plasticity models, as well as the reduced lath martensite model proposed by Maresca et al. (2014a).
Throughout the paper, the notations , , , and D denote scalars, vectors, second-order tensors and fourth-order tensors, respectively.Single and double contractions are denoted by ⋅ and ∶, respectively.The inner product between two second-order tensors is given by  =  ⋅  as   =     and the double inner product of a fourth-order tensor with a second-order tensor, written as  = A ∶ , is defined according to   =     .The action ‖ ⋅ ‖ F designates the Frobenius norm, i.e.

√
tr( ⋅  T ), where tr( ) is the trace operator, and the superscript ⋅ T indicates transposition.The tensor (or dyadic) product between two vectors is denoted by  =  ⊗ , where   =     .Finally, ̄⋅ denotes that a quantity is applied on average.

Model formulation
Fig. 1(a) shows a sketch of an infinite periodic two-phase laminate obtained by embedding infinitesimally thin soft films in a harder, continuous matrix.The elastic response of the model is governed by the homogeneous matrix phase, whereas plasticity is composed of two separate contributions: an isotropic visco-plastic part representing the matrix behavior, enhanced by a planar isotropic sliding mechanism that models the embedded films.The planar model captures a specific orientation-dependent plastic deformation mode that the composite may exhibit due to the presence of the softer thin films.In the case of shear exerted parallel to these films, see Fig. 1(b), they are expected to accommodate virtually all of the plastic strain, while applying stress in any other direction involves yielding of the harder matrix, resulting in a harder overall response.For instance, in the particular case of tension applied perpendicular to the films, Fig. 1(c), the films are inactive and the response is governed solely by the standard isotropic (matrix) part of the model.

Elasticity
The deformation of the effective laminate model is characterized by the deformation gradient tensor, , which can be split multiplicatively into its elastic contribution,  e , and a plastic contribution,  p , as follows, where  p represents an isochoric, orientation preserving plastic deformation.The elastic deformation and rigid-body rotation are encoded in  e .The stress induced by the elastic deformation is expressed via a Saint Venant-Kirchhoff type model as: where C denotes the fourth-order isotropic elasticity tensor and   = ( e − )∕2 is the elastic Green-Lagrange strain tensor, in which  e =  T e ⋅  e is the elastic right Cauchy-Green deformation tensor. e is the push-forward of the second Piola-Kirchhoff stress tensor  to the intermediate (plastic) configuration, i.e.  e =  p ⋅  ⋅  T p .The elastic response of the model is entirely governed by the matrix, and the films contribute only to the plastic part of the deformation.Accordingly, C represents the elastic stiffness of the matrix.

Plasticity
The effective laminate model incorporates, implicitly, extra discrete planar modes into the plastic response of the matrix.For this purpose, the plastic velocity gradient in the model is additively split into two separate contributions, in which  m p represents the contribution to the plastic velocity gradient by the matrix and  f p is the sliding-induced contribution to the plastic velocity gradient.

Visco-plasticity of the matrix
According to the standard isotropic visco-plasticity formulation with isotropic hardening model proposed in Anand (1985), the associated plastic velocity gradient for the matrix phase,  m P , is given in terms of the Mandel stress tensor,  =  e ⋅  e as where  dev is the deviatoric part of  and  eq = ‖ dev ‖ F ∕ √ 2 the von Mises equivalent stress.Note that we are using the shear-equivalent definition of the latter, for consistency with the planar plasticity in the film which is introduced in the next section.The plastic shearing rate ̇m is based on the widely adopted power-law relation (Hutchinson, 1976) where ̇m 0 is a reference plastic strain rate,  m y the current flow resistance (yield stress) of the matrix and  its strain rate sensitivity exponent.The evolution of the flow resistance, from the initial value,  m 0 , to the saturation value,  m ∞ , follows a phenomenological isotropic hardening law (Brown et al., 1989;Bronkhorst et al., 1992) with ℎ m 0 a modulus which characterizes the initial hardening, and  the hardening shape factor.The evolution of the hardening in the matrix, given by Eq. ( 6), indicates that only self-hardening of the matrix is accounted for.Latent hardening between matrix and film could in principle be included here in cases where such interactions are expected to be relevant.

Visco-plasticity of the films
Given the fact that the soft films are considered to be infinitesimally thin, we formulate their response as a planar (true) sliding, i.e. a relative displacement along a set of equidistant discrete planes which represent the films.This sliding is assumed to be driven solely by the shear stress acting along the planes and to follow the direction of the shear traction vector.Equilibrium between the films and the surrounding matrix material implies that this shear traction equals that in the matrix.Let  0 be the normal to the embedded sliding planes in the reference configuration.The shear stress  acting on these planes in the orientation preserving intermediate configuration can be computed from the full Mandel stress tensor, , by projecting it onto the planes as follows, in which  is the second-order identity tensor and | ⋅ | is the Euclidean norm.The sliding direction along the planes is given by The amount of sliding along an individual plane is characterized by its relative displacement .Its evolution is defined in rate form via a power-law relation similar to Eq. ( 5): where  f y is the flow resistance of the planar mode and v0 is a reference sliding velocity.Note that for simplicity we have assumed the same rate sensitivity exponent  as for the matrix.The flow resistance has an initial value  f 0 and evolves according to (cf. ( 6)) in which the constant  0 characterizes the initial hardening of the sliding mechanism, and  f ∞ is a saturation value.The exponent  is once more taken equal to that of the matrix for simplicity.No latent hardening due to activity of the matrix is considered.
The contribution of the sliding mechanism is accounted for in (3) in a homogenized sense via with  being the spacing of the films.This expression is reminiscent of the contribution of a single slip system in crystal plasticity, with the quotient v∕ of the relative sliding velocity per plane and the spacing of the planes as the mean (homogenized) plastic shear rate parallel to the planes.An important difference with crystal plasticity, however, is that the direction of shearing,  0 , is not pre-determined by the crystallography of the grain, but follows the maximum shear traction -similar to isotropic plasticity, but in the present case restricted to directions parallel to the sliding plane.Accordingly, the plastic response of the model is isotropic in the plane of the films, but depends on the angle between the plane and the applied loading.
It is furthermore worth mentioning that the two contributions to the overall plastic distortion rate, cf.Eq. ( 3), are ultimately driven by the same Mandel stress tensor  (of the matrix): the isotropic contribution by the equivalent shear stress  eq , via (4), ( 5), and the planar isotropic part due to film sliding by the resolved shear stress  via ( 9), ( 11).This also implies that the numerical implementation of the model quite closely follows that of standard visco-plasticity or crystal plasticity.

Reference model
The formulation of the model given in Section 2 leads to an overall anisotropic response, and the degree of anisotropy depends on the degree of activation of the planar sliding mode.Here, we aim at characterizing the extent of inherited anisotropy and assessing the impact of the assumptions made with respect to a reference model in which the soft films are modeled explicitly and with a finite (but small) thickness.The reference model consists in an infinite, periodic laminate of two alternating phases: comparatively thin, soft films embedded in a harder matrix.It is subjected to a homogeneous far-field applied deformation and/or loading.Taking benefit of the periodicity of the problem and the homogeneity of the individual phases, we may restrict our analysis to a periodic cell as sketched in Fig. 2 and adopt spatially constant deformation gradients  m and  f in the matrix and films, respectively, as well as spatially piecewise constant stresses.Note that, even though the sketch of Fig. 2 is two-dimensional, all components of the threedimensional deformation and stress tensors are taken into account.The constitutive law for both phases of the laminate is standard isotropic visco-plasticity, and thus their behavior separately follows Eqs. ( 1)-( 6) without the second term in (3), i.e. with  f p = .The effective, mesoscale deformation gradient, F, and the first Piola-Kirchhoff stress, P, are related to those in the two phases by volume averaging: where  represents the volume fraction of the soft films, which we consider to be small.Eqs. ( 12) and ( 13) ensure satisfaction of the Hill-Mandel condition, stating the equivalence between the virtual internal work on the microscopic and macroscopic scale, i.e. ⟨ ∶ ⟩ = P ∶  F, with ⟨⋅⟩ indicating the volume average over the periodic cell (Hill, 1963;Mandel, 1966).To obtain a closed system, they must be complemented by equations which ensure compatibility and equilibrium across the interface between matrix and films, as follows.Let  0 be the normal to the interface between the two phases.Kinematic compatibility requires the in-plane deformation gradient to be continuous across the interface, Furthermore, equilibrium requires traction continuity across the interface, Note that the kinematic compatibility and equilibrium conditions ( 14), ( 15) are necessary in the reference model, which has a finite film thickness, but are not explicitly enforced in the effective laminate model of Section 2, in which the films have no thickness.Accordingly, they have no resistance against in-plane deformation and experience the same out-of-plane stress components as the matrix.
To be able to compare the responses of the reference model and the laminate model, the film thickness  in the reference model must be taken into account.The sliding velocity v in Eq. ( 9) may be related to the shear strain rate ̇f in the reference model via Accordingly, the parameters v0 and  0 are related to ̇f 0 and ℎ f 0 via v0 =  ̇f 0 and  0 = ℎ f 0 ∕().The volume fraction of the films is taken to be  = ∕ = 0.05.The material parameters of the two phases in the  1.The mechanical phase contrast between the matrix and the films is  m 0 ∕ f 0 = 2.For simplicity, ideal plasticity is assumed, by taking ℎ m 0 = ℎ f 0 = 0; this renders the values of  ∞ and  irrelevant.

Comparison of shear and tensile responses
The responses of the models are now investigated under two applied load cases; (i) simple shear parallel to the films, and (ii) uniaxial tension applied perpendicular to the films.Tension applied parallel or perpendicular to the films is expected to result in an equal overall response for the reduced laminate model, as both load cases impose no shear on the films, and the response is hence fully determined by the matrix.In the shear load case a shear strain rate of 0.01s −1 was applied, and in the tensile load case an axial strain rate of 0.01s −1 .We represent the overall response of the models, here and in what follows, by the shear-equivalent von Mises true stress and strain, i.e. τeq = ‖‖ F ∕ √ 2 and γeq = √ 2 ‖ ln ‖ F , where  and  are the Cauchy stress tensor and left stretch tensor, respectively.Fig. 3(a) compares the overall stress-strain response of the reference model, plotted as solid lines with square markers, and of the effective laminate model, plotted as solid red lines, for the two load cases.For both load cases, the response of the laminate model follows that of the reference model with an error of less than 2%.It is shown that both models are softer by a factor of two in simple shear compared to the tension case.As the shear is applied parallel to the soft films, all the plasticity in the system is carried by the films (remember we are considering ideal plasticity here) and the matrix remains elastic.Therefore, the system yields at approximately the yield strength of the films, i.e.  f 0 = 200 MPa, and stays at the same stress level due to the absence of hardening.In tension perpendicular to the films (as well as parallel to the films -not shown here), on the contrary, no shear stress acts on the films, which are therefore not activated.Hence, the properties of the matrix determine the behavior of the system in this direction, and an equivalent yield stress of approximately  m 0 = 400 MPa is observed.
The above explanation is supported by Fig. 3(b) and (c), which compare the contributions of the two plastic modes in the laminate model and reference model under the applied shear and tension loading, respectively.It is shown that the plastic activity of the matrix in both models practically vanishes in simple shear, and hence all deformation is carried by the sliding mechanism, see Fig. 3(b).In tension perpendicular to the films, Fig. 3(c), the applied deformation is absorbed virtually completely by the matrix.In the effective laminate model, the films are inactive because they do not experience any shear stress.In the reference model, in which the films are modeled by isotropic plasticity and they have a lower yield stress than the matrix, they do activate -even before the matrix.Compatibility however requires that the matrix deforms to approximately the same degree.As soon as it also yields, the two phases hence start to contribute according to their volume fraction, implying that the matrix is dominant.

Influence of the film orientation
Next, we investigate in more detail the influence of the film orientation with respect to the applied deformation.The yield limit τ0 of the effective laminate model is defined as that macroscopic equivalent stress τeq at which where ̇c , is arbitrarily set to 0.002 s −1 .The advantage of this criterion over conventional criteria is that it is insensitive to the specific material parameters, such as the reference shear rate or hardening parameters.
It is important to note that the precise definition of the yield limit and the value of ̇c used do not affect our conclusions.For the reference model, the corresponding criterion reads Here one should notice the factor (1 − ) in front of the matrix plastic strain rate, as it is relevant for the reference model but not for the proposed effective laminate model, in which the films have no thickness.
To explore the anisotropy of the effective yield surface, we apply uniaxial tension at an arbitrary angle  with respect to the normal of the films,  0 .For the particular case in which the tensile direction coincides with the -direction (see Fig. 2), i.e. parallel to the films and hence corresponding with  = ∕2, the following overall rate of deformation gradient and first Piola-Kirchhoff stress tensor are applied: where λ = 0.01 s −1 is the stretch rate, and * indicates that the particular component of the tensor is free to evolve.For other loading directions, these conditions are rotated according to the loading angle  considered, whereas the orientation of the microstructure is kept fixed.The polar plot of Fig. 4 shows the yield surfaces obtained for the loading condition and the yield criteria discussed above for the effective laminate model and the reference model.In the diagram, the radius indicates the yield limit τ0 in MPa, and the angular coordinate indicates the angle between the normal to the films and the applied tensile loading direction, .The yield surface of the effective laminate model is shown as a red solid line and that of the reference model as a black solid line overlaid with square markers.The initial yield stresses of the individual phases in the laminate, the matrix and films, are shown as dashed lines, for comparison.The stress at which the matrix phase starts to yield,  m 0 = 400 MPa, is represented by the purple dashed circle.The green dashed curves represent the points at which the resolved shear stress acting on the films equals the initial sliding resistance of the films,  f 0 ; it is given by τ0 = 2 f 0 ∕( It is observed, first of all, that the yield surface of the effective laminate model follows that of the reference model almost perfectly.Both models show an extremely anisotropic response.The yield stress is maximum, τ0 ≈ 400 MPa, at  = 0, and  = ∕2.This is, as reported above, due to the fact that tension along or perpendicular to the films does not induce any shear stress on the soft films.Therefore, in the effective laminate model the films are inactive and yielding is governed solely by the isotropic matrix, as evidenced by the coincidence with the purple dashed circle in these regions.The reference model has a marginally lower yields stress due to the activation of plasticity in the (isotropic) films -cf.Fig. 3 and its discussion.The yield surface for both models follows that of the isotropic matrix (the purple dashed circle) for angles  in the range of approximately  < ∕8, as well as  > 3∕8.At these angles, the resolved shear stress on the films is  too small for them to be activated significantly.However, as the angle between the film normal and the applied load reaches approximately  = ∕8 (or  = 3∕8), a transition occurs to a regime in which the softer films govern the onset of plasticity and the yield surface hence follows the dashed green curve, which represents yielding of the films.
As  → ∕4 the applied load becomes more and more effective in activating the film sliding mechanism, resulting in a significant drop of the effective yield stress.The lowest level is observed for  = ∕4, where τ0 ≈ 2 f 0 ∕(  regime, the two numerical models agree perfectly, while the analytical estimate underestimates the yield stress marginally due to the fact that it neglects the matrix.
To assess the dependence of the yielding anisotropy on the mechanical properties of the two plastic mechanisms, the strength ratio,  m 0 ∕ f 0 , is varied.This is done by increasing the initial yield strength of the matrix to  m 0 = 600 MPa and 800 MPa, while keeping the initial yield strength of the films constant.Fig. 5 compares the effective yield surfaces predicted by the effective laminate model for these different levels of contrast.In shear dominated load cases, i.e. around  = ∕4, the initial yield of the model is independent of the strength of the matrix, whereas in the region where the matrix is dominant the yield strength of the laminate scales with the strength of the matrix.The contrast ratio  m 0 ∕ f 0 hence directly controls the degree of anisotropy of the yield surface.

Influence of the film spacing and hardening
Next, the influence of the film spacing, , on the response of the laminate is investigated.A simple shear deformation parallel to the films is applied, similar to Section 3.2, and the stress-strain response is computed for two different values of the film spacing , while keeping the properties of the films constant.We characterize the film spacing by the ratio  = ∕, where  is the film thickness in the reference model and  denotes the volume fraction of the films.The values considered are  = 0.01 and 0.05, corresponding to  = 100  and 20 , respectively.For each case, two sets of simulations have been done: (i) without hardening, based on the Table 1, and (ii) with hardening incorporated in the response of the matrix films.In the latter case, the hardening modulus of the films,  0 , and the matrix, ℎ m 0 , were set to  0 = 2 f 0 ∕ and ℎ m 0 = 2 m 0 , respectively.For both phases the saturation stress was set to  ∞ = 3 0 and  = 1.5 was used.The response of the effective laminate model and the reference model practically coincides for each of the cases considered.(ii) The initial yield point of the laminate is insensitive to the volume fraction and to the degree of hardening.This is because for this particular loading, i.e. shear applied parallel to the films, the planar film deformation mode is the softer mechanism in the system, and it starts to yield at τeq ≈  f 0 , irrespective of the volume fraction or subsequent hardening.(iii) A minor increase in the initial yield stress is noticeable when the film volume fraction is reduced from  = 0.05 to 0.01.It is caused by the fact that the strain rate experienced by the films is amplified by a factor of 5 due to the reduction of volume fraction.(iv) When hardening is introduced, a secondary yield point occurs when the films have hardened sufficiently for the matrix yield strength  m 0 to be reached.From this point onwards the matrix flows as well and the overall hardening hence is slower.(v) As the volume fraction of the films is reduced from  = 0.05 to 0.01, the overall hardening rate increases, and the secondary yield point is reached faster.This effect has the same origin as mentioned above: for a smaller volume fraction of the films, experience a higher rate of deformation for the same overall deformation rate applied.

Application to lath martensite
Lath martensite is one of the main phases in advanced high strength steels, and plays a key role in achieving their mechanical properties.It has a hierarchical microstructure that forms upon quenching of the austenitic phase in low-alloy steels (Morito et al., 2006(Morito et al., , 2003)).Through a diffusionless phase transformation, FCC austenite transforms to  ′ martensite with a BCC or body-centered tetragonal (BCT) crystal structure (Bhadeshia and Honeycombe, 2017).First, at the scale of a single crystal, elongated martensite laths form.Each group of martensite laths with a particular crystallographic orientation is called a variant or sub-block.Laths of martensite with a very low misorientation constitute blocks, and several blocks sharing the same habit plane form a packet.Up to four packets of martensite can be created depending on the size of the prior austenite grain.Due to the nature of the martensitic transformation, a distinct crystallographic orientation relationship exists between the parent austenite and the transformed martensite.To describe this orientation relationship, the following pairs of relationships between planes and directions in the austenite and martensite phases have been proposed, respectively: Kurdjumov-Sachs (KS) {111}  ∥ {011}  ′ and ⟨101⟩  ∥ ⟨111⟩  ′ (Kurdjumow and Sachs, 1930), Nishiyama-Wassermann (NS) {111}  ∥ {011}  ′ and Fig. 6.Macroscopic stress-strain response of the effective laminate model and the reference model under simple shear applied parallel to the films.The solid curves show the response of the laminate model, whereas solid black lines with square markers show the behavior of the reference model.Two different volume fractions are considered, and for each case, a simulation with and without hardening was carried out.
There have been experimental reports on the existence of thin austenite films retained between the laths of martensite, due to an incomplete martensitic transformation, in martensitic steels (Morito et al., 2011;Araujo et al., 2021) and dual-phase (DP) steels (Liao et al., 2010;Yoshida et al., 2015).These retained austenite films are expected to be softer than the martensite laths.It has been argued that they may facilitate an inter-lath boundary sliding mechanism akin to the film sliding in the effective laminate model presented here (Maresca et al., 2014a(Maresca et al., , 2016b;;Liu et al., 2021).Accordingly, we apply the model in this section to lath martensite with retained austenite films.

Reference model
Our reference model in this section is a periodic two-phase laminate composed of alternating BCC single crystals (the laths) and FCC single crystals (the retained austenite).The same periodic cell shown in Fig. 2 is used to account for the interaction of the adjacent crystals.The ratio of the thickness of the austenite films to the thickness of the martensite laths is reported in the literature to be 4%-10% (Morito et al., 2011;Liao et al., 2010).Here, consistent with the previous section, we take the volume fraction of austenite to be  = 0.05.Unlike the previous section, the crystallography of the two phases is accounted for by modeling their mechanical response using single crystal plasticity.The KS orientation relationship for variant 1, {111}  ∥ {011}  ′ parallel plane, and ⟨ 101⟩  ∥ ⟨ 11 1⟩  ′ parallel direction is taken into account (Morito et al., 2003).This is done by aligning the direction of the variant parallel to the -direction and perpendicular to the -direction.
The plastic velocity gradient of the crystal is computed according to the standard crystal plasticity formulation (Peirce et al., 1983;Bronkhorst et al., 1992) where   0 =   0 ⊗   0 is the Schmid tensor related to the th slip system, obtained by the dyadic product of   0 , the slip direction, and   0 , the normal to the slip plane in the intermediate configuration.Similar to Eq. ( 9), the plastic slip rate on the active slip systems, ̇ , is given by a rate-dependent power-law relation, in which ̇ 0 denotes the reference slip rate and  is the rate sensitivity exponent.The plastic deformation is governed by the resolved shear stress,   = ( e ⋅  e ) ∶   0 .In order to incorporate hardening in the model, an evolution equation for the slip resistance,   , is formulated as together with an initial resistance   0 .The hardening modulus ℎ  evolves due to self hardening of the slip system  and latent hardening induced by other systems, , according to Here, ℎ  0 denotes the reference hardening modulus,   ∞ the saturation resistance,  the latent hardening ratio, and   the Kronecker delta.
For the martensite lath, two cases with either one {110}  ′ or two {110}  ′ and {112}  ′ slip families are simulated.The same parameter set is used for both BCC slip families.For the FCC slip, the usual {111} family is considered, and a (111) plane is oriented parallel to the interface (habit plane).The BCC slip systems are rotated to satisfy the KS orientation relationship with the FCC films (Liu et al., 2022).The material parameters of the BCC and FCC crystals are taken arbitrarily and given in Table 2.No hardening is considered in what follows and the values of   ∞ and  are hence irrelevant.

Effective laminate model -orientation factors
In comparing the predicted response of the effective laminate model with that of the crystal-plasticity based reference model, once should realize that the latter takes into account the full plastic anisotropy due to crystallographic orientation, whereas the effective laminate model assumes isotropy for the matrix and planar isotropy for the films.For a fair comparison, the constitutive parameters used for the effective laminate model should reflect the average response of the respective single crystals characterized by the slip system level parameters of Table 2 for all relevant orientations.
For the isotropic matrix, we can resort to the well-known Taylor orientation factor (Rosenberg and Piehler, 1971).It links the overall uniaxial yield stress of a polycrystal to the slip resistance of an individual slip system.For BCC crystals, this factor is reported to be 2.45, based on a CPFEM estimation with one slip family subjected to uniaxial tensile load (Zhang et al., 2019).However, the formulation of the effective laminate model is based on the shear-equivalent von Mises stress rather than that equivalent to uniaxial tension.Accordingly, we employ an orientation factor equal to  m = 2.45∕ √ 3 ≈ 1.41.We use this factor as is -no attempt has been made to optimize it for the cases considered below.The Taylor factor  m is incorporated in the isotropic part of the effective laminate model by defining the plastic properties of the matrix phase as In case hardening is taken into account, the same factor  m should be applied to  BCC ∞ and ℎ BCC 0 , analogous to (24)a.For similar reasons, an orientation factor must be taken into account for linking the FCC slip system properties for one slip system to the planar isotropic plasticity model of the films.Note, however, that the situation is slightly different here than for the matrix, because the orientation of the slip plane associated with the relevant slip systems is identical to the orientation of the films, and hence the effect of random slip plane orientation does not need to be accounted for in a Taylor-like factor.What remains to be accounted for is the random orientation of the slip directions in that plane.Without an orientation factor, the planar film model would assume that a slip system is always available in exactly the direction of maximum shear traction.In reality, as reflected by the full crystal plasticity model, the 'best available' slip direction may be at certain finite angle with respect to the shear traction, resulting in a slightly harder response.The orientation factor for the film,  f , accounts, on average, for this misalignment.It has been determined by computing the average response of an FCC slip plane with 3 slip systems oriented at 60 • with respect to each other under an applied simple shear load parallel to the plane.Consequently, the average orientation factor has been computed to be  f = 1.1.This factor enters the film plasticity properties of the effective laminate model via and conditions akin to (25)a for the hardening parameters.

Single packet response and comparison to the fully resolved reference model
The overall stress-strain response of the laminate model, plotted using solid red lines, is compared to the response of the reference model, plotted as solid black lines with square markers, in Fig. 8. Three different load cases are applied.First, simple shear along the habit plane in -direction, i.e.   , parallel to the direction of the first variant, i.e. ⟨ 101⟩  ∥ ⟨ 11 1⟩  ′ .Second, tension in -direction, perpendicular to the habit plane in the direction of {111}  ∥ {011}  ′ .Third, due to the anisotropy of the reference model when it is loaded in tension along and perpendicular to the {011}  ′ plane, its tensile response along the habit plane is plotted as well.The laminate model gives the same response for the two tension load cases, i.e. the single red curve in Fig. 8 for both.
Fig. 8 shows that under shear applied parallel to the austenite films (in the direction of the in-plane slip system in the austenite films) the responses of the two models coincide.However, under tension along and perpendicular to the films, the responses are inconsistent.The main origin of the deviations is the assumption of isotropy in the matrix, requiring the Taylor factor in the laminate model.This factor has been calculated based on the average response of a sufficiently large polycrystal, i.e. it accounts for the average of all directions.In the case of the reference model subjected to a tensile load along the film direction (and {011}  ′ of BCC), the FCC films are inactive, and the Schmid factor for each direction of applied load for the BCC crystal is different than the incorporated averaged value,  m = 2.45∕ √ 3, which hence leads to significant quantitative differences.For the value of  m adopted here, the response of the laminate model is close to (within 5% of) that of the crystal plasticity based reference model for out-of-plane tension, but it underpredicts the flow stress under in-plane tension by approximately 25%.Adapting the value of  m would allow one to compromise between these approximations, but the effective laminate model will never be able to capture the anisotropy exhibited by the reference model.Note, however, that it does capture the main source of anisotropy, i.e. that due to flow of the retained austenite films, as evidenced by the factor of approximately 3 in the flow stress under tension vs that under shear.
Fig. 9 compares the yield surface prediction of laminate model to that of the reference model.The BCC crystal in the reference model is simulated with one and two slip system families.It is observed that the effective laminate model captures the main anisotropy of the material, stemming from the activity vs inactivity of the FCC films.In shear dominant loads, there is a close agreement between the prediction of the laminate model and the reference model.However, in directions close to tension, consistent with our observations in Fig. 8, differences of up to one third occur.The accuracy of the laminate model in predicting tension-dominated loading increases with a larger number of slip systems, which is a logical outcome, given that an isotropic material may be viewed as a crystalline material possessing an infinite number of slip systems, in all directions.

Mesoscale simulation of lath martensite in DP steel using the effective laminate model
We now apply the laminate model at the scale it is intended for: the polycrystalline mesoscale.At this scale, many individual grains and packets can be distinguished, yet fully resolving their substructure is computationally prohibitive.To carry out the mesoscale RVE simulations, the laminate model has been implemented in the DAMASK software (Roters et al., 2019).DAMASK is a unified multi-physics simulation package that incorporates a variety of constitutive models and homogenization approaches.Our simulations employ DAMASK's spectral solver.This implies that, instead of by the more common finite element shape functions, the problem domain is discretized by Fourier basis functions (i.e.sine and cosine functions), and the Fast Fourier Transform and its inverse are employed to render the equilibrium iterations highly efficient.The reader is referred to Roters et al. (2019) for more details on the solution methodology.
A synthetic two-dimensional DP steel microstructure with randomly distributed ferrite and martensite grains, as shown in Fig. 10, is considered.The volume fraction of the martensite grains is 40%, and each martensite grain is considered to be a single crystallographic packet, i.e., having a unique habit plane orientation.The isotropy assumed for the (martensite lath) matrix and the planar isotropy of the (retained austenite) films in the effective laminate model imply that the model cannot distinguish individual variants within a martensite packet.The orientation relationship between neighboring packets stemming from the same prior austenite grain is therefore not taken into account.Three simulations are carried out wherein the martensite grains are modeled and compared using three different constitutive models: (i) isotropic visco-plasticity as defined in Eqs. ( 1)-( 6), without the second term in Eq. ( 3) and a Taylor factor of  m = 2.45∕ √ 3, (ii) crystal plasticity (1 slip family) as defined in Eqs. ( 20)-( 23), and (iii) the effective laminate model.The crystal plasticity case is considered to investigate the effect of martensite texture on the anisotropy of the microstructure, and compare the predicted BCC-only response to that of the effective laminate model.The ferrite grains are modeled using the isotropic visco-plasticity model, as defined in Eqs. ( 1)-( 6), without the second term in Eq. ( 3), in all of the simulations.The material parameters for ferrite, retained austenite and lath martensite are adopted from Maresca et al. (2014bMaresca et al. ( , 2016a)).The volume fraction of retained austenite within the martensite grains is taken to be  = 0.05, which is in agreement with reports on lath martensite in DP steels (Yoshida et al., 2015;Liao et al., 2010).A uniaxial tensile load in the -direction, similar to Eq. ( 19), is applied to the periodic RVE, in which the stretch factor  is incremented from zero to  = 1.1 at a rate of λ = 0.01 s −1 .
The diagram on the right in Fig. 10 shows the macroscopic stressstrain response computed in the three simulations done on the DP steel microstructure shown on the left.The responses of the cases in which martensite is modeled via CP and isotropic plasticity are almost identical.However, an approximately 10% softer response is observed for the simulation with the effective laminate model.The thin soft films parallel to the habit plane of lath martensite entail softer martensite grains and hence a softer material response.Note that this mechanism is absent in both the isotropic model and the CP model and we hence consider the laminate model to be a better model in this sense.The negligible difference between the isotropic and CP models gives us reason to believe that the isotropy assumptions made in the effective laminate model are inconsequential at the present scale.
While the effect of the soft austenite films on the global response is already significant, an even more profound influence is observed on the local deformation in the microstructure.The equivalent strain maps of the ferrite and martensite grains at a global stretch factor of  = 1.1 are compared for the three simulations in the top and bottom row of Fig. 11, respectively.It is observed in the top row that the strain distributions in the ferrite are qualitatively similar for all the cases, i.e. the same localized areas are observed.This pattern also matches in a quantitative sense for the simulations in which the isotropic and CP models (both without retained austenite films) are used in the martensite.However, it is shown that by incorporating the effective laminate model in the martensite grains, the strain localization in the ferrite grains is reduced.An example of this behavior is highlighted by an ellipse around a ferrite channel in the top row of Fig. 11.As a result of activating the sliding mode in lath martensite, this phase can accommodate more plasticity in the direction of the films, leading to a decrease in the amount of localization in the ferrite grains.This observation is supported by the bottom row of diagrams in Fig. 11, where the spatial distribution of the strain in the martensite is compared for the three cases.The highlighted region (dashed circle in the bottom row) is one of the areas exhibiting a higher contribution of the martensite grains to plasticity when the effective laminate model is used.
The above observations are confirmed in Fig. 12, where the distributions of stress and strain levels observed in the three simulations are compared.It is seen that when the effective laminate model is used for the martensite grains, shown by solid lines, the stress and strain partitioning is reduced, i.e. closer mean values are observed for the ferrite and martensite grains in both diagrams.The highest scatter of the strain distribution in the martensite grains is observed when the effective laminate is used for the martensite.This agrees with the observation made by Morsdorf et al. (2016), based on their Fig. 3, that the heterogeneity of the strain distribution in the martensite grains is not correlated with the heterogeneity of the corresponding crystallographic orientation distribution, i.e. the Taylor factor map.Moreover, Fig. 12 shows that the average local stress is the lowest for martensite described by the effective laminate model.This results from the embedded soft films, even though the heterogeneity of the stress distribution is higher for martensite grains modeled with CP.The significant scatter of the stress for the CP case is due to variations in the Schmid factor, which ranges from 0.22 to 0.5 for a BCC crystal with one slip family.

Computational efficiency
Finally, the computational efficiency of the effective laminate model is compared to the other two cases studied above.The model proposed by Maresca et al. (2016b) is considered as well for the comparison of the CPU time, since, to our best knowledge, it is the only existing model in the literature accounting for the boundary sliding in the lath martensite which is computationally tractable in mesoscale simulations.In that model, the lath martensite microstructure is modeled as a periodic stack of martensite laths and retained austenite films.The model incorporates isotropy for the martensite laths and the out-ofhabit-plane systems of the retained austenite films, whereas the three in-habit-plane slip systems of the retained austenite are modeled via a standard crystal plasticity approach.
Fig. 13 shows the CPU time for all of the simulations normalized by that of the isotropic plasticity simulation.The data for the Maresca et al. (2016) model has been estimated based on simulations reported in Maresca et al. (2016a) of a DP microstructure with 33% of martensite.However, since the models are implemented in different solvers (Maresca et al. (2016b) uses a standard finite element solver, while this work uses a spectral solver), and the efficiency of the codes can be different, the CPU time comparison is only indicative, and no sharp quantitative conclusions can be made.
The effective laminate model proposed in this paper has a computation time comparable to that of isotropic plasticity while it preserves the physics of the boundary sliding and anisotropic plasticity observed in the martensite.The model proposed by Maresca et al. (2016b) has a higher computational time due to the use of a lamella rule of mixtures, which requires satisfaction of equilibrium and kinematic compatibility at each time increment at the interface between the austenite films and the martensite laths.The effective laminate model adopts an iso-stress assumption for the embedded films and hence equilibrium and compatibility are satisfied trivially.In essence, since the films can only experience in-plane shear, other out-of-plane components of deformation, such normal components, will not be relevant for the films.Moreover, the anisotropy in the habit plane which is considered in the model of Maresca et al. (2016b), is disregarded in the effective laminate model, as it has minor influence on the strains induced by the sliding mode.Compared with the conventional CP model, the effective laminate model is significantly faster, even though it incorporates the additional mechanism of inter-lath sliding.

Conclusions
A computationally efficient micromechanics based model is proposed to capture the anisotropic behavior of a particular class of lamellar microstructures in which thin soft lamellae can be condensed into discrete slip planes.These planes are embedded in a matrix which represents the harder lamellae.Accordingly, the model is summarized as an isotropic visco-plastic model which is enriched by an additional  planar plastic mode parallel to the planes.A comparison with a reference model, which is a two-phase laminate sharing the same geometric and material properties, shows that the model recovers the same response under shear applied parallel and tension applied perpendicular to the films.The yield surface of both models shows the highest yield limit for tension applied parallel or perpendicular to the soft films.The weakest response is observed when the applied load entails a maximum shear stress which is aligned with the soft films.
The formulation is applied to model the substructure of lath martensite with thin inter-lath retained austenite films.It is shown that in shear-dominated cases the yield strength of the effective laminate model coincides with the fully resolved crystal plasticity reference model.The differences observed for tension-dominated loads are due to the isotropic approximation made for the anisotropic martensite laths.The effective laminate model is used in mesoscale simulations of a DP microstructure, and the results are compared to the cases in which the martensite grains are modeled via isotropic plasticity and crystal plasticity, without the effect of the softer plastic mechanism aligned with the retained austenite films.In the case of martensite grains modeled with the effective laminate model, as a result of the embedded soft films, the contribution of the martensite grains to the plastic deformation is increased, and hence a relatively homogeneous stress and strain distribution is observed in the microstructure.For the same reason, a lower macroscopic yield stress of the DP microstructure is observed.The computational gain of the model is investigated, whereby it is demonstrated that the computationally efficiency of the model is almost at the same level of the isotropic visco-plastic model, while it still captures the physics of inter-lath sliding and extreme plastic anisotropy in the lath martensite.
The main limitation of the model in its current form is the assumption of isotropy of, in particular, the (martensite lath) matrix.We are currently developing a version of the model in which both the matrix and the film incorporate the discrete slip systems of the respective crystals.It is anticipated, however, that the computational cost of this model will be significantly higher -at least comparable to that of conventional crystal plasticity.Furthermore, the model will require the characterization (for real microstructures) or generation (for synthetic ones) of the full crystallographic orientation of all packets or subblocks.Whereas this investment may pay off for modeling and studying microscale experimental observations, it remains to be seen whether it has added value for the type of mesoscale simulations discussed in Section 5.4.Based on the results shown in this contribution, we believe that the present model may well be the more appropriate platform for developing a comprehensive understanding of the effect of mesostructure on the plasticity and, in the future, damage in multiphase Advanced High-Strength Steel grades.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: V. Rezazadeh, J.P.M. Hoefnagels, M.G.D. Geers, R.H.J. Peerlings report financial support was provided by M2i Materials Innovation Institute.V. Rezazadeh reports administrative support and equipment, drugs, or supplies were provided by Tata Steel.

Fig. 1 .
Fig. 1.(a) Sketch of a fully resolved laminate consisting of a matrix phase with embedded soft thin films.(b) Simple shear applied parallel to the films resulting in the accommodation of all slip by the films only.(c) Tension applied perpendicular to the films deforming only the matrix phase.The proposed model is able to discern the distinctive behaviors exhibited by the two-phase laminate under consideration.

Fig. 2 .
Fig. 2. The periodic cell representing the reference model -an infinite laminate with alternating layers representing the matrix and embedded thin films.The dashed lines indicate the periodicity of the geometry.The volume fraction of the film is  = ∕ = 0.05.The reference model is used to assess the response of the proposed effective laminate model.

Fig. 3 .
Fig. 3. (a)The macroscopic stress-strain response of the laminate model, plotted as solid lines, compared to that of the reference model, plotted as solid lines with square markers.(b) Activity of the two plastic modes, i.e. plastic deformation of the matrix and the film, under simple shear parallel to the films, and (c) under uniaxial tension applied perpendicular to the films.The volume fraction of the films is  = 0.05 and the applied deformation rate 0.01 s −1 .

Fig. 4 .
Fig. 4. Polar plot comparing the yield limit of the effective laminate model, shown by the red solid line, and the reference model, shown by black solid line overlaid with squares.The purple dashed circle indicates the stress at which the matrix phase starts to yield, i.e.  eq =  m 0 , whereas the green dashed curves indicate when the shear stress acting on the films equals the initial yield stress of the films,  =  f 0 .The radial axis is the yield stress in MPa, and the angular axis indicates the angle between applied tension and the normal to the films.The volume fraction of the films in the reference model is  = 0.05.
230 MPa -almost a factor of two lower than for  = 0 or  = ∕2.Throughout this film shear dominated V.Rezazadeh et al.

Fig. 5 .
Fig. 5. Polar plot comparing the yield stress of the effective laminate model for three values of the strength ratio,  m 0 ∕ f 0 .The volume fraction of the films is  = ∕ = 0.05.

Fig. 6
compares the macroscopic stress-strain responses of the laminate model, all plotted as solid red lines, with the responses of the reference model, plotted as solid black lines overlaid with square markers.Four different responses are shown for each model, as we consider two different film spacings (volume fraction in the reference model) and for each case a simulation with and without hardening is carried out.The observations made from Fig. 6 are as follows.(i)

Fig. 7 .
Fig. 7. Schematic of an individual prior austenite grain and the lath martensite substructure formed due to the martensitic phase transformation.An illustration of six KS orientation variants is given for one packet.

Fig. 8 .
Fig.8.The stress-strain response of the laminate model, plotted using solid red lines, and of the reference model, plotted as solid black lines with square markers.Three load cases are shown: simple shear along the habit plane and tension along and perpendicular to the habit plane.The volume fraction of the films is  = 0.05 for all cases.

Fig. 9 .
Fig. 9. Comparison of the yield surface of the effective laminate model with that of the crystal plasticity based reference model of lath martensite.The BCC crystal is simulated with only the {110} slip system family, marked by black square markers, and with two slip families ({110} and {112}), marked by yellow diamonds.The volume fraction of the films is  = 0.05 for all cases.

Fig. 10 .
Fig. 10.(left) Synthetic two-dimensional DP steel microstructure with martensite grains shown in red and ferrite in green.The grain boundaries are shown in black.The martensite volume fraction is 40%.A random crystallographic orientation is assigned to the martensite grains, and each individual grain is considered to be a single packet, i.e. with one particular orientation of the film as used in the effective laminate model.(right) The macroscopic stress-strain response of the microstructure computed with three different constitutive models for the martensite grains.

Fig. 11 .
Fig. 11.Computed equivalent plastic strain in the ferrite and martensite grains at an applied global deformation of  = 0.1.The three columns show the different constitutive models used in the martensite grains.The top row shows the ferrite grains, the bottom row the martensite.In all of the simulations ferrite is modeled by isotropic visco-plasticity.

Fig. 12 .
Fig.12.Probability density  of (a) the equivalent true strain  eq , and (b) the equivalent true stress  eq in the ferrite and martensite grains plotted for each of the three models used for lath martensite.The distributions were obtained for the final deformation, at  = 1.1.

Fig. 13 .
Fig. 13.Normalized CPU time for the DP microstructure shown in Fig. 10 using different constitutive models implemented for the martensite grains.The data for the Maresca et al. (2016) model is an estimate based on work reported inMaresca et al. (2016a).The CPU time reduction due to the use of a different solver, spectral versus finite elements, is not accounted for in this comparison.

Table 1
Model parameters used for simulating the response of the reference model and the proposed effective laminate model.The elastic constants of the film are relevant only for the reference model.

Table 2
Crystal plasticity parameters used in the lath martensite simulations for the retained austenite films (FCC) and the martensite laths (BCC).