In this paper, a recently developed higher-order microplane (HOM) model for softening and localization is implemented within a isogeometric finite-element framework. The HOM model was derived directly from a three-dimensional discrete particle model, and it was shown to be associated with a high-order continuum characterized by independent rotation and displacement fields. Furthermore, the HOM model possesses 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; the second associated with the size of these flaws and related to the micropolar character of the continuum. The displacement-based finite element implementation of this type of continua requires C1 continuity both within the elements and at the element boundaries. This motivated the implementation of the concept of isogeometric analysis which ensures a higher degree of smoothness and continuity. Nonuniform rational B-splines (NURBS) based isogeometric elements are implemented in a 3D setting, with both displacement and rotational degrees-of-freedom at each control point. The performed numerical analyses demonstrate the effectiveness of the proposed HOM model implementation to ensure optimal convergence in both elastic and softening regime. Furthermore, the proposed approach allows the natural formulation of a localization limiter able to prevent strain localization and spurious mesh sensitivity known to be pathological issues for typical local strain-softening constitutive equations.

## Introduction

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 [911] 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 [1216] 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 C1 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 C1 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, C0 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 [3133]. 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 C0 finite element.

## Higher-Order Microplane Model

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 [4143].

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 [4547] 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.

Following Cusatis and Zhou [50], one can calculate microplane strains and curvatures as
(1)

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; ui, $φi$ are displacement and rotation, respectively, at a generic position xi (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 r0, assumed to be a material property, can be interpreted to be the average half spacing of weak spots in the material internal structure [50].

By using appropriate vectorial constitutive equations, the microplane stresses and microplane couple stresses can be calculated from the microplane strains and microplane curvatures: formally, one can write $σβ=Fβ(εN,…,χN,…)$ and $μβ=Gβ(εN,…,χN,…)$, which can be used to obtain stress, couple stress, and high order stress tensors as
$σij=34π∫Γ1σβPijβdS; μij=34π∫Γ1μβPijβdS; Σijk=3r04π∫Γ1σβPijkβdS$
(2)

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

Furthermore, stress, couple stress, and high order stress tensors must satisfy the conservation of linear and angular momenta which reads
$πji,j+bi=0; μji,j+εijkπjk=0$
(3)

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

The weak form of equilibrium that is used for the finite element formulation read
$∫Ω(σijδγij+ΣijkδΓijk+μijδκij)dV=∫ΩbiδuidV+∫∂Ω(miδφi+tiδui+hiDδui)dS+∮CciδuidL$
(4)

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 ni 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, rj = εpqjspnq, and si 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].

## Isogeometric Finite-Element Implementation

### Isogeometric Description of Geometry.

Isogeometric analysis (IGA), introduced by Hughes et al. [51], combines computer-aided design (CAD) and finite-element method (FEM) technologies. In IGA, the basis functions used in CAD to describe exactly the geometry of a certain volume of material are also used to discretize field variables describing, for example, its deformation and mechanical behavior. CAD technologies commonly employ non-uniform rational B-splines (NURBS) basis functions which are built from B-splines [52]. B-splines are piecewise polynomials composed of a linear combination of basis functions defined recursively as
$Nip(ξ)=ξ−ξiξi+p−ξiNip−1(ξ)+ξi+p+1−ξξi+p+1−ξi+1Ni+1p−1(ξ)$
(5)

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].

B-splines cannot represent some common geometrical shapes, such as circular and ellipsoidal sections, exactly. This can be achieved by using non-uniform rational B-splines (NURBS), which can be defined in terms of B-spline basis functions as follows:
$Rip(ξ)=Nip(ξ)∑j=1nNjp(ξ)wj$
(6)

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

A NURBS curve can in general be defined as linear combinations of the functions defined in Eq. (6). Hence, in 3D, the generic position vector x=[x1x2x3]T can be expressed as
$x(ξ,η,ζ)=∑i=1n∑j=1m∑k=1lRijkpqr(ξ,η,ζ)xijk=∑I=1NcpRIpqr(ξ,η,ζ)xI$
(7)
where
$RIpqr(ξ,η,ζ)=Rijkpqr(ξ,η,ζ)=Nip(ξ)Mjq(η)Pkr(ζ)wijk∑î=1n∑ĵ=1m∑k̂=1lNîp(ξ)Mĵq(η)Pk̂r(ζ)wîĵk̂$
(8)

$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; Ncp = m × n × l is the total number of control points.

### Implementation of Isogeometric Element.

In this study, a isogeometric finite element is implemented in 3D with interpolation of displacement and rotational fields assumed to be independent. This element features six degrees of freedom (three displacement and three rotations) at each control point. Following the isoparametric concept, the displacement field u and rotation field $φ$ are formulated using the same NURBS basis functions of the geometry and, for a given element, they can be approximated in terms of element's control point displacements and rotations as follows:
$u=∑I=1NcpRIuIe Φ=∑I=1NcpRIϕIe$
(9)
where Ncp is the number of control points that are supported by one element and $ue=[uxi,uyi,uzi]$ denotes the nodal displacement vector and $ϕe=[φxi,φyi,φzi]$ the nodal rotation vector. It is worth nothing that in this study the adopted shape functions in Eq. (9) are quadratic which can be obtained from Eq. (8) with p = q = r = 2 and Ncp = 27 control points per element. In Eq. (9) and throughout the rest of the paper, the shape function superscripts are dropped for clarity of notation. Low- and high-order strain fields are obtained computing the spatial derivatives of the displacement field
$e=Lde=LRde=Bde$
(10)
where $e=[γij,κij,Γijk]T$ is the strain vector collecting all the strains, curvature, and high-order strain components (45 in total), d is the displacement vector, $d=[ui,ϕi]$ collecting all degrees of freedom of the element (6 × Ncp in total), L is a differential operator, R is the shape function matrix, and B is the strain–displacement matrix obtained by applying the differential operator L to the shape functions. The B matrix can be subdivided into submatrices for low-order and high-order terms and the separate contribution of different nodes
$B=[B1B2…BNcp]; BI=[BIγBIκBIΓ]T$
(11)
in which
$BIγ=[RI,x000000RI,y000000RI,z0000RI,x000−100RI,y−100RI,z000−10RI,y000010RI,z010000RI,x010], BIκ=[000RI,x000000RI,y000000RI,z0000RI,x000000RI,y000RI,z00000RI,y000000RI,z000000RI,x]$
(12)

$BIΓ=[RI,xx0000RI,zxRI,yx00RI,xy0000RI,zyRI,yy00RI,xz0000RI,zzRI,yz000RI,yx0RI,xx000RI,zx00RI,yy0RI,xy000RI,zy00RI,yz0RI,xz000RI,zz000RI,zx0RI,yx000RI,xx00RI,zy0RI,yy000RI,xy00RI,zz0RI,yz000RI,xz0000−RI,x00RI,x00000−RI,y00RI,y00000−RI,z00RI,z000000−RI,x00RI,x00000−RI,y00RI,y00000−RI,z00RI,z000−RI,x00RI,x00000−RI,y00RI,y00000−RI,z00RI,z00]T$
(13)
Once the strain displacement matrix, B, is obtained, the element stiffness matrix and the internal force vector can be calculated easily by employing appropriate quadrature rules
$K=[K11K12⋯K21K22⋮⋱], fT=[f1T,f2T…], KIJ=∫ΩBITDBJdΩ, fI=∫ΩBITsdΩ$
(14)

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 [5456], 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.

## Numerical Results

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.

### Strain Gradient Elasticity Solution of Cantilever Beam Problems.

The exact analytical solution [57] for a cantilever beam subjected to a concentrated load at the free end is used here as a reference for comparison with the strain gradient version of the high-order microplane model. To obtain the macroscopic constitutive equations used to derive the exact solution, the microplane constitutive equation for gradient elasticity are formulated as follows:
(15)
and $μN=μM=μL=0$, where , δij is the Kronecker delta, $εDγ=εNγ−εV$, and EV, ED, ET are the volumetric, deviatoric, and shear elastic moduli, respectively. As already shown by Cusatis and Zhou [50], Eq. (15) lead to classical elastic constitutive equations if $r0=0, EV=E/(1−2ν)$ and $ED=ET=E/(1+ν)$ where E = Young's modulus and ν = Poisson's ratio. With the formulation in Eq. (15), the microplane integration (see Eq. (2)) can be done analytically and explicit form of stress, couple stress, and high-order stress tensor can be obtained as follows:
$σij=EVεVVij+[EDAijkl+ET(Bijkl+Cijkl)]γkl, μij=0Σijk=r02[EVDijklmp+ET(Eijklmp+Fijklmp)]Γlmp$
(16)

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 σ11 and Σ111. 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 2r0 = 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, r0. 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 2r0 = 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 (uz) 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.

### Cosserat Solution of a Stress Concentration Problem.

The presented high-order microplane theory reduced to a Cosserat continuum when the material characteristic length, r0, is equal to zero, and the microplane constitutive equations are formulated as follows:
$σN=EVεV+EDεD; σM=ETεM; σL=ETεL$
(17)

$μN=WVχV+WDχD; μM=WTχM; χL=WTχL$
(18)
where εV is the volumetric strain, $εD=εN−εV, χV=Vijχij, χD=χN−χV$. Also in this case, the microplane integration (Eq. (2)) can be done analytically and explicit form of stress and couple stress tensors are obtained as follows:
$σij=EVεVVij+[EDAijkl+ET(Bijkl+Cijkl)]γklμij=WVχVVij+[WDAijkl+WT(Bijkl+Cijkl)]κkl$
(19)

where the tensors Vij, Aijkl, Bijkl, and Cijkl 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 K0 and K1 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 Nc = 1 the couple-stress theory is recovered, and for Nc = 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 mm3 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.

### Localization Behavior of a Softening Bar in Tension.

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].

The adopted softening constitutive equation follows the formulation proposed by Cusatis et al. [63]. The microplane couple stresses are assumed to be zero μN = μM = μL = 0. The microplane stresses are computed with damage-like constitutive equations as follows:
$σN=σεεN σM=ασεεM σL=ασεεL$
(20)

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: σ = E0ε, in which E0 is the microplane normal modulus. E0 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 Ht is the softening modulus in pure tension (ω = π/2) expressed as $Ht=2E0/(ℓt/ℓ0−1); ℓt=2E0Gt/σt2; ℓ0=2r0$; and Gt 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 rst = σs/σt is the ratio of shear strength to tensile strength.

By assuming E = 25 GPa, ν = 0.2, r0 = 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 mm2 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 r0 = 0 are also performed. As one might expect, for r0 = 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 r0 ≠ 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.

Cusatis and Zhou [50] conducted a 1D spectral analysis of localization and proposed to regularize the solution by computing the high-order stresses as
$Σijk=3r04π∫Γ1(σβ+S0ψβ)PijkβdS$
(21)
where $ψβ=r0PijkβΓijk$, and S0 is the localization limiter. For the 1D case and for linear softening, one can obtain a localization band equal to 2r0, if S0=(1 + 1/π2)Ht [50]. For exponential softening, in which the softening modulus is not constant, the localization limiter must be a function of the current value of strain. In this case, and again with reference to a 1D setting, one can write $S0=(1+1/π2)Ht exp(−Ht〈ε−εt〉)$ with εt = σt/E0. In order to use this 1D result in the microplane formulation one can write
$S0(ε,ω)=(1+1/π2)H0(ω)exp(−H0(ω)〈ε−ε0(ω)〉)$
(22)

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 r0 ≠ 0 and S0 = 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.

The stress locking phenomenon is due to the fact that the additional elastic high order stresses are added during the entire simulation, since the beginning when the material has yet to soften. This is clearly unnecessary. An alternative regularization scheme which avoids this shortcoming can be formulated by introducing the regularization term on the stress increments as opposed to total stress. In this case, Eq. (21) is substituted by the following equations:
$dΣijk=3r04π∫Γ1(dσβ+S0dψβ)Pijkβds; Σijk=∫dΣijk$
(23)

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

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 2r0. 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]).

## Conclusions

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.

## Acknowledgment

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).

### Appendix: Voronoi Integration Scheme

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$.

For the discretized sphere, the energetical equivalence, that relates microplane stresses, microplane couple stresses, and microplane high-order stresses to their macroscopic counterparts, can be written as
$∑s=1np(σijδγij+μijδκij+ΣijkδΓijk)Vs=∑s=1npR(σβsδεβs+μβsδχβs)As$
(24)
where np is the number of microplanes, R = 1, As is the area of the generic microplane, and Vs is the volume of a pyramid of unit height and base area equal to As. Since the macroscopic quantities are uniform inside the unit sphere and by using the kinematic constraint in Eq. (1), one obtains
$σij=3∑s=1npwsσβsPijs; Σijk=3ℓ0∑s=1npwsσβsPijks$
(25)

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

## References

References
1.
Bazant
,
Z. P.
, and
Pijaudier-Cabot
,
G.
,
1988
, “
Nonlocal Continuum Damage, Localization Instability and Convergence
,”
ASME J. Appl. Mech.
,
55
(
2
), pp.
287
293
.
2.
Bazant
,
Z. P.
,
Belytschko
,
T. B.
, and
Chang
,
T.-P.
,
1984
, “
Continuum Theory for Strain-Softening
,”
J. Eng. Mech.
,
110
(
12
), pp.
1666
1692
.
3.
Rodríguez-Ferran
,
A.
,
Bennett
,
T.
,
,
H.
, and
Tamayo-Mas
,
E.
,
2011
, “
A General Framework for Softening Regularisation Based on Gradient Elasticity
,”
Int. J. Solids Struct.
,
48
(
9
), pp.
1382
1394
.
4.
Chang
,
C.
,
,
H.
, and
Sluys
,
L.
,
2002
, “
Higher-Order Strain/Higher-Order Stress Gradient Models Derived From a Discrete Microstructure, With Application to Fracture
,”
Eng. Fract. Mech.
,
69
(
17
), pp.
1907
1924
.
5.
Cosserat
,
E.
, and
Cosserat
,
F.
,
1909
,
Théorie des corps Déformables
, Scientific Library A. Hermann and Sons,
Paris, France
.
6.
Toupin
,
R. A.
,
1962
, “
Elastic Materials With Couple-Stresses
,”
Arch. Ration. Mech. Anal.
,
11
(
1
), pp.
385
414
.
7.
Eringen
,
A. C.
,
1965
, “
Linear Theory of Micropolar Elasticity
,” Office of Naval Research, Washington, DC,
Technical Report No. 29
.
8.
Eringen
,
A. C.
,
1992
, “
Vistas of Nonlocal Continuum Physics
,”
Int. J. Eng. Sci.
,
30
(
10
), pp.
1551
1565
.
9.
Mindlin
,
R. D.
,
1964
, “
Micro-Structure in Linear Elasticity
,”
Arch. Ration. Mech. Anal.
,
16
(
1
), pp.
51
78
.
10.
Mindlin
,
R. D.
,
1965
, “
Second Gradient of Strain and Surface-Tension in Linear Elasticity
,”
Int. J. Solids Struct.
,
1
(
4
), pp.
417
438
.
11.
Mindlin
,
R.
, and
Eshel
,
N.
,
1968
, “
On First Strain-Gradient Theories in Linear Elasticity
,”
Int. J. Solids Struct.
,
4
(
1
), pp.
109
124
.
12.
Altan
,
S.
, and
Aifantis
,
E.
,
1992
, “
On the Structure of the Mode III Crack-Tip in Gradient Elasticity
,”
Scr. Metall. Mater.
,
26
(
2
), pp.
319
324
.
13.
Aifantis
,
E. C.
,
1992
, “
On the Role of Gradients in the Localization of Deformation and Fracture
,”
Int. J. Eng. Sci.
,
30
(
10
), pp.
1279
1299
.
14.
Altan
,
B.
, and
Aifantis
,
E.
,
1997
, “
On Some Aspects in the Special Theory of Gradient Elasticity
,”
J. Mech. Behavior Mater.
,
8
(
3
), pp.
231
282
.
15.
Aifantis
,
E.
,
1999
, “
Strain Gradient Interpretation of Size Effects
,”
Int. J. Fract.
,
95
(
1–4
), pp.
299
314
.
16.
Aifantis
,
E.
,
2003
, “
Update on a Class of Gradient Theories
,”
Mech. Mater.
,
35
(
3
), pp.
259
280
.
17.
,
H.
, and
Aifantis
,
E. C.
,
2011
, “
Gradient Elasticity in Statics and Dynamics: An Overview of Formulations, Length Scale Identification Procedures, Finite Element Implementations and New Results
,”
Int. J. Solids Struct.
,
48
(
13
), pp.
1962
1990
.
18.
Misra
,
A.
, and
Yang
,
Y.
,
2010
, “
Micromechanical Model for Cohesive Materials Based Upon Pseudo-Granular Structure
,”
Int. J. Solids Struct.
,
47
(
21
), pp.
2970
2981
.
19.
Yang
,
Y.
,
Ching
,
W.
, and
Misra
,
A.
,
2011
, “
Higher-Order Continuum Theory Applied to Fracture Simulation of Nanoscale Intergranular Glassy Film
,”
J. Nanomech. Micromech.
,
1
(
2
), pp.
60
71
.
20.
Pijaudier-Cabot
,
G.
, and
Bazant
,
Z. P.
,
1987
, “
Nonlocal Damage Theory
,”
J. Eng. Mech.
,
113
(
10
), pp.
1512
1533
.
21.
Fleck
,
N.
, and
Hutchinson
,
J.
,
1997
, “
,”
,
33
, pp.
296
361
.
22.
Fleck
,
N.
, and
Hutchinson
,
J.
,
2001
, “
A Reformulation of Strain Gradient Plasticity
,”
J. Mech. Phys. Solids
,
49
(
10
), pp.
2245
2271
.
23.
Gao
,
H.
,
Huang
,
Y.
,
Nix
,
W.
, and
Hutchinson
,
J.
,
1999
, “
Mechanism-Based Strain Gradient Plasticity—I: Theory
,”
J. Mech. Phys. Solids
,
47
(
6
), pp.
1239
1263
.
24.
Huang
,
Y.
,
Gao
,
H.
,
Nix
,
W.
, and
Hutchinson
,
J.
,
2000
, “
Mechanism-Based Strain Gradient Plasticity—II: Analysis
,”
J. Mech. Phys. Solids
,
48
(
1
), pp.
99
128
.
25.
Yun
,
G.
,
Hwang
,
K.
,
Huang
,
Y.
, and
Wu
,
P.
,
2005
, “
A Reformulation of Mechanism-Based Strain Gradient Plasticity
,”
Philos. Mag.
,
85
(
33–35
), pp.
4011
4029
.
26.
Huang
,
Y.
,
Qu
,
S.
,
Hwang
,
K.
,
Li
,
M.
, and
Gao
,
H.
,
2004
, “
A Conventional Theory of Mechanism-Based Strain Gradient Plasticity
,”
Int. J. Plast.
,
20
(
4–5
), pp.
753
782
.
27.
Jiang
,
H.
,
Huang
,
Y.
,
Zhuang
,
Z.
, and
Hwang
,
K.
,
2001
, “
Fracture in Mechanism-Based Strain Gradient Plasticity
,”
J. Mech. Phys. Solids
,
49
(
5
), pp.
979
993
.
28.
Papanicolopulos
,
S.-A.
,
Zervos
,
A.
, and
Vardoulakis
,
I.
,
2009
, “
A Three-Dimensional C1 Finite Element for Gradient Elasticity
,”
Int. J. Numer. Methods Eng.
,
77
(
10
), pp.
1396
1415
.
29.
Zervos
,
A.
,
Papanastasiou
,
P.
, and
Vardoulakis
,
I.
,
2001
, “
A Finite Element Displacement Formulation for Gradient Elastoplasticity
,”
Bifurcation and Localisation Theory in Geomechanics
,
H. B.
Mühlhaus
,
A.
Dyskin
and
E.
Pasternak
, eds.,
Balkema
, Perth, Australia, pp.
177
187
.
30.
Fischer
,
P.
,
Klassen
,
M.
,
Mergheim
,
J.
,
Steinmann
,
P.
, and
Müller
,
R.
,
2011
, “
Isogeometric Analysis of 2D Gradient Elasticity
,”
Comput. Mech.
,
47
(
3
), pp.
325
334
.
31.
Zervos
,
A.
,
Papanicolopulos
,
S.-A.
, and
Vardoulakis
,
I.
,
2009
, “
Two Finite-Element Discretizations for Gradient Elasticity
,”
J. Eng. Mech.
,
135
(
3
), pp.
203
213
.
32.
Amanatidou
,
E.
, and
Aravas
,
N.
,
2002
, “
Mixed Finite Element Formulations of Strain-Gradient Elasticity Problems
,”
Comput. Methods Appl. Mech. Eng.
,
191
(
15
), pp.
1723
1751
.
33.
Phunpeng
,
V.
, and
Baiz
,
P.
,
2015
, “
Mixed Finite Element Formulations for Strain-Gradient Elasticity Problems Using the Fenics Environment
,”
Finite Elem. Anal. Des.
,
96
, pp.
23
40
.
34.
Ru
,
C.
, and
Aifantis
,
E.
,
1993
, “
A Simple Approach to Solve Boundary-Value Problems in Gradient Elasticity
,”
Acta Mech.
,
101
(
1–4
), pp.
59
68
.
35.
Bažant
,
Z. P.
, and
Oh
,
B. H.
,
1983
, “
Microplane Model for Fracture Analysis of Concrete Structures
,” Northwestern University, Evanston, IL,
Technical Report No. ADP001715
.
36.
Bažant
,
Z. P.
, and
Oh
,
B. H.
,
1985
, “
Microplane Model for Progressive Fracture of Concrete and Rock
,”
J. Eng. Mech.
,
111
(
4
), pp.
559
582
.
37.
Caner
,
F. C.
, and
Bažant
,
Z. P.
,
2012
, “
Microplane Model M7 for Plain Concrete—I: Formulation
,”
J. Eng. Mech.
,
139
(12), pp. 1714–1723.
38.
Bažant
,
Z. P.
, and
Zi
,
G.
,
2003
, “
Microplane Constitutive Model for Porous Isotropic Rocks
,”
Int. J. Numer. Anal. Methods Geomech.
,
27
(
1
), pp.
25
47
.
39.
Brocca
,
M.
,
Bažant
,
Z. P.
, and
Daniel
,
I. M.
,
2001
, “
Microplane Model for Stiff Foams and Finite Element Analysis of Sandwich Failure by Core Indentation
,”
Int. J. Solids Struct.
,
38
(
44
), pp.
8111
8132
.
40.
Brocca
,
M.
,
Brinson
,
L. C.
, and
Bažant
,
Z. P.
,
2002
, “
Three-Dimensional Constitutive Model for Shape Memory Alloys Based on Microplane Model
,”
J. Mech. Phys. Solids
,
50
(
5
), pp.
1051
1077
.
41.
Beghini
,
A.
,
Bažant
,
Z. P.
,
Zhou
,
Y.
,
Gouirand
,
O.
, and
Caner
,
F. C.
,
2007
, “
Microplane Model M5f for Multiaxial Behavior and Fracture of Fiber-Reinforced Concrete
,”
J. Eng. Mech.
,
133
(
1
), pp.
66
75
.
42.
Cusatis
,
G.
,
Beghini
,
A.
, and
Bazant
,
Z. P.
,
2007
, “
Spectral Stiffness Microplane Model for Quasibrittle Composite Laminates—I: Theory
,”
ASME J. Appl. Mech.
,
75
(
2
), p.
021009
.
43.
Salviato
,
M.
,
Esfahani
,
S. E. A.
, and
Cusatis
,
G.
,
2015
, “
Spectral Stiffness Microplane Model for Quasibrittle Textile Composites
,”
arXiv
:1509.02501.
44.
Bažant
,
Z. P.
, and
Oh
,
B. H.
,
1983
, “
Crack Band Theory for Fracture of Concrete
,”
Matér. Constr.
,
16
(
3
), pp.
155
177
.
45.
Bažant
,
Z. P.
, and
Ozbolt
,
J.
,
1990
, “
Nonlocal Microplane Model for Fracture, Damage, and Size Effect in Structures
,”
J. Eng. Mech.
,
116
(
11
), pp.
2485
2505
.
46.
Bažant
,
Z. P.
, and
Di Luzio
,
G.
,
2004
, “
Nonlocal Microplane Model With Strain-Softening Yield Limits
,”
Int. J. Solids Struct.
,
41
(
24
), pp.
7209
7240
.
47.
Di Luzio
,
G.
,
2007
, “
A Symmetric Over-Nonlocal Microplane Model M4 for Fracture in Concrete
,”
Int. J. Solids Struct.
,
44
(
13
), pp.
4418
4441
.
48.
Kuhl
,
E.
, and
Ramm
,
E.
,
1999
, “
Simulation of Strain Localization With Gradient Enhanced Damage Models
,”
Comput. Mater. Sci.
,
16
(
1
), pp.
176
185
.
49.
Kuhl
,
E.
,
Ramm
,
E.
, and
de Borst
,
R.
,
2000
, “
Anisotropic Gradient Damage With the Microplane Model
,”
Comput. Methods Appl. Mech. Eng.
,
183
(1–2), pp.
87
103
.
50.
Cusatis
,
G.
, and
Zhou
,
X.
,
2013
, “
High-Order Microplane Theory for Quasi-Brittle Materials With Multiple Characteristic Lengths
,”
J. Eng. Mech.
,
140
(7), p. 04014046.
51.
Hughes
,
T. J.
,
Cottrell
,
J. A.
, and
Bazilevs
,
Y.
,
2005
, “
Isogeometric Analysis: CAD, Finite Elements, Nurbs, Exact Geometry and Mesh Refinement
,”
Comput. Methods Appl. Mech. Eng.
,
194
(
39
), pp.
4135
4195
.
52.
Farin
,
G. E.
,
1995
,
NURB Curves and Surfaces: From Projective Geometry to Practical Use
,
AK Peters
, Wellesley, MA.
53.
Nguyen
,
V. P.
,
Anitescu
,
C.
,
Bordas
,
S. P.
, and
Rabczuk
,
T.
,
2015
, “
Isogeometric Analysis: An Overview and Computer Implementation Aspects
,”
Math. Comput. Simul.
,
117
, pp.
89
116
.
54.
Auricchio
,
F.
,
Calabrò
,
F.
,
Hughes
,
T.
,
Reali
,
A.
, and
Sangalli
,
G.
,
2012
, “
A Simple Algorithm for Obtaining Nearly Optimal Quadrature Rules for NURBS-Based Isogeometric Analysis
,”
Comput. Methods Appl. Mech. Eng.
,
249–252
, pp.
15
27
.
55.
Schillinger
,
D.
,
Hossain
,
S. J.
, and
Hughes
,
T. J.
,
2014
, “
Reduced Bézier Element Quadrature Rules for Quadratic and Cubic Splines in Isogeometric Analysis
,”
Comput. Methods Appl. Mech. Eng.
,
277
, pp.
1
45
.
56.
Hiemstra
,
R.
,
Calabro
,
F.
,
Schillinger
,
D.
, and
Hughes
,
T.
,
2016
, “
Optimal and Reduced Quadrature Rules for Tensor Product and Hierarchically Refined Splines in Isogeometric Analysis
,”
ICES Rep.
,
16
(
11
), pp. 1–56.
57.
Papargyri-Beskou
,
S.
,
Tsepoura
,
K.
,
Polyzos
,
D.
, and
Beskos
,
D.
,
2003
, “
Bending and Stability Analysis of Gradient Elastic Beams
,”
Int. J. Solids Struct.
,
40
(
2
), pp.
385
400
.
58.
Lam
,
D. C. C.
,
Yang
,
F.
,
Chong
,
A.
,
Wang
,
J.
, and
Tong
,
P.
,
2003
, “
Experiments and Theory in Strain Gradient Elasticity
,”
J. Mech. Phys. Solids
,
51
(
8
), pp.
1477
1508
.
59.
Eringen
,
A.
,
1967
, “
Theory of Micropolar Elasticity
,” Office of Naval Research, Washington, DC,
Technical Report No. 1
.
60.
Mindlin
,
R.
,
1963
, “
Influence of Couple-Stresses on Stress Concentrations
,”
Exp. Mech.
,
3
(
1
), pp.
1
7
.
61.
Pothier
,
A.
, and
Rencis
,
J.
,
1994
, “
Three-Dimensional Finite Element Formulation for Microelastic Solids
,”
Comput. Struct.
,
51
(
1
), pp.
1
21
.
62.
Sachio
,
N.
,
Benedict
,
R.
, and
Lakes
,
R.
,
1984
, “
Finite Element Method for Orthotropic Micropolar Elasticity
,”
Int. J. Eng. Sci.
,
22
(
3
), pp.
319
330
.
63.
Cusatis
,
G.
,
Pelessone
,
D.
, and
Mencarelli
,
A.
,
2011
, “
Lattice Discrete Particle Model (LDPM) for Failure Behavior of Concrete—I: Theory
,”
Cem. Concr. Compos.
,
33
(
9
), pp.
881
890
.
64.
Caner
,
F. C.
, and
Bažant
,
Z. P.
,
2013
, “
Microplane Model M7 for Plain Concrete—II: Calibration and Verification
,”
J. Eng. Mech.
,
139
(12), pp. 1724–1735.
65.
O'Rourke
,
J.
, and
Goodman
,
J. E.
,
2004
,
Handbook of Discrete and Computational Geometry
,
CRC Press
, Boca Raton, FL.
66.
Renka
,
R. J.
,
1997
, “
Algorithm 772: Stripack: Delaunay Triangulation and Voronoi Diagram on the Surface of a Sphere
,”
ACM Trans. Math. Software (TOMS)
,
23
(
3
), pp.
416
434
.