INTRODUCTION

In quasibrittle materials such as concrete, rocks, soils, ice, tough ceramics, or certain composites, failure is preceded by a gradual development of a nonlinear fracture process zone and by localization of strain. Realistic failure analysis of quasibrittle structures requires the consideration of progressive damage due to microcracking, which is usually modeled by a constitutive law with strain softening. In the context of a standard Boltzmann continuum, the boundary value problem becomes ill-posed due to the loss of ellipticity of the governing differential equation. Numerically obtained results are not objective with respect to the discretization, and the total energy consumed by the fracture process tends to zero as the computational grid is refined.

Simple remedies are based on the adjustement of the softening part of the stress-strain law depending on the discretization parameters, e.g., on the size, type and orientation of finite elements and on the number of integration points. The aim of these techniques, known as the crack band approach, composite fracture model, or mesh-adjusted softening modulus, is to enforce the correct energy dissipation by keeping the product of the width of the computationally resolved process zone and the local dissipation density constant. In many cases, this leads to globally realistic responses, but no information about the strain and damage distribution within the fracture process zone is obtained. Moreover, the expression for the softening modulus is usually derived for the idealized case of failure under uniaxial tension, and its validity under general triaxial stress states is questionable.

Advanced regularization methods introduce an additional material parameter -- the characteristic length, which is related to the size and spacing of major inhomogeneities and controls the width of the numerically resolved fracture process zone. Such methods can be based on nonlocal averaging, higher-order gradient theories, continua with microstructure, or viscous regularization. They serve as localization limiters, in the sense that they impose a nonzero limit on the minimal width of the process zone and prevent strain localization into a set of measure zero.

In this contribution, attention is focused on nonlocal continuum models of the integral type. Such models replace a suitably chosen local quantity by its nonlocal counterpart, obtained as a weighted average over a certain finite neighborhood of each material point, called the domain of influence. The size of this domain and the shape of the nonlocal weight function determine the characteristic length imposed by the models, which should be related to the material microstructure. The nonlocal quantity is then inserted into the original constitutive equation. Besides acting as a localization limiter, the nonlocal approach reduces the mesh-induced directional bias.

The papers by Jirasek and Patzak [4,5] present a general methodology of the derivation of tangent stiffness matrices for nonlocal damage models, based on the consistent linearization of the incremental stress-strain relations. Due to the nonlocal interactions, the stiffness matrix is nonsymmetric and its bandwidth increases compared to the ``local'' stiffness matrices. The nonlocal character of the constitutive law makes the element stiffness matrix depend not only on the degrees of freedom associated with the given element but also on the degrees of freedom associated with the nodes of all the elements whose integration points are within the radius of interaction. This leads to nonstandard contributions that must be taken into account in the assembly procedure. From the computational point of view it is therefore essential to use an averaging function with a limited support, in order to keep the global stiffness matrix sparse.

The performance of the developed procedure is evaluated and compared to the existing approaches that use the elastic or secant stiffnesses. Examples in two- and three-dimensions are presented. It is shown that the procedure is robust and the equilibrium iterations converge at a quadratic rate, which confirms that the linearization is indeed consistent. In terms of efficiency, the proposed technique is quite competitive, especially when higher accuracy is required. For larger problems, especially in 3D, the efficiency can be increased by combining the consistent linearization with an iterative solver of the system of linearized equilibrium equations. The best results have been obtained using the Generalized Minimum Residual Method solver (GMRES, see, for example, [11]) with preconditioning based on the incomplete LU decomposition.