On Ductile Damage Modelling of Heterogeneous Material Using Second-Order Homogenization Approach

Jurica Sorić,Tomislav Lesičar and Zdenko Tonković

Faculty of Mechanical Engineering and Naval Architecture,University of Zagreb,Zagreb,10000,Croatia

ABSTRACT The paper deals with the numerical modelling of ductile damage responses in heterogeneous materials using the classical second-order homogenization approach.The scale transition methodology in the multiscale framework is described.The structure at the macrolevel is discretized by the triangular C1 finite elements obeying nonlocal continuum theory,while the discretization of microstructural volume element at the microscale is conducted by means of the mixed type quadrilateral finite element with the nonlocal equivalent plastic strain as an additional nodal variable.The ductile damage evolution at the microlevel is modelled by using the gradient enhanced elastoplasticity.The macrolevel softening is governed by two criterions expressed by the increase in homogenized damage variable and the threshold of the local equivalent strain.The softening at each material point at the macrolevel is detected by the critical value of the homogenized damage,where homogenization of the damage variable is performed only within softening area.Due to the nonlocal continuum theory applied,a realistic softening behaviour is demonstrated after the damage initiation,compared to the widely used first-order homogenization approach.All algorithms derived have been embedded into the finite element code ABAQUS by means of the user subroutines and verified on the standard benchmark problems.The damage evolution at both microlevel and macrolevel has been demonstrated.

KEYWORDS Ductile damage;second-order homogenization;multiscale approach; C1 finite element

1 Introduction

It is well-known that the heterogeneities appearing at the material microscale,such as porosities,inclusions and various kinds of defects may lead to undesired phenomena like damage and material softening,which can result in loss of structural load-carrying capacity as well as structural integrity.Derivation of numerically efficient and accurate algorithms connecting both microstructural and macrostructural deformation responses still remains a challenge in the scientific and engineering community.

The conventional numerical approaches employing softening constitutive behaviour are unable to adequately model the microstructural deformation[1–3].This is a direct consequence of the basic assumptions of the classical local continuum theory concerning the material point concept and homogeneity.On the other hand,from physical and numerical standpoint,the nonlocal models have shown very good regularization abilities[4].Based on the definition and implementation into constitutive model,nonlocal gradient theories can be divided into explicit and implicit.In explicit methods,the nonlocal field is defined as an extension of the corresponding local one by high-order members[5].Implicit approaches are dealing with additional differential equations relating between local and nonlocal variables which should be solved simultaneously[6].Furthermore,the most recent approach to fracture modelling is the phase-field framework[7]which seems to be very effective for handling complex fracture processes concerning with crack branching or curvilinear crack paths[8]and it is well-suited for problems such as brittle and ductile fracture[9].The variational phase-field approach to fracture can be considered as a special type of the gradient damage models where the sharp crack discontinuity is approximated by a smeared surface using the phase-field variable which continuously varies over the domain and differentiates between fully broken and intact material phases[10].

From physical standpoint,explicit modelling of the complete material microstructure,and considering of all relevant microstructural governing mechanisms is the most accurate approach for modelling of material response.Unfortunately,modelling with such high fidelity still represents burdensome task for engineering practices.Therefore,efficient computational strategies based on multiscale procedures have been developed,incorporating physical understanding of material behaviour at the lower scales[11–15]in order to describe the response of coarse scale problems.Therein,the boundary value problem(BVP)driven by macroscopic quantities has been solved over a sample of material called representative volume element(RVE).In this framework,the computational homogenization has been proven as highly versatile[16].

The formulation of the second-order homogenization relies on employment of the nonlocal continuum theory at the macrolevel,where surrounding material influences on the behaviour of a considered material point[17,18].The microstructural scale is usually described by ordinary local continuum,but the nonlocal theory can be appointed to the microscale too,as presented in the authors’previous work[19].The nonlocal theory adopted at macroscale is mostly based on the Mindlin’s theory[20,21],as it is the case in this contribution.The other approaches are also applicable,as in[22],where the Aifantis theory has been utilized.

Unfortunately,the utilization of the computational homogenization for softening materials seems to be questionable.One of the reasons is satisfaction of the scale separation principle which states that the RVE size should be much smaller than variation of the macrolevel loading[23],except when there is no clear distinction between scales or when scale separation is not applicable.For the first-order homogenization,this means that only constant strain can be imposed on the RVE boundaries[24,25].This represents a limitation when softening phenomena appears at the microscale,considering the fact that during localization,strains within macrolevel localization band are distributed in highly nonlinear manner,which should be prescribed on the RVE in constant manner(first-order scheme)or in linear manner(second-order scheme).On the other hand,it has already been mentioned that the classical continuum theory in combination with the usual numerical approaches cannot regularize the formation of strain localization.Also,it has been demonstrated in[26,27]that the RVE loses its representativeness during the localization.Therefore,in the last few years,there are many proposals how to bridge material constitutive behaviour between scales when softening is initiated.Considering constitutive behaviour,the brittle damage and linear elasticity simulations have been more fruitful compared to the modelling of ductile damage and elastoplastic material responses.For example,the multigrid methods employ overlay discretization of the microstructure over the macroscopic model at the hotspots,where damage is anticipated[28]or disassemble the coarse scale model in the interest zones[29].But,in this case the localization formation should be a priori known.An another approach introduces macroscale discontinuity enrichments,governed by the localization formation at the MVE(microstructural volume element)[30].However,this method relies on the assumption that the constitutive behaviour of surrounding materials is linear elastic only.Regarding ductile damage,the multiscale method for description of the softening in polycrystalline materials has been proposed in[31],where the integral nonlocal terms are introduced in order to preserve objectivity of the results.In[32],the ductile damage in heterogeneous materials is described by means of the mean field homogenization,with an addition of isotropization procedure.As explained in the paper,the method is suited only for mild damage.The “Failure Oriented Multiscale Formulation”for consistent upscaling of the ductile softening behaviour has been presented in[33],with an emphasis on the RVE boundary conditions during the localization.A multiscale approach which employs strong coupling between scales has been presented in[34],based on the principle of the operator split on a two-phase material.

An efficient multiscale method for modelling of damage responses at microlevel still remains an open question.A first-order homogenization multiscale scheme for modelling of ductile damage has been proposed by the authors in[35].Therein,the homogenization at the microlevel has been performed over the two microstructural samples representing two boundary value problems.One is without damage,where only elastoplastic response has been computed,and another is with the embedded nonlocal ductile damage model to compute softening evolution.Besides,the additional microstructural boundary value problems had to be solved in order to numerically compute the derivatives of damage variable with respect to three macrostrain components,which are necessary for the complete macrostructural constitutive matrix.It means,the five microstructural volume element(MVE)computations should be performed at every macrolevel material point.The standard displacement based finite elements have been used for discretization at the macrolevel,while the microstructural homogenization including damage has been carried out employing the mixed finite element formulation with nonlocal equivalent plastic strain interpolation.The benchmark examples show that the softening evolution has been captured correctly,and physically realistic structural responses have been modelled.However,it is easy to conclude that the microstructural homogenization procedure,employing the multiple MVE microscale computations,is complex and time demanding.

In another contribution by the authors[36],a two scale procedure for quasi brittle damage modelling has been shown,where the macrolevel localization,modelled by means of the nonlocal continuum theory,is governed by the microlevel damage responses.Therein,using the standard microstructural averaging procedure,an accurate and mesh independent damage evolution at the macrolevel has been obtained.On the other hand,in the case when the classical continuum theory andC0continuity finite elements have been applied at the macroscale,the wrong results have been computed.Therefore,the authors’idea is to use a standard averaging procedure with only one MVE at the microscale and to employ the nonlocal continuum theory withC1continuity finite element discretization at the macrolevel in order to simplify the computation procedure,rather than complicate ductile damage modelling,where the five MVE computations have to be performed.

Hence,in this paper a second-order homogenization procedure using standard averaging approach is employed for the computation of ductile damage in heterogeneous materials.The nonlocal theory embedded into the triangularC1finite element formulation is used for discretization at the macrolevel.The constitutive matrices are upscaled from the microlevel,where they are computed by the homogenization.Using one MVE,the implicit gradient-enhanced elastoplasticity,employing the von Mises yield function is applied for the consideration of softening behaviour.The microstructural discretization is performed by using the quadrilateral mixed finite element formulation including the nonlocal equivalent plastic strain interpolation.Therein,the homogenization procedure is much simpler than in the previous contribution[35].Instead of five MVEs,only one is considered.However,the disadvantage of the newly proposed approach is the computation of state variables for greater number of material points at macrolevel,as a consequence of using high order triangular finite element satisfyingC1continuity.Instead of 4 integration points associated to the classical quadrilateral finite elements,13 integration points are used in the finite element formulation obeyingC1continuity.It is important to mention,that theoretical parts of the paper have already been considered in the authors’previous works.Therefore,they have been presented here only shortly.The original contribution of the paper is a new methodology for computation of ductile damage evolution in heterogeneous materials.

The paper is organized as follows.In Section 2,the triangular finite element satisfyingC1continuity is shortly discussed.Section 3 gives an insight into the mixed finite element formulation employed at the microscale,ensuring objective results for the ductile softening phenomena.The second-order homogenization methodology is also briefly discussed.In Section 4 the scale transition strategy with the emphasis on the upscaling procedure is presented.Section 5 deals with the usual benchmark problems considering an academic sample of heterogeneous material,where the physically realistic results are demonstrated.In Section 6 some concluding remarks are given.

2 Macrolevel Finite Element Formulation

As explained before,the second-order homogenization methodology relies on application of the nonlocal continuum theory at the macroscale.Accordingly,the macrolevel discretization has been performed by using theC1triangular finite element which has been derived in the previous authors’contributions[19,37]and it is not repeated here.Therefore,the main expressions are only presented very briefly.The element is shown in Fig.1,adopting small strain assumption and plane strain condition.It consists of three nodes,each with twelve degrees of freedom(DOF),and the displacement field is approximated by the condensed fifth order polynomial.The nodal degrees of freedom are the displacements and their first and second spatial derivatives with respect to the Cartesian coordinates.The formulation described could be extended to the 3D formulation which requires much more complicate tetrahedralC1finite elements.

Figure 1: C1 triangular finite element

According to the nonlocal theory,the element equations are derived from the principle of virtual work which may be expressed as:

whereσandμare the stress and double stress,whileεandηare the strain and double strain tensors,respectively.u is the displacement vector,and the values t and T represent the traction and double surface traction,respectively.The integration is performed over areaAas well as local boundarys.According to Mindlin’s theory[21],in the discretization procedure the strain tensors are expressed in terms of DOF,as:

Here Bεand Bηare the strain matrices containing corresponding derivatives of the interpolation functions,and v represents the vector of DOF.The computational homogenization concept requires no assumption on the form of constitutive behaviour because the constitutive operators are computed on-the-fly during the numerical simulation.Due to the nonlinear formulation applied,all state variables should be written in the incremental form.Therefore,the following incremental constitutive relations are adopted for the updates of the stress and double stress tensors

where the material constitutive matrices Cσε,Cση,Cμεand Cμηare computed at the microlevel,using the homogenization procedure over the MVE.

Using the standard procedure in the finite element method,the linearized finite element equation KΔv=Fe−Fiis derived,where K is the element stiffness matrix,and Feand Fiare the external and internal nodal force vectors which are expressed in[37].

3 Nonlocal Ductile Damage Model and Computational Homogenization

The nonlocal implicit ductile damage model and computational homogenization have been presented separately in more detail in[6,35].Here the main expressions are presented briefly,and the connections between the governing equations of two boundary value problems are highlighted.

To describe ductile damage responses,the gradient elastoplastic formulation,performed at microlevel,employs the damage evolution law proposed in[6].Accordingly,the yield function is displayed as:

whereσerepresents the equivalent von Mises stress,whiledenotes the local equivalent plastic strain measure.The linear isotropic hardening is modelled by.The damage response is described by the following evolution law

whereβrepresents the material dependent softening parameter,andis the nonlocal equivalent plastic strain measure.

The mixed finite element formulation employed at the microlevel is based on the weak form of both the standard equilibrium equation and an additional partial differential equation of the Helmholtz type written as:

using the regularizing microstructural parameterl2,which controls the width of softening band.Furthermore,this parameter ensures objectivity of numerical results and it is related to the characteristic size of microstructural constituents.It is a fact that the microstructural parameter should be determined experimentally.However,since only academic examples are considered in this contribution,it is taken from the literature considering similar problems.As can be seen,Eq.(6)gives relation among the local and nonlocal counterpart of the equivalent plastic strain.The local counterpart determines the amount of the plastic yielding in material,while the nonlocal term determines softening.Therein,theC0continuous 4-node mixed quadrilateral finite element under plane strain condition has been derived.According to the mixed formulation,the nodal variables are the two displacement components and the nonlocal equivalent plastic strain.After the derivation procedure described in[6,35],the finite element equation may be expressed as

For the computational homogenization procedure which is in detail described in the authors’contributions[19,38],the following finite element relation is extracted from the equation system(7)

which is further partitioned in the form

Here Δfbis the nodal force vector increment comprising all boundary nodes on the MVE,while Δfavanishes in the convergence state.

According to the homogenization approach,the following condensed MVE stiffness matrix is computed from equation system(9)

which builds the constitutive relations at the macrolevel,where the incremental stress and double stress are expressed as:

Eq.(11)yields the constitutive matrices in the following form

where D and H are the coordinate matrices[37]containing the Cartesian coordinates of the all boundary MVE nodes

All formulations,presented and expressed by the relevant element matrices,have been implemented into the FE software ABAQUS via user subroutine UEL[39].It is to note that the undesired phenomenon such as volumetric locking in elastoplastic analysis does not appear because of mixed quadrilateral element formulation and a high interpolation polynomial inC1element formulation as well.

4 Scale Transition Methodology

As known in multiscale schemes,boundary value problems of multiple scales(at least two)are solved simultaneously.In this paper,the two scales are considered,where the MVE is associated to every material point at the macroscale.At the microlevel,the classical local theory concept is preserved,which governs macroscale material response,where the nonlocal continuum theory is employed.Since a single macrolevel material point is now replaced by a small part of the microstructure,all macroscale state variables are expressed as an average of the corresponding microscale values.The computational scheme of two-scale algorithm is presented in Fig.2.As evident,all multiscale simulations are run in online fashion,i.e.,during the macroscale simulation.Using the macrolevel displacement increment Δv,the displacement gradients ΔεMand ΔηMare calculated.From those displacement gradients,the increment of displacement field Δubis expressed

which is imposed over the MVE boundaries obeying the periodic boundary conditions.For the classical second-order homogenization,the averaging of Cauchy stressσMand double stressμMare performed over the complete MVE,where the well-known Hill–Mandel condition is used.Accordingly,the homogenized stress and double stress tensors are computed and may be written as:

The homogenized constitutive response is obtained by performing the static condensation procedure over the whole MVE as shown above.

Figure 2:Computational scheme of the two-scale algorithm

As written in the previous section,the averaged quantities computed at the microstructural level are expressed in terms of the damage variable,describing the softening response.It means,with the development of the localization zone on the MVE,the homogenized stress tensors and constitutive behaviour experience decrease in their values,leading to the softening at the macroscale.Considering the localized character of damage evolution,it is evident that the averaging performed in the classical homogenization manner overestimates material response.For mild damage this issue is not so apparent,but for high localization the homogenization through the complete MVE leads to the overstiff behaviour.As known in the standard homogenization,the damage variable can appear only at the microscale,contributing to the overall average material response.It means,all state variables,computed at the microlevel and upscaled using averaging procedure to the macrolevel,are expressed in terms of damage.To avoid the overestimation,a controlling parameter at the macroscale,expressed by the homogenized damage variable,is proposed.Accordingly,when the homogenized damage reaches some appropriately chosen critical value,the corresponding macrolevel material point should be excluded from subsequent computation by setting its stiffness to small magnitude,but sufficiently large to preserve numerical stability.As explained in[27],to obtain realistic homogenized material behaviour during softening,the computation of the homogenized variables should be conducted only in the MVE material points within the localization zone.The procedure for determination of the localization zone at microscale has been derived and verified in[35].The proposed method relies on the computation of the local equivalent strainεeqthrough the MVE area during loading history.The equivalent strain is taken according to[1],and it is defined as:

Herein,ε1,ε2andε3denote the principal strains.Since the plane strain is considered,ε3is set to zero.On the onset of softening,the highest equivalent strainεeqmaxis usually occurred in the localization band.However,only a single material point exhibits maximum value of equivalent strain.To capture a set of material points which form the localization band,the following threshold criterion based on the equivalent strain is introduced

This expression is used in order to make distinction whether the particular MVE material point is insideor outsidethe localization zone.In that way,the threshold parameterαcontrols the width of the band,where the homogenization of damage is performed.The verification of the proposed methodology has been thoroughly discussed in[35],where it has been concluded that setting threshold parameterαto 0.3 captures the microstructural localization band and ensures that averaged value of the damage determined within the band rises up to the failure.Besides the criterion(17),to ensure that the MVE localization band is precisely captured,an additional criterion is used which is defined in terms of the microstructural damage variableDm.Namely,the MVE localization band consists of the material points which should satisfy the criterion(17)and the increase in damage is exhibited,

In this way,the material points associated with unloading are excluded from consideration.Therefore,the MVE material points which obey conditions(17)and(18)form the localization areaVd,where the controlling parameterDis evaluated using an averaging procedure,

In the numerical examples considered in this paper,the macrolevel points in which the controlling parameterDreaches value 0.95 are excluded from subsequent computations.

5 Numerical Examples

5.1 Homogeneous Strip Subjected to Tensile Loading

The algorithm presented above is firstly tested by computing a simple homogeneous strip subjected to tensile loading,as shown in Fig.3.As evident,a weakened zone with lower yield stress is placed in the middle of the strip in order to initiate softening.

Figure 3:Macrolevel strip

Accordingly,a “homogeneous MVE” is considered which is discretized by the 4 mixed quadrilateral finite elements,as shown in Fig.4,employing the softening algorithm.

Figure 4:“MVE” of a homogeneous material

The MVE side length isL=0.5 mm.The elastic properties of material are Young’s modulus 210 GPa and Poisson’s ratio 0.3.The material exhibits linear isotropic hardening at the yield stress of 250 MPa with the hardening modulus ofh=20000 MPa.The softening is described by the damage evolution law expressed in(4)with the microstructural parameterl=0.5 mm and softening exponentβ=200.The yield stress of 240 MPa is used in the weakened zone.

Here the results obtained by the present algorithm employingC1triangular finite element and second-order homogenization are compared with the approach using the first-order homogenization procedure,where the macrolevel discretization is performed by the standardC0finite element CPE8R with eight nodes used form Abaqus element library.The first-order homogenization can be found in many references[24,40]and it is also applied in the authors publication[24].Therein,only the Cauchy stress tensor and the corresponding strain are used.It means,the classical continuum theory has been applied at macrolevel instead of the nonlocal theory considered in this paper.

Both the discretization with 70 triangular labelled as C1PE3 and the discretization with 90 quadrilateral elements CPE8R are presented in Fig.5.TheC1computational model consists of 576 DOF,while the model with quadrilateral elements has 618 DOF.To enforce the straight left and right edges,an appropriate boundary conditions are applied,as described in[19].

Figure 5:Discretizations of the strip:(a) C1 triangular finite elements,(b) C0 quadrilateral finite elements

Figure 6:Macrolevel force-displacement diagram

The macrolevel deformation response of the strip is displayed in Fig.6 for both homogenization schemes.As evident,the second-order homogenization usingC1triangular finite elements exhibits the realistic softening behaviour,which can be confirmed by comparison with the results obtained using other methods in the available literature[41–43].After reaching the peak point,continuous unloading appears,as expected.The homogenized damage distribution at macrolevel is presented in Fig.7.On contrary,the nonphysical post peak behaviour is obtained when the firstorder homogenization and the quadrilateral finite element mesh are used.The unrealistic spurious local hardening response is displayed,which is demonstrated by the increase in the macrolevel force in the post-peak behaviour.This is also associated with an unrealistic softening distribution.

Figure 7:Distribution of homogenized damage

Therein,it is shown that the second-order homogenization,based on the nonlocal continuum theory at the macrolevel,gives a realistic softening behaviour even though the standard averaging procedure over only one MVE is used in the microlevel homogenization approach displayed in Fig.2.According to the analysis performed in[36],it is expected that the finite element discretisation does not affect the results.A thoroughly mesh dependency examination is out of scope of this paper.In the following computation,some benchmark examples dealing with ductile damage of heterogeneous materials are considered.

5.2 Heterogeneous Strip Subjected to Tensile Loading

As the second example,the strip under tensile loading is again considered,but a heterogeneous material is here used.The geometry of the macromodel is presented in Fig.3.To reduce the computational cost the number of elements is decreased in the computational model.The macrostructural discretization is performed by 42 triangular finite elements using the mesh shown in Fig.8,and the boundary conditions are the same as in the previous example.The MVE is shown in Fig.9.

Figure 8:Triangular finite element mesh

The material of the MVE is an academic porous steel with the unchanged material parameters.The microstructural parameter is set tol=0.1 mm,and the damage exponent isβ=400.The MVE side length isL=0.2 mm,and it consists of 13% of voids with the average radius of 0.043 mm.The discretization is performed by 508 finite elements using the mesh as presented in Fig.9.

Figure 9:RVE of heterogeneous ductile steel

The macrolevel deformation response presented by the load-displacement diagram is displayed in Fig.10.

Figure 10:Force-displacement diagram of the strip subjected to tensile loading

As observed,after increasing of macrolevel force and reaching the peak point,the softening is exhibited which is in accordance with the solutions in previous example considering homogeneous material.Due to the porosities acting as strain concentrators,where the damage initiates and spreads through the MVE,here very fast softening occurs.It could be considered as physically correct.Unfortunately,the results cannot be compared to results from literature,because it is not possible to find the same academic MVE.An experimental study is possible,but it is out of scope of this contribution.Figs.11 and 12 present the distribution of damage variable at the micro- and macrolevel at loading stage I,where the macrolevel softening is initiated and the loading stage II at the structural collapse,as labelled in Fig.10.

Figure 11:Damage distribution at micro- and macrolevel for the loading stage I in Fig.10

Figure 12:Damage distribution at micro- and macrolevel for the loading stage II in Fig.10

According to the boundary conditions,the softening is concentrated over the weakened zone and slightly spreads to other parts of the structure,as expected.The MVEs appointed to the points A and B in Figs.11 and 12 have very similar mechanical response,but the point A demonstrates a slightly stronger softening,which is additionally confirmed by the homogenized damage distribution at the macromodel.Very similar mechanical response of the strip has been found in[35]by utilization of the multiple homogenization.However,the numerical values cannot be compared due to the differences in material properties and boundary conditions.

5.3 Plate Subjected to Compressive Loading

The last example in the paper is again the standard benchmark problem of a plate subjected to compression expressed by the displacementv.This problem is usually used for testing accuracy and numerical efficiency of algorithms proposed and may be found in extensive literature[6,44,45].The macrolevel model with discretization by 70 triangular elements is shown in Fig.13.The bottom and top edge of the plate have boundary conditions ensuring straight edge.The plate geometry parameter is set toH=50 mm.The middle plate zone with dimensions 0.1 H×0.1 H has decreased yield stress in order to initiate softening response.The material used is the same as in the previous problem.The MVE is presented in Fig.9 with the unchanged material parameters.

Figure 13:Plate subjected to compression

The force-displacement diagram presenting softening response is shown in Fig.14,where the loading levels I and II are again marked.

Figure 14:Force-displacement diagram of the plate boundaries subjected to compression

In comparison to the thin strip deformation response displayed in the previous sections,the softening in this problem is developing slightly slower.Although the material is identical in both examples,difference in the slope of the curve during the softening can be explained by the fact that in the previous example the weakened zone is spreading through entire cross-section of the macromodel.However,in this particular problem,the weakness is initiated only in a small portion of the model,and thereafter the localization propagates towards the opposite boundary.The spreading of the softening at the microscale as well as the distribution of the homogenized damage at the macrolevel for the loading stages I and II are displayed in Figs.15 and 16.

Figure 15:Damage distribution at micro- and macrolevel for the loading stage I in Fig.14

As evident,the realistic behaviour has been again obtained.The physical consistency is additionally confirmed by plotting the damage contours over the MVEs assigned to the material points A,B and C on the macromodel.The macroscale localization firstly initiates at the point A and propagates to the points B and C.The comparison of the results obtained here with the softening responses of the same sample presented in the authors’work[35],where the multiple homogenization has been used at the microlevel,can also prove accuracy of the newly proposed computational procedure.It is known that computational homogenization scheme is time demanding,as a consequence of the “online” microscale computations in every macrolevel material point.The simulations within this paper were performed on a workstation with 8 cores and clock rate 2 GHz.The numerical examples performed on the homogeneous material and presented in this paper were completed in two days.The examples comprising heterogeneous material were finished in approximately five days.Improvement of the numerical efficiency is clearly needed,through parallelization and speeding of the micro-macro scale transition procedure,as well.In the code developed by the authors,speedup of the performance was accomplished at the macrolevel,where all material points within single macrolevel element were computed in parallel.The authors’ further research is concerned with the development of much more efficient homogenization procedures using clustering approaches[46,47].

Figure 16:Damage distribution at micro- and macrolevel for the loading stage II in Fig.14

6 Conclusions

The paper demonstrates ability of the second-order homogenization procedure to model ductile damage responses in heterogeneous materials.Accordingly,theC1continuity finite element discretization employing the nonlocal continuum theory has been performed at the macrolevel.The softening behaviour is modelled at the microlevel using the implicit gradient-enhanced elastoplasticity,where the damage evolution is expressed in terms of the nonlocal equivalent plastic strain which is an additional nodal variable in the mixed quadrilateral finite element formulation applied.

The homogenization governed by the displacement gradient variables and the regularization employing the gradient-enhanced elastoplasticity,where averaging over only one microstructural volume element has been performed,yields the physically realistic softening responses at the macrolevel.However,if theC0continuity displacement based finite element discretization at macrolevel is applied,which is associated with the first-order homogenization,the unrealistic results are obtained.

To overwhelm the well-known problem of overestimating homogenized behaviour in softened materials when standard averaging scheme is employed,the homogenized damage variable computed within the MVE localization zone is used as the controlling parameter.In addition,the macrolevel softening is governed by the microlevel damage variable and the threshold criterion based on the local equivalent strain.The physically realistic damage responses of the MVE expressing by the localization bands are computed,which are connected with the realistic macrolevel softening behaviour,as expected.The results obtained compare very good with the solutions obtained in the previous authors’work,where the more complicate multiple homogenization at the microlevel has been used.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.