A Multiscale Method for Damage Analysis of Quasi-Brittle Heterogeneous Materials

A novel multiscale algorithm based on the higher-order continuum at both microand macrostructural level is proposed for the consideration of the quasi-brittle damage response of heterogeneous materials. Herein, the microlevel damage is modelled by the degradation of the homogenized stress and tangent stiffness tensors, which are then upscaled to govern the localization at the macrolevel. The C1 continuity finite element employing a modified case of Mindlin’s form II strain energy density is derived for the softening analysis. To the authors’ knowledge, the finite element discretization based on the strain gradient theory is applied for the modeling of damage evolution at the microstructural level for heterogeneous materials for the first time. The advantage of the novel C1 finite element formulation in comparison with the standard finite element discretization in terms of the regularization efficiency as well as the objectivity has been shown. An isotropic damage law is used for the reduction of the constitutive and nonlocal material behaviour, which is necessary for the physically correct description of the localization formation in quasi-brittle materials. The capabilities of the derived finite element to capture the fully developed localization zones are tested on a random representative volume element (RVE) for several different loading cases. By employing the conventional second-order computational homogenization, the microstructural material constitutive response is averaged over the whole RVE area. In order to model the loss of structural integrity when sharp localization is formed across RVE, the specific conditions which detect a completely formed localization zone are developed. A new failure criterion at the microstructural level has been proposed. The derived finite element formulation, as well as the multiscale damage algorithm, are implemented into the finite element program ABAQUS. The capabilities of the presented multiscale scheme to capture the effects of the deformation localization are demonstrated by few benchmark numerical examples.


Introduction
Physically, formation of damage and macroscopic cracks is a direct result of the cascade of events happening at the microstructural level [Lemaitre (1992); Murakami (2012) ;Liu, (2008) ;Sofi, Bishay and Atluri (2018); Sane, Padole and Uddanwadiker (2018)]. Given local mesh refinement techniques are mainly efficient if only a mild and small localization zone is expected, position of which is known a priori. However, the most progress in the multiscale modelling of the microscopic localization is achieved by utilization of CH methods. Through the homogenization process, response of the heterogeneous microstructure is averaged over an RVE, whereby a new, effective homogeneous material is formed. Classical homogenization techniques are built upon the principle of separation of scales, which states that the RVE size should be much smaller, than a characteristic length over which the macroscopic loading varies in space. In other words, the uniform distribution of the macro-strain over the entire RVE domain is assumed. This, however, is violated when the first-order CH schemes, which are based on classical continuum formulation at both scales, as described in Kouznetsova et al. [Kouznetsova, Brekelmans and Baaijens (2001); Miehe and Koch (2002); Temizer and Zohdi (2007)], are used with the problems where strain softening occurs at the microlevel. When there is no clear separation of scales, capturing of the propagation of the underlying rapidly fluctuating responses can be remedied to some extent by higher-order enrichment of the macroscopic continuum. Besides, standard continuum formulation at the macroscale cannot regularize the formation of the strain localization, which in addition leads to the ill-posedness of the macrostructural BVP. As an improvement to the first-order CH, second-order CH is proposed [Kouznetsova, Geers and Brekelmans (2004b)], which is shown to be successful in treating only the mildly softening materials, specifically the materials not exhibiting the deformation beyond a quadratic nature in the displacements, as stated in Geers et al. [Geers, Kouznetsova and Brekelmans (2010)]. Classical homogenization in its essence implies the averaging of some physical phenomenon, and it is believed that, in the case of the sharp localization which is characteristic to the certain RVE, it should not be performed. Additionally, with the occurrence of the sharp strain localization, homogenized response stops being objective with the respect to the size of the RVE-by increasing the size of the micro-sample, the macroscopic structural response becomes more brittle [Gitman, Askes and Sluys (2007)]. In that case RVE stops being statistically representative for the macroscopic material point and should be called a microstructural volume element (MVE) instead, as stated in Coenen et al. [Coenen, Kouznetsova and Geers (2011)]. Another class of homogenization methods which deal with the strain softening problems is based upon the enrichment of the macroscale continuum with a discontinuity, where the microscale strain localization band is lumped into a macroscale cohesive crack. Such way of thinking emanates from the standpoint that the failure of the MVE cannot be extrapolated to the neighbouring material, resulting in the introduction of the equivalent discontinuity at the macroscale integration point, rather than continuous display of localized deformation. Taking into account the techniques used for the extraction of the equivalent discontinuity from the localized MVE and formation of the corresponding macrostructural effective constitutive relation, several different procedures can be found in the literature. For example, this can be done by assuming that continuum which is not comprehended by the localization behaves as a liner-elastic material [Massart, Peerlings and Geers (2006) ;Verhoosel, Remmers and Gutiérrez (2009);Nquyen, Lloberas-Valls, Stroeven et al. (2011)], by employing the additional RVE for the homogenization of the bulk material [Belytschko, Loehnert and Song (2007); Nquyen, Stroeven and Sluys (2012)], or by employing the same MVE for the extraction of both the discontinuity and the stress-strain response of the neighbouring material [Coenen, Kouznetsova and Geers (2012)]. The existence of an RVE for softening quasi-brittle materials undergoing localized damage has been confirmed in Nquyen et al. [Nquyen, Lloberas-Valls, Stroeven et al. (2010)], where a new averaging technique based on extraction of the deformation of just a localization band is proposed. By using this technique, a CH scheme for discrete macroscopic crack modelling that is objective with respect to the size of the RVE is presented in Nquyen et al. [Nquyen, Lloberas-Valls, Stroeven et al. (2011);Nquyen, Stroeven and Sluys (2012)]. It is clear from the given overview that development of the multiscale models which deal with the damage localization problems still needs to be improved in order to enable their employment in practical, engineering considerations. When it comes to advanced multiscale schemes, a new second-order CH scheme is derived recently in Lesičar et al. [Lesičar, Tonković and Sorić (2017)], where the C 1 continuous finite elements are employed at both macro-and microlevel. Employment of the nonlocal theory at the microscale has shown better efficiency compared to available homogenization schemes, additionally offering an advanced frame for damage modelling. A new damage model employing the strain gradient theory embedded into C 1 continuous finite element is recently presented in Putar et al. [Putar, Sorić, Lesičar et al. (2017)], where the exceptional regularization capabilities of such model are demonstrated. Not only the numerical results obtained are independent on the discretization, but also the spurious damage growth reported in Geers et al. [Geers, de Borst, Brekelmans et al. (1998);Simone, Askes, and Sluys (2004); Poh and Sun (2017)] is completely eliminated. As shown in Poh et al. [Poh and Sun (2017) ;Putar, Sorić, Lesičar et al. (2017)], due to decrease in the size of the microstructural interaction domain with increasing loading, a physically acceptable formation of the localization of the deformation can be described, where a macrocrack comes into existence from the initially scattered network of microcracks. The topic of this paper is the consideration of utilization of the damage model presented in Putar et al. [Putar, Sorić, Lesičar et al. (2017)] at the microstructural level of the twoscale scheme presented in Lesičar et al. [Lesičar, Tonković and Sorić (2017)], and the main objective is to demonstrate the algorithm's ability to describe the formation of the localization zone at the macrostructural level. Both scales are described as a higher-order continuum, which necessitates the employment of the C 1 continuity. Since the nonlocality is intrinsically included in the theory, neither objectivity issues related to discretization nor spurious damage growth are expected, i.e., a full regularization of the localization problems at both scales should be achieved. A following approach is employed for the consideration of the microstructural damage. Due to nature of the macrocrack formation in the quasi-brittle materials which is preceded by a diffuse network of microcracks, an assumption can be made that the material locally has a similar structural response. In that sense, the formation of the sharp localization at microlevel does not necessarily mean the initiation of the macrocrack, but a situation where the macrostructural material point loses its stiffness while still remaining a part of the continuum. By doing so, the conventional homogenization can be applied to obtain the averaged stiffness behaviour and stresses from the observed damaged RVE until the full formation of the localization zone, when the material conditions that account for the loss of material integrity have to be applied at the appropriate macrolevel integration point. During the evolution stage of the microstructural localization, the homogenization process can be interpreted as the means for transformation of the localized deformation band and the rest of elastically unloading area of an RVE, into the equivalent RVE where damage is homogeneously dispersed, as represented graphically in Fig. 1. The topic regarding the preservation of the objectivity with respect to the size of the RVE is not included in this consideration.

Figure 1: Homogenization of the localized deformation
The paper is organized as follows. In Section 2 the finite element developed in Putar et al. [Putar, Sorić, Lesičar et al. (2017)] is modified in order to coincide with the underlying continuum theory as presented in Lesičar et al. ]. Its capabilities in describing the localization are tested on an RVE example for three different loading cases. In Section 3 the basic relations of the existing two-scale framework that employs the gradient continuum at both structural levels are given. In addition, RVE cracking conditions, i.e., conditions that indicate the formation of the localization band at the RVE level and consequential change of the macrostructural stiffness are derived. The algorithm's capability to provide mesh objective results are given in Section 4 by employing the homogenous microstructure. Afterwards, a set of benchmark examples is presented, where the macrostructural localization of the deformation is achieved by utilizing the heterogeneous microstructure. Section 5 is dedicated for some concluding remarks.

C 1 continuity finite element for softening analysis of microstructure 2.1 Weak formulation
In the work presented in Lesičar et al. ] both macroand microlevel are discretized by C 1 continuity finite elements based on the modified Mindlin's form II strain energy density as described in Altan et al. [Altan and Aifantis (1992); Ru and Aifantis (1993)]. In order to make use of this setting, a modification of the finite element based on strain gradient theory derived in Putar et al. [Putar, Sorić, Lesičar et al. (2017)] and presented in Fig. 2 has to be made. The element consists of three nodes and 36 degrees of freedom with the displacement field approximated by the condensed fifth order polynomial. The nodal degrees of freedom are the two displacements and their first-and second-order derivatives with respect to the Cartesian coordinates. The physical interpretation of the mentioned nodal degrees of freedom is comprehensively described in Lesičar et al. ]. In comparison to the constitutive model described in Lesičar et al. [Lesičar, Tonković and Sorić (2014)], in this manuscript only the classical constitutive tensor σε C is employed, while the higher-order variables are expressed as the gradients of the strain field in form of   x , while v is the vector of the nodal degrees of freedom. Interpolation matrix ε B contains adequate first derivatives of the element shape functions N , needed for the description of the element displacement field u in form The strain tensor ε , defined as is coupled with its work conjugate Cauchy stress tensor, given as The work conjugates of strain gradients The principle of virtual work is written in terms of strain gradient tensors as with s and A representing the perimeter and the surface of the element, while t and T stand for the traction and the double traction tensor, respectively. In addition to Eq. (8), the boundary conditions expressed by the displacement and the normal derivative of displacement ( ) ∇ ⊗ ⋅ u n should be prescribed to solve the boundary value problem. After the application of discretization given by Eqs. (1)-(4), Eq. (8) adopts the following form Considering a nonlinear problem described by (9), the displacement vector u, the stress tensor σ and the double stress tensors 1 x μ and 2 x μ are updated according to where the exponent ( ) 1 i − refers to the last converged equilibrium state, and the symbol ∆ indicates an incremental change and mathematically acts as a differential operator. Considering Eqs. (10)-(13), the incremental form of the principle of virtual work given in matrix notation reads as Linearized incremental constitutive relations are expressed by where l represents an internal length scale parameter which introduces the nonlocal material behaviour in the model. In the same manner as shown in Putar et al. [Putar, Sorić, Lesičar et al. (2017)], softening behaviour is implemented by employing the isotropic damage law given by where D is the scalar damage variable [Lemaitre and Chaboche (1990)], and eff C and C denote the effective and the elastic stiffness tensors, respectively. By inserting Eq. (18) into the non-linearized form of relations represented by Eqs. (15)-(17), damage enhanced constitutive model can be written in the following way: Next, by employing the following updates Implementation of the constitutive relations given by Eqs.
with the particular element stiffness matrices defined as and the external and internal forces e F and i F described by the first and second term on the right-hand side of Eq. (14), respectively. The algorithm derived above is implemented into the commercial finite element software ABAQUS/Standard [Abaqus (2014)] via user-defined subroutine UEL.

Numerical test
Applicability of the finite element for softening analysis derived in previous section is tested on an RVE of the side length of 0.5 mm representing a high-strength porous steel, depicted in Fig. 3. In order to examine the damage response, three elementary loading cases are applied on the RVE, namely, tensile, compressive and shear. The straight edges at the left and right RVE sides are enforced by suppressing the following degrees of freedom: 1,11 u , 1,22 u , 1,12 u , 2,12 u , 1,2 u and 2,1 u . The matrix material is considered homogeneous with the elastic behaviour described by the Young's modulus 5 2 2.1 10 N/mm E = × and the Poisson ratio 0.3 ν = , while the microstructural nonlocal behaviour is described by the internal length scale parameter 0.01 l = mm. For the description of the damage evolution, the exponential softening law, defined in Peerlings [Peerlings (1999)], is employed as where κ is the monotonically increasing scalar history parameter which governs the damage state and can be determined as the largest value of the equivalent strain eq ε reached in a material point during the loading history. 0 κ denotes the threshold equivalent strain needed the start of the damaging process, while the parameters α and β determine the stress decrease and the damage rate, respectively. For the description of the equivalent strain eq ε , modified von Mises' definition is used, given as where the parameter k represents the ratio between uniaxial compressive and tensile strength of the material, and 1 I and 2 J are the first invariant of the strain tensor and the second invariant of the deviatoric strain tensor, respectively, defined in de Vree et al. [de Vree, Brekelmans and van Gils (1995)]. Following parameters are employed for the description of material damage behaviour: 0 0.0002 mm is applied on the appropriate loading edges. The computational models and deformed shapes of the RVEs with damage contour plots observed in the failure stages for tensile, compressive and shear loading case are given in Figs. 4-6, respectively.   The load-displacement diagrams for the all three loading cases are depicted in Fig. 7. As evident, the compressive structural response is associated with much higher values of the reaction forces reached due to postponed initiation of damage. This is a consequence of the utilization of parameter 10 k = in the modified von Mises' equivalent strain definition given by Eq. (34), causing the elastic equivalent strain in compressive loading case to increase with the ten times slower rate compared to tensile loading case. Furthermore, as obvious from the presented figures, the derived finite element is capable of capturing the localization of the deformation over the RVE for different loading situations. In order to assess the influence of different internal length scales on the deformation responses, an additional analysis is made with a tensile loading case on the RVE employing two length scales, 0.01 l = mm and 0.02 l = mm. The structural responses for both internal length scales are shown in Fig. 8. It can be concluded that smaller nonlocal parameter results in more brittle softening, which is in accordance with the nonlocal theory and numerical results presented in Putar et al. [Putar, Sorić, Lesičar et al. (2017)]. In Figs. 9 and 10, evolutions of the damage localization zones are given for several different loading steps for internal length scales 0.01 l = mm and 0.02 l = mm, respectively.  A different formation of the localization zones can be noted in above figures, which is obviously a consequence of utilization of different internal length scales that define different nonlocal behaviours. When the nonlocal effect is wider, i.e., when a bigger area of neighbouring material is in the interaction, one predominant and wider localization zone is formed, as evident from Fig. 10. In contrast, for the weaker nonlocal effect, two narrower localization bands are formed, as shown in Fig. 9. This phenomenon would probably not exist in reality, since a crack would be initiated in a dominant localization zone, across which the fracture would eventually happen. The constitutive model presented in Eqs. (19)-(21) does not consider the crack formation, allowing in this way for the loading to be carried over the localization zone. As discussed in Putar [Putar (2018)], the definition of the damage variable should probably be written in terms of both strain and strain-gradient tensors in order to reduce the internal forces in the localization zone and overall reaction force response. However, in both examples the localization zones are fully formed, where no additional damage growth outside of it can be noticed. In other words, RVEs still have the load-carrying capacity, with the change of strain level allowed only in the middle of the localization band. This trait is important because it allows the simple assessment of the RVE failure, which will be incorporated into the multiscale scheme in the next section. From the presented numerical tests it can be concluded that the proposed C 1 continuity finite element formulation can successfully predict the damage evolution and it is suitable for the implementation in the multiscale algorithm.

Nonlocal multiscale scheme for damage analysis
Employment of the nonlocal continuum at the macrolevel instead of the classical one is beneficial for two main reasons. First, since both strain and strain gradients are used in formulation of the RVE boundary conditions, the more complex and realistic RVE deformation modes can be obtained, and consequently a more realistic prediction of the damage growth can be made. The second reason is concerned with the regularization capabilities of the strain gradient theory, which enables the transfer of the nonlocal parameters from the microstructure by homogenization of additional higher-order constitutive tensors. The C 1 -C 1 multiscale scheme is given in Fig. 11 and can be briefly described in the following way. Starting from the converged global macrolevel nonlinear boundary value problem (BVP), described by ∆u can be then formulated by the appropriate micro-to-macro scale transition relations, which will be discussed later. The microlevel BVP is formed and solved afterwards, followed by the homogenization of the resulting microlevel variables needed for the formation of the macrolevel constitutive behaviour. Here, the Cauchy stress tensor M σ , double stress tensor M μ and nine constitutive tensors M C have to be homogenized in the most general case. Once computed, the homogenized stresses and stiffness are transferred back to the macrolevel finite element integration point, where they contribute in the calculation of the finite element stiffness matrix k and internal forces vector i f . When the homogenized response is obtained in all integration points for all macrolevel finite elements, the global stiffness matrix K and internal force vector i F are formulated, and the updated macroscale BVP can then be established. From this point the process repeats in the loop until the convergence for the given loading conditions at the macrolevel is reached. The presented routine has already been described in more detail in Lesičar et al. ] for the consideration of linear-elastic material behaviour, where the nonlinear phenomena have not been considered. The contribution of this work is the implementation of the damage algorithm presented in Section 2 at the microstructural level of the multiscale scheme, while all other procedures remain the same. Figure 11: Scheme of C 1 -C 1 multiscale algorithm

Basic macro-micro scale transition relations
All relevant relations regarding the transition of variables from macro-to microscale are presented in Tab. 1 and will be briefly discussed here for the clarity reasons only. Again, a detailed derivation of presented relations can be found in Lesičar et al. [Lesičar, Tonković and Sorić (2017)]. Formation of the RVE displacement field is made by a Taylor series expansion of the macrolevel strain M ε and macrolevel strain gradient strain M ∇ ⊗ ∆ε as shown in Eq. (35), where x is the microlevel spatial coordinate vector, while r represents the microfluctuation field, i.e., the microlevel contribution to the RVE displacement field, as stated in Lesičar et al. [Lesičar, Tonković and Sorić (2014)]. Eqs.
(36) and (37) represent the microlevel strain m ε and strain gradient tensor m ∇ ⊗ m ε , respectively, which can be easily obtained by applying the microlevel gradient operator ∇ m on Eq. (35). The integral constraints given by Eqs. (38) and (39) have to be satisfied by the appropriate boundary conditions, which are chosen in this contribution as periodic, where the equality of the microfluctuation fields on the opposite sides of the RVE is prescribed. The displacement values are directly prescribed in the RVE corner nodes by using Eq. (40), while the displacements of the opposite RVE sides are connected via periodicity Eqs. (41) and (42). In order to completely define the microstructural BVP, the microfluctuation integrals represented by Eqs. (43) and (44) have to be prescribed, in which the integration is performed by trapezoidal rule. i D , 1i H and 2i H represent the socalled coordinate matrices. Detailed derivation of the coordinate matrices can be found in Lesičar [Lesičar (2015)].
Microlevel strain tensor: Microlevel strain gradient tensor: Constraints on the microfluctuation field: RVE displacement of the ith node on the RVE boundary: Periodicity equations: Microfluctuation integrals: The most important relations and variables related to the homogenization process are given in Tab. 2. The derivation of the consistent micro-to-macro transition relations for the nonlocal multiscale scheme starts with the Hill-Mandel energy condition, which correlates the work variation done at the RVE level with the work variation at the macrostructural material point, in form of Eq. (48). After some straightforward mathematical manipulations in Eq. (48), the relations for the homogenized stress and double stress tensors can be obtained, shown in Eqs. (49) and (50). The constitutive relations given by Eqs. (51)-(53) contain nine constitutive matrices, which take into account the contribution of all the heterogeneities and the interactions taking place at the microstructural level. Since the constitutive tangents are dependent on the microstructure, a proper regularization of the damage effects at the macrolevel is expected. The tangent stiffness matrices shown in Eq. (54) are formulated by the employment of only external RVE boundary nodes through the condensed stiffness matrix bb K  , derived in Lesičar et al. [Lesičar, Tonković and Sorić (2014)]. The condensed RVE stiffness matrix is constructed from the stiffness response of the whole RVE, which is, on the other hand, composed of the finite element stiffness matrices, obtained here by Eqs. (30)-(32). As obvious from the given relations, the finite element stiffness matrices are dependent on the damage variable D, which implies that the condensed RVE stiffness matrix is implicitly dependent on the damage variable, i.e.,

( )
Generally, in the second-order computational homogenization procedure, the macrostructural stiffness is obtained as a function of the condensed stiffness matrix, which can be written in the form of ( ) By insertion of Eq. (45) in Eq. (46), the following relation holds from which it is clear that the homogenized macrostructural stiffness is directly influenced by the damage growth observed at the microstructural level. Due to this degradation effect, the capability of the proposed multiscale scheme in capturing the macrostructural localization of the deformation is expected.

RVE failure conditions
With the formation of sharp localization zone inside an RVE, as depicted in Fig. 12, the macrolevel material point represented by that particular RVE should not be able to carry the load anymore. The presented numerical model considers that material remains the part of the continuum, where in reality cracks should form, which is done by keeping the stiffness values close to zero by the application of the isotropic damage law. When the homogenization of the tangent stiffness is performed over an RVE, both the points inside and outside the localization zone are included in the calculation. When the localization zone is formed, there is practically no contribution to homogenized stiffness from this area, since the damage is maximum there. On the other hand, the material points excluded from the localization zone still possess an intact or a slightly degraded stiffness, which makes a significant contribution to the homogenized RVE stiffness tangents. In order to prevent this spurious contribution, the conditions used for the detection of the occurrence of fully formed localization zone are developed here.

Figure 12: Localization zone inside an RVE
First condition that has to be fulfilled is given as where d A represents the localization area, i.e., area of all integration points which experienced the damage growth. If the ratio of the localization area from the last converged increment where α stands for a threshold value, usually taken very close to the critical value of damage. Two described conditions are necessary for the evaluation of the fully formed localization zone, since one without the other could provide the misleading information in some loading cases.

Numerical examples 4.1 Plate subjected to tensile load, homogeneous microstructure
In the following example plate with an imperfect zone subjected to tensile load, already analysed in Putar et al. [Putar, Sorić, Lesičar et al. (2017)] for numerical testing of the one-scale damage model, is utilized for the damage analysis by employing the nonlocal multiscale scheme, whereby the homogeneous microstructure is considered. By elimination of all effects resulting from the microstructural heterogeneities, clear conclusions can be made regarding the macrolevel structural response. The computational model of the problem is shown in Fig. 13, while the two different discretizations for the C 1 continuous macrolevel are given in Fig. 14. In order to show the shortcomings of employing the C 0 continuous macrolevel discretization, three additional discretizations are considered where second-order rectangular finite elements with reduced integration provided by ABAQUS/Standard are used, as depicted in Fig. 15. Regarding the scale transition relations between the macro-and microstructural levels in the C 0 -C 1 multiscale scheme, they become much simpler once the higher-order variables at the macrolevel vanish. Only the macrostructural deformation tensor is involved in the formulation of the RVE boundary conditions, which in turn cannot induce complex RVE deformation modes. Hence, the prediction of the microstructural damage evolution might not be as accurate as in the case of the C 1 continuous macrolevel discretization. The major deficiency of the model employing the C 0 continuous macrolevel discretization is the inability to include the nonlocal effects at the macrostructural level, since only the classical tangent stiffness and Cauchy stress are present in the constitutive model. The detailed derivation of all scale transition relations for the C 0 -C 1 multiscale scheme can be found in Putar [Putar (2018)]. The homogeneous microlevel is considered as a rectangular MVE with uniform C 1 continuity discretization, as shown in Fig. 16. The periodic boundary conditions are considered at the MVE boundaries during the calculation. In the example, Mazars' equivalent strain measure defined in Mazars and Pijaudier-Cabot [Mazars and Pijaudier-Cabot (1989)] as is used together with the damage evolution governed by the linear softening law given in Peerlings [Peerlings (1999)] as In Eq. . The horizontal displacement of 0.0325 u = mm is prescribed at the right edge of the coarse scale model. In order to trigger the localization, the Young's modulus is reduced by 10% in the 10 mm wide hatched area of the plate. Along the vertical edges the second-order derivatives 1,11 u , 1,22 u , 1,12 u and 2,12 u are suppressed. The first-order derivatives 1,2 u and 2,1 u , are also set to zero. These boundary conditions yield the straight vertical edges. Here, the indices 1 and 2 refer to the Cartesian coordinates x and y, respectively. In the following multiscale analyses, the microstructural nonlocal parameter of micro 0.6608 l = mm is considered at the MVE level, while the MVE size is taken as 2.6 L = mm. The same numerical example is then calculated by the one-scale model based on the modified case of Mindlin's form II strain energy density described in [Altan and Aifantis (1992); Ru and Aifantis (1993)], while an adequate microstructural parameter according to the consideration in ] is taken into account. Herein, 1 l = mm is used. The structural responses for the three different macroscale discretizations shown in Fig. 15 and used in the C 0 -C 1 multiscale scheme, together with the structural response obtained by the converged one-scale damage model are depicted in Fig. 17. Obviously, three different mesh densities provide three different structural responses for C 0 -C 1 multiscale scheme. Increase in the mesh density, i.e., decrease in the finite element width in the damage process zone, where energy dissipation takes place, leads to more brittle response. Such behaviour can be attributed to the local continuum model where the strain tends to localize in the smallest possible volume, or a narrowest band of finite elements in this case. With this example it is proven that it is not possible to obtain the objective results in multiscale damage analysis by the employment of C 0 continuity at macrolevel. The reaction force diagrams obtained by the two multiscale analyses with the macroscale discretizations depicted in Fig. 14 and by the one-scale damage model are given in Fig.  18. As evident, the results obtained by the multiscale analyses are now basically identical, which confirms that, generally, in the multiscale problems which involve the damage analysis at the microlevel, the macrolevel has to be discretized by the numerical scheme which enables the regularization of the strain localization phenomenon. A slight deviation from the results obtained by the one-scale damage model can be ascribed to the treatment of the nonlocality. While in the one-scale model the initial nonlocal material behaviour is degraded by the isotropic damage law as shown in Eqs. (20) and (21), in multiscale analyses the nonlocality continuously changes with the evolution of the microstructure. It can be seen that the reaction forces are closest at the onset of softening and then start to deviate as the nonlinearity progresses.

Plate subjected to tensile load, heterogeneous microstructure
Herein, the same example analysed by both C 0 -C 1 and C 1 -C 1 multiscale procedures and the consideration of the homogeneous microstructure, is analysed for the heterogeneous microstructure, included via the RVE presented in Fig. 19. In order to get intense localization zones, the microstructural internal length scale is taken as micro 0.025 l = mm, while the size of the RVE is left 2.6 L = mm. The computational model of the macroscale analysis remains the same as shown in Fig. 13, and all material parameters and softening characteristics are kept the same as described previously. The discretization of the macrostructural model is made by employing 48 triangular finite elements, as it is depicted in Fig. 14(a). Considering that such finite element mesh density provided converged results in case of the homogenous microstructure, as can be seen from Fig. 18, it can be assumed that accurate solutions will also be reached in case of the heterogeneous microstructure. Besides, a finer discretization would lead to a significant slowing down of the already very time-consuming computational process. This is a consequence mainly of the computationally very expensive matrix inverse calculations needed for the formation of the condensed stiffness matrix bb K  , as explained in Lesičar et al. [Lesičar, Tonković and Sorić (2014)]. The load-displacement curve obtained by the presented model is depicted in Fig. 20, where also the previously obtained structural response for the homogeneous microstructure is plotted for the comparison purpose. Obviously, the softening is initiated at a smaller reaction force level for the heterogeneous material. After the peak is reached, a very steep drop of the reaction force is noticed. In the peak stage of the analysis, the area of the localization zone is already formed at the macrostructural level, and few subsequent converged incremental steps lead only to the rise of the deformation level in the centre of the zone. It should be mentioned that, in order to cross the peak, the convergence criteria have to be relaxed slightly for the ABAQUS/Standard solver. However, the continuation of the analysis is not possible anymore at some point after the peak when the problems with the convergence increase, as shown in Fig. 20.

Figure 20:
Structural responses of the plate subjected to tensile load obtained by the C 1 -C 1 multiscale scheme for homogeneous and heterogeneous microstructure The deformed shape of the plate at the final stage of the analysis with the distribution of the strain tensor component in the direction of the x axis is depicted in Fig. 21. Due to very sparse mesh, smooth visualization of the variables at the macrolevel is difficult to obtain. The physically acceptable localization can clearly be seen from the figure, which indicates that for the heterogeneous microstructure the damage initiation and subsequent position of the macrocrack can successfully be captured by the proposed nonlocal multiscale algorithm. The distribution of the damage over several characteristic RVEs in the macrostructural localization zone is presented in Fig. 22. It can be observed that the most intense damage bands are formed in the middle of the plate, where the localization is the strongest. By moving away from the localization, the damage bands at the microlevel are becoming milder, until they eventually become negligible and disappear. Obviously, the material behaviour can be interpreted as realistic at both scales.

Shear band problem
In the final example C 1 -C 1 multiscale scheme is employed for the analysis of the shear band problem, as considered in Putar et al. [Putar, Sorić, Lesičar et al. (2017)] for numerical testing of the one-scale damage model. The computational model for the analysis of the heterogeneous macrostructure, presenting a quadrilateral plate with the dimension h = 60 mm, is given in Fig. 23 Fig. 23(a). Since the loaded edges have to remain straight during the analysis, the boundary conditions for the straight edge are enforced there by supressing the degrees of freedom 2,11 u , 2,22 u , 1,12 u , 2,12 u , 1,2 u and 2,1 u .
The finite element mesh of the whole model is depicted in Fig. 23(b). Homogenous microstructure is modelled by the MVE of side length 2.6 L = mm, as presented in Fig.  16, same as in the previous example. On the other hand, for the modelling of heterogeneous microstructure the RVE given in Fig. 19 is used, where the side length of 2.6 L = mm and the microstructural size of the nonlocal interactions micro 0.025 l = mm are considered. The load-displacement curves for both the homogeneous and the heterogeneous microstructure are given in Fig. 24. The softening initiation happens at a smaller reaction force level for the heterogeneous material, while the subsequent drop, although pronounced, is very short due to emergence of convergence problems. The distribution of the equivalent elastic strain at the onset of the softening for the heterogeneous material can be seen in Fig. 25. Obviously, the localization zones have started to develop as expected for the shear band problem. The distribution of the damage over several characteristic RVEs in the macrostructural localization zone at the onset of the softening is presented in Fig. 26. The direction in which the damage has evolved so far suggests that the lower shear band would eventually become dominant. The similar situation is already predicted by using the one-scale damage model where the microstructural evolution is excluded from the analysis, as shown in Putar [Putar (2018)]. As evident from the results of the considered numerical examples, C 1 -C 1 multiscale scheme is able to successfully predict the initiation of the localization at the macrostructural level. By taking into account that the constitutive behaviour of the macrostructure is obtained directly from the analysis of the evolving heterogeneous microstructure, it can be said that such algorithm represents a step forward with respect to the one-scale damage model presented in Putar et al. [Putar, Sorić, Lesičar et al. (2017)]. Although the structural responses can be qualitatively interpreted as physically acceptable, the additional comparison with the experimental results is desirable. It should also be stressed that, once the formation of the initial localization zone at the macrolevel is reached, convergence issues can be noticed in the computational process. It is not yet clear why this happens, and further research is needed to reach some reasonable conclusions.  It is to note that the computation performed is time consuming. To increase the computational efficiency, the parallelization of the macroscale computations by computing several elements at the same time has been used, which is provided within the Abaqus software infrastructure. An additional parallelization inside of the macroscale element has been applied too, where the RVEs appointed to material points of a single macrolevel element are computed simultaneously. In this manner, the multiscale computation of one complete macrolevel element using the RVE depicted in Fig. 19 can be finished relatively fast. Due to the linear elastic material employed, the homogenization of the constitutive behaviour can be done only once as a pre-processing step, since it remains constant until development of the damage, which also contributes to the faster multiscale computation. Hence, prior to damage initiation inside of the RVE, the multiscale computation can be significantly boosted. However, after damage initiation at specific macrolevel material points, the complete homogenization of the RVE behaviour is required due to the nonlinearities involved, which results in the decrease of the macroscale loading increment and increased computational time.

Conclusion
The presented paper is concerned with the multiscale modelling of materials exhibiting the strain localization phenomena, whereby the evolution of damage is observed at the microstructural level described by the appropriate representative volume element. In order to simplify the problem at the microscale, the theory based on the modified case of Mindlin's form II strain energy density is employed, where only a classical constitutive tensor and the internal length scale parameter are required in constitutive relations. The triangular C 1 finite element for softening analysis is derived in the similar manner as presented in Putar et al. [Putar, Sorić, Lesičar et al. (2017)], by the application of the isotropic damage law to the constitutive relations of the gradient elasticity theory. This finite element formulation employing microstructural length scale parameter and damage variable is applied for the modeling of damage evolution at the microstructural level of heterogeneous materials for the first time. Applicability of the finite element is tested on the RVE example, where several loading conditions are applied, and the results show the expected material behaviour. It is found that the change in the size of the nonlocal interaction zone leads to the formation of completely different localization paths. Implementation of the derived finite element at the microstructural level of the multiscale model is carried out by the user defined subroutine UEL in scope of the FE software ABAQUS/Standard. In order to test the applicability of both C 0 and C 1 continuous discretization at the macroscale, a standard benchmark example with different finite element meshes is considered, where a homogeneous microstructure is analysed. The objectivity of the novel finite element formulation has been discussed and its advantage over the standard finite element procedure is proved. As expected, the results obtained by the local-to-nonlocal multiscale scheme are dependent on the discretization and the convergence to physically acceptable solution cannot be reached. In contrast, the nonlocal multiscale scheme shows almost identical structural responses for different finite element mesh densities, and the results are comparable with the one-scale damage model solution. In order to assess the complete formation of the localization zone across the RVE, and to emulate the formation of the crack in the macrostructural integration point, the so-called RVE failure conditions are derived. Two examples including heterogeneous microstructure are considered, where it is shown that the presented multiscale damage model can successfully describe the initiation of the macrostructural localization.