Isogeometric Implementation of High-Order Microplane Model for the Simulation of High-Order Elasticity, Softening, and Localization

The nonlinear and fracturing behavior of quasi-brittle, heterogeneous materials such as fiber reinforced composites, toughened ceramics, and cementitious composites, among many others, is governed by weak spots and defects of their internal structures which cause initiation of cracks and localization of damage. In order to model such behavior, one can use discrete particle models or continuum models. Each type of model has its own advantages and disadvantages.

In discrete particle models, the internal structure of materials is directly simulated through its geometric approximation and accounts for the actual intrinsic randomness specific of the selected length scale. Particles replicate the effect of major heterogeneities, and their interaction is formulated by contact algorithms or through interparticle lattice struts whose constitutive behavior simulates specific physical mechanisms. As such, discrete models have the advantage that many microscopic effects can be accounted for in a straightforward manner. However, they often lead to large computational systems, especially to simulate real structures, and require a significant, often overwhelming, amount of computational resources.

Continuum models are computationally less expensive. But most of the continuum models are based on the Cauchy continuum which cannot take into account features of the internal structure and does not possess any internal characteristic length. These type of models have some well-known drawbacks. First of all, they are not able to produce a mathematically well-posed problem in the softening regime causing strong mesh sensitivity in the numerical results [1,2]. Second, the dispersive nature of wave propagation, typical of many heterogeneous materials, cannot be reproduced. Furthermore, conventional continuum formulations predict unrealistic stress singularities at the tip of sharp cracks and fail to simulate size effect phenomena correctly [3,4].

In order to overcome these drawbacks, nonlocal continuum theories have been proposed in the literature. A popular approach consists in enriching continua with high-order strain and stresses. Cosserat brothers' study [5] was one of the initial attempt in this direction. They considered independent rotational degrees of freedom, modeling smaller scale kinematics, in addition to translational degrees of freedom. Later, the couple-stress theory was developed by Toupin [6]. In this theory, the only independent field is the displacement field that also determines the rotations. Eringen developed the so-called micropolar continuum theory [7] and nonlocal elasticity [8]. Another type of high-order theory incorporate higher-order gradients of displacements, mostly second gradient (strain gradient) and sometimes higher, into the constitutive model [6]. Mindlin [9–11] proposed an enhanced elastic theory with microstructural effects by considering the potential energy density as a quadratic form of gradient of strain in addition to quadratic form of strains. Aifantis and coworkers [12–16] introduced a simplified gradient elasticity theory with one internal length scale which incorporates Laplacian of strain to enrich the continuum. They applied it for investigation of crack tips, localization of shear bands, and to simulate size effect in torsion and bending. Various forms of gradient elasticity and their performance in static and dynamic applications are discussed by Askes and Aifantis [17]. They also provided a procedure for the identification of the relevant characteristic lengths. In this work, the finite-element implementation of gradient elasticity was also discussed.

Chang et al. [4] developed a high-order gradient model based on a discrete particle structure, which included both high-order strains and stresses. They investigated the possibility of formulating high-order constitutive equations with and without high-order stresses. They found that, contrarily to the model with high-order stresses, the one without is unstable. Misra and Yang [18] developed a model for cohesive materials based on microstructural concepts. They considered a representative volume of material as a collection of grains interacting through intergranular force displacement relations. Intergranular force displacement relations are formulated based on atomistic-level particle interactions. These force–displacement relationships are then used to derive the incremental stress–strain relationship of a representative volume of material. Later, Yang et al. [19] applied this type of microstructural granular based higher-order continuum theory to model the failure behavior of nanophased ceramics. They compared results obtained from the ab initio simulations with the result of higher-order continuum theory with good agreement.

Nonlinear nonlocal constitutive equations were also proposed without enrichment of the elastic part and by using integral or gradient approaches to enrich plasticity and damage theories. Integral type models are based on replacing a local variable with its weighted spatial averaging over a certain volume of material. Pijaudier-Cabot and Bažant [20] pioneered research in this area. On the contrary, gradient plasticity or gradient damage models are based on enriching state variables with high-order terms. A phenomenological strain gradient plasticity theory was developed by Fleck and Hutchinson [21] in the context of couple-stress theory where only rotational part of strain gradients is taken into account. They extended the classical J2 plasticity to account for strain gradient effects. Later, they reformulated their model and included also stretch gradients in addition to rotational ones [22]. They investigated their model performance by applying it to the simulation of size effect on torsion of wires, sheet bending, indentation, and void growth. The so-called mechanism-based theory of strain gradient plasticity (MSG) [23,24] was proposed on the basis of a multiscale framework linking statistically stored and geometrically necessary dislocations to plastic strains and strain gradient. A reformulation of mechanism-based strain gradient plasticity was given by Ref. [25] that involves the third-order tensor of higher-order stress to a much simpler version within Fleck and Hutchinson's theoretical framework of strain gradient plasticity theory. Huang et al. [26] provided a conventional theory of mechanism-based strain gradient plasticity where the plastic strain gradient appears only in the constitutive model, and the equilibrium equations and boundary conditions are the same as the conventional continuum models. Jiang et al. [27] studied on fracture employing mechanism-based strain gradient plasticity.

Finite-element implementation of material models that includes second gradient of displacements as degrees-of-freedom, require C¹ continuity at the interelement boundaries or mixed-element formulation in which both kinematic and static variables are defined as degrees-of-freedom. Papanicolopulos et al. [28] and Zervos et al. [29] developed C¹ type finite element for gradient elasticity. Fischear et al. [30] implemented isogeometric analysis for gradient elasticity in two dimensions. On the contrary, in mixed formulations, C⁰ continuity element can be used but in this case first derivatives should be interpolated in addition to displacements which leads to formulations with a large number of degrees-of-freedom. These types of elements can be nonconforming and are required to fulfill the patch test [31–33]. Furthermore, Ru and Aifantis [34] developed a strategy based on operator split method for their simplified strain gradient model, in which the fourth-order partial differential equations of gradient elasticity are split into two separate second-order partial differential equations, and they implement it using C⁰ finite element.

Classical microplane model formulations, pioneered by Bažant and Oh [35,36], assumed that nonlinear phenomena such as, but not limited to, fracture, shearing, and plastic deformations occur on specific orientations within the internal structure of the material. Therefore, microplane models simulate these phenomena through vectorial constitutive laws formulated on the so-called “microplanes” representing generic orientations in space. The strain vector on a generic microplane is obtained by projecting the strain tensor onto a local system of reference, and the stress tensor is computed from tractions on all possible orientations in space through an energetic equivalence.

The initial formulations of microplane models focused on concrete materials and the most recent one, labeled as M7 [37], has demonstrated a remarkable ability to reproduce typical experimental data relevant to a large variety of loading conditions, from tension to confined compression; from quasi-static to highly dynamic conditions. Microplane models have also been successfully developed for other materials such as rock [38]; rigid foam and shape memory alloys [39,40]; fiber reinforced concrete [41]; as well as composite laminates [41–43].

The aforementioned microplane formulations are “local” in the sense that the stress tensor at a given point is the only function of the strain tensor at the same point. Unfortunately, such constitutive equations suffer from mesh sensitivity and spurious energy dissipation [44] when softening behavior is simulated. To overcome this problem, nonlocal microplane models were proposed by various authors [45–47] by exploiting integral methods, in which the nonlinear behavior at one point depends on a weighted average of the strain in the neighborhood of that point. These approaches while preventing mesh sensitivity quite effectively lead to a significant increase of the computational cost of simulations.

An alternative to integral microplane models is to adopt strain-gradient formulations such as the ones proposed by Kuhl et al. [48,49] and, more recently, by Cusatis and Zhou [50]. In the latter, starting from a discrete particle model, they derived a high-order microplane (HOM) theory which includes Cosserat theory and strain gradients elasticity as a special cases. The HOM model possess two characteristic lengths: the first associated with the spacing of flaws in the material internal structure and related to the gradient character of the continuum; and the second associated with the size of these flaws and related to the micropolar character of the continuum.

where $γij=uj,i−εijkφk$ is the strain tensor, ε_ijk is the Levi-Civita permutation symbol $κij=φj,i$ is the curvature tensor; $Γijk=γij,k$ is the strain gradient tensor also referred to as high order strain tensor; u_i, $φi$ are displacement and rotation, respectively, at a generic position x_i (i = 1, 2, 3) in a 3D continuum. The second order tensor $Pijβ=eiNejβ$ and the third order tensor $Pijkβ=eiNejβekN$ denote projection operators for strain and strain gradient, respectively, in which $eiβ(β=N,M,L)$ are unit vectors defining a local system of reference on each microplane (see Fig 1). The characteristic length r₀, assumed to be a material property, can be interpreted to be the average half spacing of weak spots in the material internal structure [50].

where the integrals are calculated over the surface of a unit sphere, $Γ1; β=N,M,L$ ; and summation rule applies over β.

where $πij=σij−Σijk,k$ is the so-called effective stress tensor; b_i denote body forces.

Finally, the differential boundary value problem presented above must be completed with appropriate boundary conditions. With reference to Fig. 2, one can write, on the volume boundary, $πjinj−Dj(Σijknk)+(Dlnl)hi=t¯i$ or $ui=u¯i, Σijknjnk=h¯i$ or $Dui=Du¯i, (μji−εpkiΣpkj)nj=m¯i$ or $φi=φ¯i$⁠, and $[[rinkΣpkj]]=c¯i$⁠, where n_i is the unit vector orthogonal to the boundary, $Di≡(δij−ninj)∂/∂xj$ is the surface gradient differential operator, δ_ij = Kronecker delta, $t¯i$  = applied tractions, $ui¯$ = applied displacements, $hi¯$  = applied moment energetically conjugate to $Dui, Dui=ui,jnj$ is the displacement normal derivative, $Du¯i$  = applied displacement normal derivative, $m¯i$  = applied moments energetically conjugate to rotations, $φi¯$  = applied rotations, $c¯i$  = applied edge traction acting along a sharp edge, $[[·]]$ is the jump operator defined on a sharp edge, r_j = ε_pqjs_pn_q, and s_i is the unit vector tangent to a sharp edge (See Fig. 2). Detailed derivation of Eq. (3) and boundary conditions can be found in Cusatis and Zhou [50].

where i = 1,…, n; n and p are the number and maximum order, respectively, of the basis functions. The knots vector, $X={ξ1,ξ2,…,ξn+p+1}$ is a set of nondecreasing real numbers called knots and representing coordinates in the parameter domain. If all knots are equally spaced, the knot vector is called uniform, otherwise they are nonuniform. For p = 0, $Ni0(ξ)=1$ if ξ is contained in the interval [ξ_i, ξ_i+1], and $Ni0(ξ)=0$ otherwise. For p = 1, B-spline basis functions are identical to basis function of classical constant-strain finite-element formulations. B-spline basis functions are linearly independent, have partition of unity property and their support is compact. However, they, in general, do not satisfy the Kronecker delta property unless coincident knots exists in X. A B-spline is said to be open if its first and last knots appear p + 1 times. Open knot vectors are interpolatory at the end knots which means they satisfies the Kronecker delta property at those knots [53].

where w_j are weighting factors that are selected depending upon the type of curve to be represented exactly [51,53].

$Nip, Mjq, Pkr$ are the basis function defined according to Eq. (5) and with reference to given knot vectors for each direction $X={ξ1,ξ2,…,ξn+p+1}, Y={η1,η2,…,ηm+q+1}$ and $Z={ζ1,ζ2,…,ζk+r+1}; xI≡xijk$ is the position vector of a control point I (a.k.a. node, per the usual computational mechanics terminology) defined on a 3D grid (i = 1,…,n, j = 1,…,m, k = 1,…,l) associated to the knot vectors X, Y, and Z; N_cp = m × n × l is the total number of control points.

where $s=[σij,μij,Σijk]T$ is the stress vector collecting all components of the stress tensor, couple stress tensor and high-order stress tensor; and D is the stiffness matrix, $s=De$⁠. Even though there are recent studies attempting to derive optimal quadrature rules for IGA [54–56], in this study, Gauss quadrature is used and number of Gauss points in each direction is assumed to be 1 plus the order of basis function in each direction, that is three for quadratic element which leads to 27 Gauss points per element. Gauss quadrature was originally developed for polynomial function, but it has been employed successfully for IGA element as well [51,53].

It is important to point out that the proposed element formulation implicitly assumes that the high order surface moments h conjugate to the normal gradient of the displacement field (see Sec. 2 and Fig. 2) are zero. This is due to the fact that only displacement and rotation degrees-of-freedom are introduced and no kinematic boundary conditions can be imposed on the normal gradient of displacements.

The performance of the proposed isogeometric implementation of the high-order microplane model is demonstrated in this section with reference to classical examples for gradient elasticity, Cosserat elasticity, and strain softening.

where $Aijkl=3∫Γ1PijDPklDds/4π, Bijkl=3∫Γ1PijMPklMds/4π, Cijkl=3∫Γ1PijLPklLds/4π, PijD=PijN−Vij, Dijklmp=3∫Γ1PijkNPlmpNds/4π, Eijklmp=3∫Γ1PijkMPlmpMds/4π$ and $Fijklmp=3∫Γ1PijkLPlmpLds/4π$⁠.

By assuming $ED=ET=E0/(1+ν), EV=E0/(1−2ν), ETΓ=0$ and $ENΓ=E$ then the gradient elasticity formulation in Beskou et al. [57] is obtained for the components σ₁₁ and Σ₁₁₁. It must be observed, however, that this microplane formulation leads to other nonzero high-order stress components which are absent in the analytical solution.

A cantilever beam, 1000 mm long with 100 × 25 mm cross section, is modeled with quadratic isogeometric solid elements. The beam is discretized with four different meshes with 10 × 1 × 1, 20 × 2 × 1, 40 × 4 × 1 and 80 × 8 × 2 elements, respectively. Essential boundary conditions are applied at both ends of the beam. E = 25000 MPa, ν = 0.2, and 2r₀ = 0.2L are assumed in the calculations. In Fig. 3(a), the relative displacement (normalized with displacement value obtained from the Cauchy solution) as a function of position along the beam is given for different meshes and compared with analytical result given by Beskou et al. [57]. It can be seen from the figure that with mesh refinement the tip displacement converges to a certain value. The overall calculated elastic curve is in general good agreement with the analytical solution, to which, however, it does not converge exactly. The main reason for this discrepancy is that in the analytical solution, they only consider one high-order term Γ_xxx = γ_xx,x, but the proposed HOM model is a 3D model and it features automatically other high order terms.

Figure 3(b) reports the obtained response for different values of the internal length scale, r₀. As expected, with the increase of the internal length scale, the beam response becomes stiffer. Also, for the internal length that tends to zero, the response converges to the classical Bernoulli beam solution.

The second beam example is relevant to an experimental study carried out by Lam et al. [58] on epoxy beams. Tests were carried out in plane strain conditions on end-point loaded cantilever beams with 20, 38, 75, and 115 μm depth and length-to-depth ratio equal to 10. The elastic properties of the epoxy are reported as E = 1.44 GPa and ν = 0.38 [58]. The microplane simulations were performed with $ETΓ=ET$ and $ENΓ=EV$⁠, and the material characteristic length was optimized according to the experimental result for the 20 μm beam. The best fit is obtained with 2r₀ = 33.2 μm. Based on these material parameters, the other beams are simulated. Three-dimensional quadratic isogeometric elements are used, and the beams are discretized with 40 × 4 × 1 elements. In order to provide plane strain condition, the lateral displacement (u_z) of nodes on external surfaces (x-y plane) are restricted. Figure 3(c) illustrates the change of normalized bending stiffness, $D=PL3/(3wbh3)$⁠, predicted by the proposed model and obtained in the experiments. In the definition of D, P is the applied load; w is the corresponding displacement; b, h, and L are beam thickness, depth, and length, respectively. The model can reproduce very well the beam behavior exhibiting size effect: the smaller the thickness the larger the normalized bending stiffness [58]. This behavior cannot be produced by classical beam, which always gives constant value for the bending stiffness.

where the tensors V_ij, A_ijkl, B_ijkl, and C_ijkl are the same ones introduced in Sec. 4.1. Cusatis and Zhou [50] derived the relationship between the microplane parameters and the ones of Cosserat elasticity as $EV=3λ+2μ+χ, ED=5μ+χ, ET=χ$ and $WV=3π1+π2+π3, WD=π2+4π3, WT=π2−π3$⁠.

In order to demonstrate the efficiency of the isogeometric element for the integration of the Cosserat theory, this section discusses a well-known stress concentration problem consisting of an infinite plate with a circular hole subjected to a uniaxial far field tension. The Cosserat continuum has a characteristic length which arises from the unit mismatch between stresses and couple stresses. This affects the stress concentration factor around the hole. Analytical solutions of this problem for classical Cosserat continuum were given by Eringen [59] and for Couple Stress theory by Mindlin [60] in 2D case. In both references, the stress concentration factor is defined as $SCF=σx(0,r)/σx∞=(3+F)/(1+F)$⁠, where $F=[8lb2/lc2(1−ν)]/[4+(r/lc)2+2(r/lc)K0(r/lc)/K1(r/lc)]$⁠, r is the radius of the hole, $lb=[π3/2(2μ+χ)]1/2, lc=[(π2+π3)/(2μ+χ)]1/2$ [61] are the material characteristic lengths, and K₀ and K₁ are the second kind modified Bessel functions of order zero and one, respectively. Mindlin gives exactly the same equation for the case of $lb2/lc2=1$⁠. Furthermore, it is convenient to define [61] the coupling factor, $Nc=[χ/2(μ+χ)]$⁠, which is a nondimensional quantity taking values between 0 and 1. For N_c = 1 the couple-stress theory is recovered, and for N_c = 0, the Cauchy continuum is obtained [61].

This example was also simulated by Pothier and Rencis [61] and Sachio et al. [62]. In the finite-element simulations discussed in this section, the material parameters are taken from Ref. [61], and they are reported in Table 1. Four simulations with different coupling factor are performed by assuming plate dimensions of 300 × 300 × 10 mm³ and hole radius of 10 mm. Considering the symmetry along the x- and y-axes, only one quarter of the specimen is modeled as shown in Fig. 4. Applied boundary conditions are also specified in Fig. 4. A 3D finite-element model was setup with quadratic isogeometric solid element, and the IGA meshes was obtained by using IGAFEM [53]. Two different meshes are used with 1024 (32 × 32) and 4096 (64 × 64) elements. For both meshes, only one element is used along the thickness. The coarser mesh is shown in Fig. 4, where a zoom-in highlights the geometry of one element and the corresponding control points. Two more layers of control points are present through the thickness (for a total of 27 control points) but they are not visible in the figure. Table 2 presents the comparison of the obtained stress concentration factors with the analytic solution. The results show very good agreement with analytical values with error less than 3.7% and 0.9% for mesh-1 and mesh-2, respectively.

This section presents the simulation of softening behavior with the objective of evaluating the effectiveness of the regularization strategy proposed by Cusatis and Zhou within the presented high-order isogeometric implementation [50].

where α is a material parameter, $σ=(σN2+σT2/α)1/2$ is the effective stress, $σT=(σM2+σL2)1/2$ is the total shear stress, $ε=(εN2+αεT2)1/2$ is the effective strain, and $εT=(εM2+εL2)1/2$ is the total shear strain. In the elastic regime, the effective stress is proportional to the effective strain: σ = E₀ε, in which E₀ is the microplane normal modulus. E₀ and α can be related to Young's modulus, E, and Poisson's ratio, ν, through the following expressions: $E0=E/(1−2ν)$ and $α=(1−4ν)/(1+ν)$⁠. The nonlinear behavior is imposed by computing the effective stress incrementally, $σ˙=E0ε˙$⁠, and satisfying the inequality $0≤σ≤σbt(ε,ω)$ where $σbt=σ0(ω)exp [−H0(ω)〈ε−ε0(ω)〉/σ0(ω)], 〈x〉=max{x,0}, ε0(ω)=σ0(ω)/E0$⁠, and $tan(ω)=εN/αεT=σNα/σT$⁠. The post peak softening modulus is defined as $H0(ω)=Ht(2ω/π)nt$⁠, where H_t is the softening modulus in pure tension (ω = π/2) expressed as $Ht=2E0/(ℓt/ℓ0−1); ℓt=2E0Gt/σt2; ℓ0=2r0$⁠; and G_t is the mesoscale fracture energy. In the previous equations, the effective strength is formulated in such a way to have a smooth transition between pure tension (ω = π/2) and pure shear (ω = 0) with parabolic variation in the σ_N versus $σT/α$ plane: $σ0(ω)=σtrst2(−sin(ω)+ sin2(ω)+4α cos2(ω)/rst2)/[2α cos2(ω)]$⁠, where r_st = σ_s/σ_t is the ratio of shear strength to tensile strength.

By assuming E = 25 GPa, ν = 0.2, r₀ = 5 mm, σ_t = 3 MPa, σ_s/σ_t = 4 MPa, and ℓ_t = 100 mm, a prismatic bar subject to uniaxial tension is modeled with 3D quadratic isogeometric solid elements. The bar is 100 mm long and has 75 × 75 mm² rectangular cross section. One end of the bar is fixed in all directions, and a displacement is prescribed at opposite end in the longitudinal direction. Four different meshes are used with 10 × 1 × 1, 20 × 1 × 1, 40 × 1 × 1 and 80 × 1 × 1 elements. The microplane stress integration is performed by a Voronoi integration scheme (see the Appendix) with 66 microplane over the entire unit sphere. In order to control the location of the localization zone, the central one tenth of the bar is weakened by reducing the tensile strength by 10%.

For comparison, simulations with r₀ = 0 are also performed. As one might expect, for r₀ = 0 the bar response is mesh dependent. The load–displacement curve is increasingly more brittle (Fig. 5(a)) as the mesh becomes more refined and damage localizes always in one single element (see in Fig. 5(e)), where the longitudinal strain profile is plotted for a bar end displacement of 0.2 mm. Figures 5(b) and 5(f) show that the mesh sensitivity is only partially mitigated for r₀ ≠ 0. In this case, the bar response is increasingly brittle upon mesh refinement immediately after the peak. The solution is less dependent on the mesh in the softening branch. The strain profile (Fig. 5(f)) does not show the extreme strain localization of the local simulation but still there is not a clear convergence upon mesh refinement.

where ω accounts for the presence of shear on the microplanes (see the adopted constitutive equations above).

From the physical point of view, this regularization strategy corresponds to adding, in parallel with the softening microplane stresses, additional nonsoftening microplane stresses which ensure that the phase velocity of a propagating uniaxial wave is always real. It is worth noting that such parallel coupling of softening and elastic components is also intrinsically captured by the microplane model by virtue of the kinematic constraint (Eq. 1) and the micro–macro stress relation (Eq. 2). This is the reason why partial regularization is attained with r₀ ≠ 0 and S₀ = 0 (Figs. 5(b) and 5(f)).

Figures 5(c) and 5(g) show the obtained response by adopting the formulation in Eqs. (21) and (22). As one can see, the initial softening branch is very well regularized. However, the rest of the softening curve shows an oscillatory behavior, the convergence upon mesh refinement is not clear, and a residual load appear even for relatively large strains (stress locking). The strain profile does not localize in one element but, again, clear convergence cannot be demonstrated.

where S₀ = 0 for ε < ε₀(ω) and Eq. (22) holds for ε ≥ ε₀(ω).

Figures 5(d) and 5(h) shows the results obtained with this incremental formulation. As one can see, both the load–displacement curve and the strain profile along the bar are fully regularized and both show a clear convergence upon mesh refinement. The localization band is larger than the theoretical 1D value of 2r₀. This is due to the intrinsic 3D formulation of the microplane model which, even under uniaxial tensile macroscopic conditions, feature tensile and shear strains and stresses at each microplane orientation. The correct estimate of the localization band requires a fully 3D spectral analysis, which, however, is intractable from the analytical point of view (see discussion in Ref. [50]).

This paper presents the isogeometric finite-element implementation of a recently developed high-order microplane (HOM) model. The HOM model was originally derived based on a discrete particle model and the resulting theory includes gradient elasticity and Cosserat theory as a special cases. In addition, the HOM model allows for an effective regularization of softening constitutive equations by means of a simple modification of the high-order stresses. The numerical simulations are carried out with an isogometric finite element characterized by 27 control points and six degrees of freedom, three displacement components and three rotation components, at each control point. The spatial integrals were performed with a gauss quadrature scheme with 27 gauss points, and the microplane integration was carried out by a novel integration scheme based on the Voronoi discretization of the unit sphere with 66 microplanes.

Based on the results presented in the paper the following conclusions can be drawn:

1. (1)

   The implemented finite element with adoption of the HOM model as constitutive equation performs very well in the numerical simulations of classical high-order elasticity problems. This was verified by simulating bending of strain-gradient elastic cantilever beams and computing the stress concentration due to a circular hole in a Cosserat elastic plate subject to tension.

2. (2)

   The HOM formulation with softening constitutive laws is only partially regularized, and it shows mesh dependence in the softening branch close to the peak load.

3. (3)

   The regularization technique previously proposed by Cusatis and Zhou [50] solves the mesh dependency in the initial post peak but leads to stress locking and lack of clear convergence in the far post peak.

4. (4)

   An optimal regularization of the postpeak response is obtained by introducing the regularizing term on the high-order stress increments as opposed to the total high-order stresses. With this approach a clear convergence can be observed for both the load versus displacement curve and the strain profile along the bar.

Financial support from the U.S. National Science Foundation (NSF) under Grant No. CMMI-1435923 is gratefully acknowledged. The work of the first author was also partially supported by the Scientific and Technological Research Council of Turkey (TUBITAK).

For inelastic constitutive equations, Eq. (2) can only be computed numerically. In the previous microplane model works [35,36,64], gaussian optimal integration formulas were developed for the integration over the unit hemisphere. However, for the high-order microplane model adopted in this paper, integration over the entire unit sphere is required. Hence, a new integration scheme was developed based on the voronoi tessellation of the unit sphere [65,66]. Figure 6 show the discretization of a unit sphere with 66 microplanes. The corresponding areas and spherical angles are reported in Table 3.

The local system of reference, attached to each microplane, can be calculated as $n1=sin ϕ cos θ, n2=sin ϕ sin θ, n3=cos ϕ;m1=cos ϕ cos θ, m2=cos ϕ sin θ, m3=−sin ϕ$⁠; and $l1=−sin θ, l2=cos θ, l3=0$⁠.

where $ws=As/∑s=1npAs$ is the area ratio of microplane s over the total area of all microplanes.