A thermodynamically nonlocal damage model using a surface-residual-based nonlocal stress

In this research, a surface-residual-based nonlocal stress was introduced into nonlocal damage theory to describe the long-range actions among microstructuresthatwereexcludedinthedefinitionofCauchystress.Byusingthesurface-residual-basednonlocalstresstensor,athermodynami-callyconsistentnonlocalintegraldamagemodelwasestablishedtosimulatethestrainlocalizationbehaviorforelastic-brittledamageproblems.Inthismodel,boththestrainandthedamageweretakenasnonlocalvariablesinthefreeenergyfunction,andtheintegral-typedamageconstitutiverelationshipsandtheevolutionequationwerederivedviathermodynamiclawsinordertoensuretheself-consistencywithinthethermodynamicframework.Basedonthenonlocaldamageformulationsusingarealnonlocalstressconcept,wesimulatedthestrainlocalizationphenomenoninanelasticbarsubjectedtouniaxialtension.Theresultsshowedclearlocalizingandsofteningfeaturesofstraininthedamagezone,andtheboundaryeffectsarisingfromthenonlocalsurfaceresidualwereilluminated.Furthermore,thestrainlocalizationbehaviorsfordifferentinternalcharacteristiclengthsweresimulated,throughwhichwefoundthatthecharacteristiclengthwascomparabletothesizeofthestrainlocalizationzone.


INTRODUCTION
Strain localization refers to a phenomenon in which the deformation of a solid body localizes into a narrow region (localization zone) in which a fracture first occurs. In the damage process of materials, strain localization is a macroscopic behavior closely related to microscale dissipative mechanisms (e.g. nucleation, coalescence and the development of microcracks and microdefects) [1,2] that essentially reflects the macro-micro correlation in damage mechanisms and that is of much importance to the failure prediction. However, since classical continuum theory fails to describe the micromechanical properties of a material, applying the classical continuum damage model to simulating strain localization frequently leads to some nonphysical features [3][4][5], such as zero energy dissipation when an element is refined to vanishing size and spurious mesh sensitivity. From a mathematical perspective, these incorrect phenomena can also be regarded as the result of the ill-posedness of the boundary value problem caused by strain softening [6][7][8]. To overcome the blemishes of classical damage theory, various approaches have been proposed to regularize the boundary value problem and ensure the solutions to be numerically convergent and physically meaningful [9][10][11][12][13]. Among these approaches, the most general and effective approach is the nonlocal continuum approach, which is widely used to solve macro-micro correlating problems.
Nonlocal continuum theory was developed from the 1960s [14,15], having the main feature that the local hypothesis in classical continuum mechanics was abandoned [16], which meant that the physical state of a material point in the continuum had to depend on the state of the whole body. Therefore, for the nonlocal continua, the long-range interactions among microstructures within materials and the scale effects, which are neglected in classical local theory, must be considered [17,18]. In the mid-1980s, nonlocal concepts were extended into the area of continuum damage mechanics, giving rise to the development of nonlocal damage theory. Bazant et al. first proposed an imbricate nonlocal continuum model to handle the strain softening problem in finite element analysis [3,4,19]. Then, Pijaudier-Cabot and Bazant improved this nonlocal damage model by adopting a formulation in which the nonlocal concept was only used for variables that controlled damage while the stress-strain relationship remained local [9]. In their work, the damage at a given point was measured with a spatially weighted-averaging integral with respect to the local damage distributed over a representative volume centered at that point, which was always referred to as a strongly nonlocal approach [20]. In addition, a group of researchers developed a gradient-enhanced damage model by approximating the nonlocal integral format with a truncated Taylor series expansion [21][22][23][24], which was categorized as "weakly nonlocal" [20]. So far, the integral-type and gradient-enhanced damage models have played a dominant role in nonlocal damage theory. The averaging integral or high-order gradient terms are introduced into the governing equations to account for the nonlocal effects on strain softening behaviors in physics, and to eliminate the ill-posedness of the boundary value problem in mathematics [21,25]. Additionally, the intrinsic length scale of material is incorporated into the continuum model via the weight function (nonlocal kernel) or gradient terms, acting as a "localization limiter" to limit the strain localization zone to a certain size comparable to the characteristic length [4]. By means of the nonlocal regularization techniques, strain localization solutions remain convergent when the mesh is refined, and nonzero energy dissipation is achieved [19,26]. Due to the great advantage of the nonlocal approach in solving damage problems, many works have been done to improve the nonlocal damage theory. Bazant et al. [3,4,20], Vermeer and Brinkgreve [27], Stromberg and Ristinmaa [28], Polizzotto and Borino [29][30][31], Luzio and Bazant [32] and Lu et al. [33] have made great contributions in this respect. Recently, the numerical computation of nonlocal damage model has developed rapidly, and the strain localization of brittle or ductile damaged materials could be effectively simulated through finite element methods [26,[34][35][36][37].
In existing nonlocal integral-type damage formulations, the stress-strain constitutive relationships always follow Pijaudier-Cabot and Bazant's viewpoint to remain local [3,4,9,19,20], and the measure of stress in a nonlocal damage model is still the same as the Cauchy stress, which is essentially a contact force exerted by one part of a body on another part through a common interface [15]. This treatment makes the analysis rather simple, but the long-range interactions within the microstructure of a material are excluded in the definition of Cauchy stress, which is apparently incompatible with the premises of nonlocal continuum mechanics [16,18]. As is well known, this incompatibility may cause ill-posedness in nonlocal mixed boundary value problems [38][39][40]. Therefore, to refine nonlocal damage theory, some further explorations still need to be done to define the measure of stress in nonlocal damage models.
Finding a reasonable measure of stress has long been a key issue in nonlocal continuum theory. Many works have focused on this topic [5,18,[41][42][43][44], in which the nonlocal stress was formulated as a weighted average integral of the local stress, or further simplified as the addition of local stress and its second-order gradient. Such descriptions were directly given based on the concept of nonlocality that the stress at a material point depends on the stresses at all material points over the whole body. In most of these works [5,18,41,42,44], the stress boundary condition was still the same as the classical one under local hypothesis, without considering the nonlocal surface residual [18]. It means that the stress measure in the boundary condition is actually a local one, which is inconsistent with the nonlocal stress measure in the balance equation. To solve this problem, Romano and Barretta [43] proposed a nonclassical boundary condition considering the nonlocal effect, but did not correlate their model with the nonlocal surface residual.
Recently, Huang [45] proposed a new measure of nonlocal stress based on the concept of nonlocal surface residual, in which a nonlocal stress formula depending on the Cauchy stress and the nonlocal surface residual (surface-induced traction) is derived, and the latter is consequently included in the stress boundary condition to characterize the nonlocal effect. Therefore, a nonclassical stress boundary condition related to the nonlocal surface residual was provided by Huang's work, in contrast to the classical ones used in existing nonlocal models [5,18,41,42,44]. A consistency between nonlocal stress measures in the balance equation and in the boundary condition can be guaranteed, and the boundary effect induced by nonlocality, e.g. the lattice shrinkage at the surface, can be captured with the introduction of the nonlocal surface residual. Moreover, comparing with the nonclassical stress boundary condition given by Romano and Barretta [43], the present one is simpler in mathematical derivations and gives a clear characterization of the nonlocal surface residual. This achievement provided a proper way to overcome the aforementioned shortcomings by substituting the new measure of stress for the Cauchy stress in nonlocal damage formulations. The aim of this research was to establish a nonlocal integral damage model on the basis of the surface-residualbased nonlocal stress in [45] and apply the model to investigating strain localization behavior in brittle damage problems.
The paper is divided into five parts. After the introduction in the first section, the second section briefly presents the new surface-residual-based nonlocal stress tensor. The third section describes how the nonlocal form of Clausius-Duhem inequality for elastic-brittle damage mechanisms was established with the new definition of nonlocal stress. In the inequality, the strain and the damage were taken as nonlocal variables. Then, the corresponding nonlocal integral-type constitutive equations and the damage evolution equation were derived via the thermodynamic laws. In the fourth section, we describe the application of the new damage model to simulate the strain localization phenomenon in an elastic bar subjected to uniaxial tension, and the strain localization behaviors for different internal characteristic lengths were computed in order to discuss the influence of the intrinsic length scale on the size of the localization zone. Finally, conclusions are drawn in the final section.

NEW DEFINITION OF NONLOCAL STRESS
According to Pijaudier-Cabot and Bazant's study, the stressstrain relationship for linear elastic-brittle damage cases is commonly written as [4,9] where D e is the fourth-order elastic moduli tensor, σ and ε are the local stress and the strain tensor, respectively,d is the nonlocal isotropic damage scalar obtained by the averaging integral over domain V, while d represents its local counterpart, and α(x, y) is the weighting function (nonlocal kernel). As mentioned above, the measure of stress in Eq. (1) is defined as the Cauchy stress, and the long-range forces between material points are excluded, which is contradictory to the essence of nonlocal theory.
To deal with this problem, Huang focused his study on finding a new measure of stress related to the nonlocal effects among microstructures within a material. In [45], a small tetrahedron element at point x with nonlocality was analyzed (see Fig. 1). Aside from the contact tractions p i and p n applied on the oriented surfaces A i and A, respectively, the tetrahedron was also subjected to nonlocal residual arising from the long-range interactions exerted by other parts of the body. Huang argued that the nonlocal residual could be divided into two parts: the nonlocal body residual F and the surface residual (surface-induced traction) a , with a being a long-range action induced by the cutoff of the nonlocal effects at the boundary surface, independent of the orientation of the surfaces. If the nonlocal residuals are ignored, the tetrahedron in Fig. 1 will reduce to a subvolume under Cauchy postulate, which leads to the Cauchy stress formula: σ · n = p n . However, for the nonlocal case, by taking operations analogous to the derivation of the Cauchy stress, a definition of the nonlocal stress tensor is obtained: where t refers to the nonlocal stress tensor, σ is the Cauchy stress and e i is the unit normal vector on A i . From Eq. (2), the nonlocal stress formula can be further obtained: In Eqs (2) and (3), the new measure of stress concerning nonlocality is presented, which is distinguished from the Cauchy stress by containing the long-range actions induced by the boundary surfaces, i.e. the nonlocal surface residual. Huang managed to use this surface-residual-based nonlocal stress to predict the surface effects in nonlocal elasticity, proving the validity of Eq. (3) [45]. Therefore, in this research, we developed a new nonlocal damage model using the surface-residual-based nonlocal stress instead of Cauchy stress, based on which the strain localization behavior in brittle damage problems is investigated. However, since the measure of stress was defined as nonlocal, the local stress-strain constitutive relationship [Eq. (1)] might not have been applicable in this case. In the next section, we describe the re-deduction of the constitutive equations based on the fundamental constitutive axioms and thermodynamic laws in order to ensure that the nonlocal damage formulations abide by thermodynamic principles.

ESTABLISHMENT OF NONLOCAL DA MAGE MODEL 3.1 Balance equations and thermodynamic laws
We considered a domain V occupied by the volume of a solid body. It was assumed that the continuum V was locally mass and entropy closed; i.e. the nonlocal residuals were neglected. Then, the local forms of the balance equations of the mass, linear momentum, angular momentum, energy and Clausius-Duhem inequality were written as [46] ∂ρ ∂t ρη where t is the surface-residual-based nonlocal stress tensor defined by Eq. (3), and ρ, d, f, r, ε, e, q, h, η and θ represent the mass density, rate of deformation, body force, position vector, permutation tensor, internal energy density, heat vector, heat source density, entropy density and absolute temperature, respectively. F , J and π denote the nonlocal residuals of the body force, angular momentum and energy, respectively. To satisfy the invariance of Eq. (6) with respect to the rigid motion of a material, we let F and J be zero. Moreover, we confined our research only to the case of small deformation. Therefore, Eqs (4)-(8) could be simplified into Equation (11) shows that the nonlocal stress tensor was symmetric, and the deformation rate tensor d in the energy conservation law was replaced by the rate of the strain tensorε under the condition of small strain. As seen above, Eqs (10) and (12) present the balance equations and thermodynamic laws of the nonlocal damage model using the definition of nonlocal stress, based on which we could deduce the damage constitutive relationships and the damage evolution law within a rigorous framework of continuum thermodynamics.

Constitutive relationships and damage evolution law
Then, we introduced the Helmholtz free energy density φ = e + θη. Combining this equation with Eqs (12) and (13) yields which was the nonlocal form of the Clausius-Duhem inequality. The free energy density φ was a function depending on the strain, damage and temperature, in which the damage acted as an internal variable describing the intrinsic dissipative mechanisms. Since the temperature θ was always assumed to be uniform through V for nonlocal elastic solids [18,47], we took the strain and the damage as nonlocal variables in the free energy density function, similar to the treatment in previous researches [1,9,20]. According to the nonlocal constitutive axioms [16,18], the functional forms of φ, η, t, q and π were expressed as where x and x represent the material points over the entire body V. It was assumed that φ belonged to the class of additive functionals in the sense of Friedman and Katz. Through the representation theorem of Friedman and Katz [48], we obtained Since Eq. (17) had to hold for any admissible thermal and deformation processes of the material, the state equations could be obtained: Equation (18) is the stress-strain constitutive relationship with a nonlocal integral form, which was consistent with the definition of the nonlocal stress as the long-range actions between different material points were described. As a result of Eqs (18)- (20), Eq. (17) was reduced to The left-hand side in Eq. (21) is equal to the energy dissipation rate per unit volume. The damage energy release rate Y was introduced as the thermodynamic conjugate force of the damage d, i.e.
The integral expression of Eq. (22) shows that the thermodynamic force Y associated with the damage flux d was also a nonlocal variable. Then, the damage dissipation inequality was obtained: Based on Eq. (23), it was assumed that there was a damage activation function F d (Y) that depended on the variable Y [1]. With the normality condition, the evolution law of the damage was written asḋ =λ where λ is the damage multiplier related to the damage behavior of material. It was noted that the function F d in Eq. (24) also refers to the damage complex potential satisfying the Legendre transformation below [1]: in which the function F * d (ḋ) is the damage potential corresponding to F d (Y).
For the nonlocal elasticity, the effects of the temperature on the material properties were ignored. Under the assumption of no initial stress, the free energy density ψ in Eq. (16) was set to be [18,47] where C 1 , C 2 and C 3 represent the fourth-order tensor function of two points (X, X ). Their components had the symmetry shown below: Substituting Eq. (26) into Eqs (18) and (22) led to where Commonly, it could be qualitatively assumed that the damage complex potential F d (Y) was proportional to the square of energy release rate Y [1], i.e. F d ∼ Y 2 . For simplicity, we simply let F d = Y 2 /2s, where s denotes the back stress [1]. From Eq. (26) and the expression of the energy release rate Y in Eq. (28), the damage evolution law of Eq. (24) reaḋ Up to this point, the integral-type constitutive equations and damage flow law were derived through standard thermodynamic procedures, which were self-consistent in the framework of thermodynamics. In comparison with former nonlocal damage formulations, the measurement of the stress in this newly established nonlocal damage model was the nonlocality-dependent stress tensor (see Section 2), instead of Cauchy stress that excluded the long-range forces among microstructures of the material. Accordingly, the stress-strain constitutive relationship was formulated as a nonlocal averaging integral form [see Eq. (28)], which was consistent with the definition of the nonlocal stress. Therefore, in this research, it was guaranteed that the nonlocal damage model was compatible with the basic concepts of nonlocal theory.

AN E X A MPLE OF A ONE -DIMENSIONAL BAR IN TENSION 4.1 Governing equations
As described in this section, the nonlocal damage formulations obtained in Section 3 were applied to a one-dimensional example to study the strain localization behavior of elastic-brittle materials. A homogeneous elastic bar of length L subjected to a quasi-static tensile load p at two ends x = ±L/2 was considered, as shown in Fig. 2. It was assumed that the solid bar was free of body force and inertia force. According to Eqs (3), (10) and (28), the basic governing equations for the one-dimensional damage case are presented below: where u and E represent the displacement in the x direction and Young's modulus, respectively. The coordinate x ranged from −L/2 to L/2. Since the averaging integrals in Eqs (32) and (33) were performed over a finite domain V, the nonlocal kernel α(|x − x|) was expressed as [21] α which satisfied the normalizing condition ∫ L α(|x − x|)dx = 1. α ∞ was always chosen as a nonnegative bell-shaped function.
As can be seen from Eq. (34), since the nonlocality-dependent stress was used, the nonlocal surface residual (surfaceinduced traction) a was involved in the boundary condition. Equation (3) shows that in the spatial domain, a was a continuous function on the boundary surfaces. However, for the one-dimensional case, the boundaries of the bar shrank to two end points, so a could be correspondingly taken as a constant.
In addition, to determine the concrete formulation of the damage evolution law, the damage multiplier λ had to be given. We letλ =ε(1 − d).
Substituting Eq. (36) into Eq. (33) led tȯ Equation (37) is a rate-independent form of damage evolution law, through which we could get the derivative of d with respect to ε: Integrating Eq. (38) yielded the total form of the damage evolution equation: Equation (40) could be further simplified by expanding ε(x ) into a Taylor series at x, i.e.
(41) For a small deformation case, the terms with an order higher than one in Eq. (41) were neglected. Then, Eq. (41) was substituted into Eq. (40) From Eq. (35), we had the normalizing condition: ∫ L α(|x − x|)dx = 1. It should be noted that the odd derivative term could always be eliminated approximately for the finite domain [12,22]. Consequently, Eq. (40) was linearized into Equation (43) was substituted into Eq. (39) to arrive at in which the damage solely depended on the strain. After substituting Eq. (44) into Eq. (32) and operating the same simplifying procedures as those used when deriving Eq. (43), the stressstrain constitutive relationship finally became  From Eqs. (30), (34) and (45), we obtained It was assumed that the loading increment p was uniform for every load step, so the tension exerted at the kth step was p = k p. Next, we let γ = E/s, f p = p/E, f a = a /E, x = x/L andx = x /L. Equation (46) was normalized to a dimensionless form: Equation (47) is a second Fredholm integral equation with respect to ε that has a unique solution. Therefore, the ill-posedness of the governing equations in classical damage theory could be avoided by using a nonlocal regularization technique. The next subsection describes how Eq. (47) was numerically solved to obtain the strain distribution along the bar.

Results and discussion
We used the numerical integration method to solve Eq. (47). First, we divided the one-dimensional bar into n equal segments (n was an even number), each of length h; i.e. h = L/n. Hence, the dimensionless coordinate of the ith integration node was x i = -1/2 + i/n, and the strain value atx i was ε(x i ) = ε i (i = 0, 1, 2, …, n), as shown in Fig. 3. Using Simpson's rule [49], the discretized form of Eq. (47) was generally written as  0, 1, 2, . . . , n).
where l* = l/L is the dimensionless internal characteristic length. Let γ = 1 × 10 9 , ε d = 1 × 10 −3 , f a = 3 × 10 −5 , f p = 2 × 10 −4 , l* = 0.08 and n = 500. Here, γ = E/s is a dimensionless material parameter related to the Young's modulus E and the back stress s [1]; ε d denotes the critical strain at the damage initiation; i.e. when ε ≤ ε d , d = 0; f a = a/E denotes the normalized surface residual traction; f p = p/E represents the normalized load increment; l* = l/L denotes the normalized internal characteristic length; and n is the number of segments of the bar. The strain localization behaviors at the 10th, 15th and 20th load steps are depicted in Fig. 4. It was seen that the strain was concentrated in a small region in the middle of the bar whose size was comparable to the internal characteristic length. At both sides of the localization zone, there lay two unloading areas in which the strain decreased sharply. This phenomenon could also be found in the results obtained with the over-nonlocal model [32], which might have been induced by oscillations in the numerical integration. Outside the localization and unloading zones where the deformation changed strongly, the strain distribution stayed relatively gentle. Due to the action of the nonlocal surface residual, the strain rose slightly at the two ends of the bar, which reflected the boundary effects in a solid body. However, since the nonlocal surface residual was rather small as compared to the external applied load, this perturbation of the deformation near the boundaries was insignificant.
Additionally, Fig. 4 shows that as the load increased, the strain at each point of the bar grew simultaneously with the different magnitudes. The strain in the localization zone (especially near the center) apparently rose more rapidly than that of the other parts. This was consistent with the strain evolution characteristics of elastic-brittle materials [12,26,50].
The damage localization behaviors at the 10th, 15th and 20th steps are illustrated in Fig. 5. These behaviors exhibited similar patterns to that of the strain distribution. The damage climbed rapidly to the peak value in the localization zone, beyond which it almost approached zero. With the increase of the load, the damage in the localization zone developed much more quickly than that in the nonlocalizing parts. At the final step, the maximum of the damage was close to 1, indicating that the rupture first occurred in the central part of the bar.
Next, we defined the nominal stress σ = [1 − d(ε)]Eε in order to evaluate the influence of the damage on the stress-strain responses of the elastic material. For simplicity, we let σ * = σ /E. When the load step k = 20, the σ *-ε relationships at different points in the bar were as presented in Fig. 6 (due to the symmetry of the strain distribution, only the half part -0.5 ≤ x/L ≤ 0 was considered). As can be seen in Fig. 6a-c, outside the localization zone (−0.5 ≤ x/L ≤ −0.05 roughly), the stress-strain relationships remained linearly elastic in the loading process, which indicated that the damage in this part was too small (see Fig. 5) to affect the mechanical properties of the solids, so the stress-strain constitutive relationship almost complied with Hooke's law of elasticity. However, within the localization zone, the fast growth of the damage caused severe stiffness degradation that led to the weakening of the material's resistance to deformation. As a result, when the stress rose linearly to a critical value, it started to decay sharply with the increase of the strain (see Fig. 6d-f), showing clear softening features. It was found that the softening branch of the stress-strain curve gradually increased with the approach to the center of the localization zone, and the stress value almost vanished at x = 0. This phenomenon was consistent with the damage localization profile sketched in Fig. 5.
Additionally, we made a comparison between the linear elastic stress-strain responses at different points of the bar (see Fig. 7). It was observed that the stress-strain curve at the boundary point deviated slightly from the uniform ones at interior points, which further demonstrated the boundary effect induced by the surface-residual-based nonlocal stress.
Moreover, from Figs 4 and 5, we found that the width of the localization zone always kept a constant value comparable to the internal characteristic length, no matter how the external load or the number of integration points changed. This meant that the internal characteristic length involved in the nonlocal kernel acted as a "localization limiter" to directly control the size of localization zone in practical computation. Therefore, in our integral-type nonlocal damage model, there was no need to set a fix-sized localization zone or a reduced cross-sectional area inside the bar in advance [12,29,30], which made the analysis simpler and more physical.

Influence of internal characteristic length on strain localization
Based on the nonlocal damage model in this research, we calculated the strain localization behaviors for different internal characteristic lengths to investigate the effects of varying intrinsic length scales on the size of the localization zone. Letting l* = 0.04, 0.05, 0.06, 0.07, 0.08 and 0.1, the corresponding strain distributions along the bar at the 10th, 15th and 20th steps are shown in Fig. 8. It was clearly seen that the size of localization zone was proportional to the internal characteristic length. For l* = 0.04, 0.05, 0.06, 0.07, 0.08 and 0.1, the corresponding localization widths w* were roughly 0.02, 0.032, 0.052, 0.08, 0.1 and 0.17, respectively, which had a comparable order of magnitude with l*. This relationship was consistent with the conclusions from former research [12]. For smaller l*, the strain tended to be more concentrated in the localization zone, exhibiting a larger peak value and a more prominent damage. This phenomenon might have arisen from the decrease of the intrinsic length scale with the evolution of damage mechanisms in the material, which has been reported in our another research on the damage of microscale metallic materials.
As is well known, in most nonlocal damage models, the internal characteristic length is always taken as a constant for simplicity. This treatment will lead to the unchanged size of a localization zone during the whole loading process [12,33], which is contradictory to the experimental observations [50][51][52]. Therefore, it suggests that the internal characteristic length may be changeable during the loading process, acting not only as a "limiter" to confine the domain of the localization behavior but also as a "controller" to govern the variation of the size of the localization zone. In this case, the evolution law of the localization zone could be theoretically explained and simulated by introducing a varying internal length scale into the nonlocal damage model. In the damage process, the internal characteristic length should decrease with the evolution of damage. Such variation can be described by an analytical function l = (1 − d) 2 l 0 obtained in our another research on the damage of microscale metallic materials [53], where d is the damage variable and l 0 is the internal characteristic length of an undamaged material.

CONCLUSIONS
The use of Cauchy stress in nonlocal damage theory is incompatible with the essence of nonlocal continuum mechanics. To solve this problem, we developed new nonlocal damage formulations by employing the surface-residual-based nonlocal stress tensor [45] instead of Cauchy stress, based on which the strain localization behavior of elastic-brittle materials was investigated. From the analysis, some major conclusions were drawn, as listed below.
1. The Clausius-Duhem inequality for elastic damage mechanisms was presented, in which the measure of the stress followed the definition of the nonlocal stress tensor, and the strain and the damage were taken as nonlocal variables in the free energy function. Based on this, the damage constitutive equations and the evolution laws were derived via fundamental constitutive axioms and thermodynamic laws. The stress-strain relationship was formulated as a convolutional integral form to describe the influence of long-range forces among the microstructures of the material, which was consistent with the nonlocal stress definition. Additionally, the nonlocal damage variable was in conjunction with a nonlocal description of the energy release rate [see Eqs (23), (28) and (29)], so it could be ensured that the nonlocal damage formulations complied with the thermodynamic principles. In addition, since the surface-residual-based nonlocal stress was utilized, the nonlocal surface residual (surface-induced traction) was incorporated into the boundary condition, which reflected the long-range actions induced by the boundary surface. 2. By applying the thermodynamically based nonlocal damage model to a one-dimensional example, the strain and damage localization profiles were numerically simulated for elastic-brittle materials. Through a nonlocal kernel, the intrinsic length scale was introduced into the present model in a natural manner to act as a "localization limiter", which resulted in the strain and damage in the bar being localized into a small central region whose size was comparable to the internal characteristic length. Therefore, there was no need to postulate a localization zone with fixed size beforehand [12,29,30], which made the analysis both more simple and more physical. However, two major differences appeared in the strain distribution as compared to former research: (i) At both sides of the localization zone, two elastic unloading areas existed where the strain decreased sharply. This strange phenomenon might have arisen from the oscillation induced by numerical integration. (ii) The strain at the end points rose slightly due to the effects of the nonlocal surface residual, but this boundary effect was insignificant because the surface residual was much smaller than the external applied load. 3. Corresponding to the strain and damage distribution profiles, the stress-strain responses in and out of the localization zone were computed. It was clearly seen that in the nonlocalizing part, the stress-strain curves retained linear elastic patterns throughout the loading process. Within the localization zone, as the damage grew drastically, the tangential modulus of the stress-strain curves turned negative and the stress dropped from a peak value to almost zero with the increasing strain, exhibiting obvious softening features. Moreover, by comparison, we found that the stress-strain curve at the boundary point deviated slightly from those in the interior of the bar, which was in accord with the insignificant boundary effects observed in the strain localization behavior. 4. The influence of varying intrinsic length scales on the size of the localization zone was studied by simulating the strain localization behavior for different internal characteristic lengths. The results suggested that the characteristic length was proportional to the localizing width. As the internal length was reduced, the strain was more concentrated in a smaller localization zone to form a higher peak value. This correlation between the internal length scale and the size of the localization zone provided us a reasonable way to explain and simulate the evolution law of the localization zone in the loading process, which could be realized by introducing varying internal characteristic lengths into the nonlocal damage model.