Influence of nonlinear spatial distribution of stress and strain on solving problems of solid mechanics

The stress and the strain should be defined as statistical variables averaged over the representative volume elements for any real continuum system. It is shown that their nonlinear spatial distributions undermine the classical framework of solid mechanics and may cause non-ignorable errors to the solutions. With considering the high-order gradients of the stress and the strain, a two-step solution scheme is proposed to compensate for the influence. Through a revisit to three simple but typical problems, i.e., the hole size-dependence of the fracture strength of perforated plates, the indentation depth-dependence of the measured elastic modulus by micro-indentation tests, and the tensile necking of metallic materials as well as hyperelastic materials, the effect of the nonlinear spatial distribution of stress and strain on solving these problems is illustrated. The observed size effect and the instability of deformation can be quantitatively explained if the effect is properly considered by the proposed method.

1 Fundamental assumptions of the classical framework of solid mechanics As one of the most successful theories for modeling the deformation of materials and structures, the mechanics of solids has witnessed a glorious development since Augustin-Louis Cauchy introduced the concepts of stress and strain in 1882. However, the great success covers up the fact that some seemingly quite simple problems, such as the size-dependence of the fracture strength of perforated plates [1] and the tensile necking of hyper-elastic materials [2][3] , yet cannot be predicted by the well-developed numerical tools built on the classical theory of solid mechanics. A revisit of the fundamental assumptions of solid mechanics may shed some light on identifying the causes to the above problems and further on ways to improve the classical framework of solid mechanics.
As the most fundamental concept in solid mechanics, the stress and the strain are defined on idealized, infinitely small material points. However, any real solid continuum is inherently heterogeneous. To exploit the classical theory of solid mechanics, the real heterogeneous system has to be modeled as a continuous mass, i.e., a continuum, after a homogenization performed over the representative volume elements (RVEs) defined at a lower scale. Although not explicitly stated, it is assumed that the statistical averaging of the stress over each RVE is equivalent to the stress at the idealized material point. And the stress at the centroid of the RVE is assumed to be equal to the stress averaged over the RVE (see Fig. 1(a)), i.e., where Ω RVE denotes the domain of the RVE, and x c denotes its centroid. · Ω denotes the average of a variable over Ω. Similarly, the assumptions also hold for other field variables.
Another fundamental assumption is that the behavior of a material is inherent to itself and independent of the geometric structure and the loading conditions (see Fig. 1(b)). Thus, the mechanical behavior of complex structures under different loads (which can be decomposed into combinations of simple loads) can be predicted by studying the mechanical behavior of simple material samples under simple loads. It is the constitutive relationship that bridges the mechanical behaviors of the structure and the material. This assumption is an indispensible ingredient of solid mechanics. However, this is only an empirical assumption rooted in scaling experiments. It is logical to accept that the averaged response of a RVE can be inferred from the averaged response of the material sample when their loading conditions are similar. But this does not guarantee that the microscopic deformation patterns of a material remain the same under different loads. In the present study, we find that the a common phenomenon, i.e., the nonlinear spatial distribution of stress and strain, undermines the validity of the above assumptions and may cause non-ignorable errors to the solutions based on the classical framework of solid mechanics. The effect will be discussed in Section 2 in terms of simple one-dimensional analysis, and then possible remedies will be proposed. In Section 3, three very typical yet incompletely solved problems, i.e., the hole size-dependence of the fracture strength of perforated plates, the indentation depth-dependence of measured elastic modulus by micro-indentation tests, and the tensile necking of elasto-plastic materials as well as hyperelastic materials, will be reanalyzed in terms of the newly proposed measures. Further discussions will be given in Section 4, and a conclusion will be drawn in Section 5.
2 Effects of nonlinear spatial distribution of stress and strain on fundamental assumptions As mentioned above, the stress and the strain are defined as statistical variables averaged over RVEs for any real continuum system. And the statistical quantities are defined at the centroids of the RVEs. For simplicity, we choose a scalar field variable, i.e., p(x), to illustrate the effect of the nonlinear spatial distribution on the statistical quantity, i.e., p(x) ΩRVE . For conciseness of notation, p(x) ≡ p(x) ΩRVE , and p(x) ≡ p(x c ) hereinafter. By the Taylor series expansion, it is trivial to show that if the field variable follows a constant or a linear spatial distribution, the averaged quantity is equal to the counterpart defined at the centroid of the RVE, i.e., p(x) = p(x), and therefore, the field variable defined within the modeled continuum can reliably describe the field within the real heterogeneous system statistically. However, this equivalence does not hold if the spatial distribution is nonlinear. It is interesting to note that the representative field defined at the centroid of the RVE can be recovered through the Taylor series expansion of the averaged counterpart by supplementing the information from the neighboring RVEs, i.e., by incorporating the gradients of the field, where h denotes the size of the RVE. ∆ and ∆ 2 denote the Laplacian operator and the biharmonic operator, respectively. For a real heterogeneous system, by considering that the size of the RVE is finite, the high-order derivative terms in Eq. (2) cannot be omitted if the spatial variation of the field variable is noticeable. For the highly nonlinear distribution (for example, the stress close to a crack tip or at the grain boundaries) or even discontinuous distribution (for example, the stress and the strain across micro-cracks or shear bands), incorporating more terms of high-order gradients can only alleviate but cannot completely remove the discrepancy between the averaged quantity and the representative counterpart defined at the centroid of the RVE. A correction term in an integral form is more favorable for the highly nonlinear cases, i.e., However, by considering that the spatial distribution is an unknown a priori, additional efforts are needed to seek the formulation of the correction term. It may be a promising way to seek the specific form of the correction term within the framework of thermodynamics as proposed by Wang [4] .
The discrepancy discussed above also exists in other high-order tensor fields such as the stress and the strain. If their spatial distribution is only modestly nonlinear, the representative stress defined at the materials points (i.e., the centroids of the RVEs) can be recovered from the microscopically averaged counterparts by incorporating the high-order derivatives, For isotropic materials, the simple tensor algebra shows that the components of the Laplacian term and the biharmonic terms can be obtained by performing the Laplacian and the biharmonic operations only on the individual component of the stress tensor. Likewise, if the stress follows a highly nonlinear distribution within the RVE or the whole domain, it is more favorable to define the correction term in an integral form.
The above analysis clearly indicates that the nonlinear spatial distribution of the stress and/or the strain, which is rather common in practical problems, undermines the classical framework of solid mechanics. The nonlinear distribution affects not only the solution to the stress and the strain but also the constitutive modeling of materials by considering that the constitutive model is usually defined in terms of stress and strain. If such an effect is not properly considered, the measured material behavior may not be the 'true' behavior.
To further illustrate the effect and to investigate an appropriate solution technique to the problems involving nonlinear distribution of stress and strain, three simple but typical problems will be discussed in the following section. A two-step 'predict-correct' scheme is proposed to remedy the classical solution method. We limit our discussion to the modestly nonlinear cases, i.e., the correction terms can be expressed in differential forms, and the investigation of the highly nonlinear cases is left for our immediately following work.

Two-step solution scheme
For conciseness of notation, we collect all the correction terms and denote the lumped term with the superscript of 'C'. For the Cauchy stress, where σ(x) is the conventional stress, and the corrected stress σ C (X) is defined either in differential form or in integral form. Similarly, the displacement field can be corrected as follows: If a proper measure is chosen for the strain (i.e., the engineering strain in the updated Lagrangian formulation), a similar correction in linear form can be formulated, The distinction between the original configuration and the deformed configuration is left for the discussion of specific problems. We assume that the true constitutive behavior of materials still can be described in terms of stress and strain, which provided that the effect of the nonlinear distribution has been properly considered during processing the measurements of stress and strain. The corrected governing equations can thus be formulated as where b is the body force. S sdv denotes the state dependent variable. ∂ d Ω and ∂ n Ω, respectively, denote the boundaries where the displacement u and the traction t are specified. And n denotes the outer normal direction of the boundary. The present study focuses on modestly nonlinear cases, and the correction terms are defined by the high-order derivatives of the tensor variables. The high-order differential operator(s) much complicates the solution to the above governing equations. Considering that the solution to the classical governing equations, i.e., the equations without the correction terms, is a good trial solution, a 'predict-correct' scheme is proposed to solve the corrected governing equations. In the 'predict' step, the classical governing equations are solved either analytically or numerically. The second order derivatives of the stress and the strain are then evaluated. Given the length scale of the RVE, the magnitudes of the corrected stress as well as the corrected strain can be readily estimated. If the estimated corrected terms are non-ignorable, a 'correct' step should be followed to compute the corrected terms. For linear elastic problems, an auxiliary problem can be defined to solve the corrected stress and the corrected strain. The governing equations of the auxiliary problem are listed as follows: where C is the stiffness matrix of the elastic material. u C and σ C are computed from the solution obtained by the 'predict' step. For nonlinear problems, the nominal constitutive relationship between the corrected stress and the corrected strain usually cannot be formulated separately as in the elastic problem. The full problem and the auxiliary problem have to be solved iteratively in a coupled way.
3.2 Hole size-dependence of fracture strength of perforated plate under uniaxial loading The stress distribution within a perforated plate under the uniaxial loading (see Fig. 2) is a classical problem well discussed in the textbooks of solid mechanics. If the dimension of the plate is sufficiently larger than the size of the hole, the stress on the periphery of the hole is independent of the absolute size of the hole and therefore, the critical far-field stress upon the initiation of the fracture should be insensitive to the hole size. However, it was recognized that the critical stress decreases with increasing the hole size and approaches to a horizontal asymptote for large holes [1] . To explain the hole-size dependence of the fracture strength, the theory of stress gradient is proposed, and the variation of stress gradient across a length scale is included to remedy the classical solution [1] .
In the present study, the size effect of the sidewall fracture of the cavities (perpendicular to the loading direction) within rocks is analyzed [1] (see Fig. 2). The radius of the holes ranged from 1.6 mm to 31.0 mm. Both the width and the height of the samples are 10 times larger than the hole radius. Therefore, the boundary effect is negligible, and the sample can be modeled as an infinite solid. According to the classical solution, the stress intensity factor at the locations of the strain gauges is a constant, and the classical solution gives the approximate normal stress in the load direction, i.e., σ yy , along the x-axis as [5] where σ 0 is the axial load, and σ 0 = σ yy /3 at the idealized location of the strain gauges, i.e., r = a. By considering the high nonlinearity of the spatial distribution of σ yy , this solution is only an approximation to the true solution, because the corrected term of stress is neglected during solving the problem.
For the present problem, if the first correction term, i.e., the Laplacian term, is included in the stress and the strain, the analytical solution to σ yy along the x-axis can be sought analytically and is given as The corresponding axial load upon fracture initiation can be derived by substituting r = a into the above equation, The simple formula can describe the dependence of the axial load on the hole size very well except the smallest one (see Fig. 3). The best fitting parameters σ f 0 = 21.0 MPa and h = 5.07 mm.
If both the Laplacian and the biharmonic terms are included, the analytical solution to the corrected problem can be sought by following the same procedure: The corrected intensity factor of the axial load upon fracture initiation can be easily found by substituting r = a into the above equation (see Fig. 3). The best fitting parameters σ f 0 = 21.0 MPa and h = 5.02 mm, The two solutions yield almost the same results, when the radius of the hole is larger than 6 mm, and the latter solution can describe the observed size dependence slightly better.
Given the solution with consideration of the hole size, the effect of the stress distribution along the x-axis can be studied (see Fig. 4). Not surprisingly, when the radius of the hole is much larger (e.g., 10 times larger) than the length scale, the stress distribution can be reliably described by the classical solution. However, if the radius of the hole is comparable to the length scale, the stress intensity factor decreases from the well-known value of three if the length scale increases (in other words, the normalized radius of the hole decreases). This is consistent with the above experimental results. It is very interesting to note that if the radius of the hole is smaller than the length scale, there is a stress drop behind the edge of the hole. Then the stress increases again after a small distance and this possibly contributes to form the 'process zone' ahead of a crack tip [6][7] .  Fig. 4 Influence of hole size on stress distribution (color online)

Depth-dependence of reduced modulus extracted from micro-indentation tests
Nano-and micro-indentation tests are possibly the most commonly applied methods to characterize the mechanical behavior of small-sized materials and thin films. Unlike the conventional uniaxial tensile test from which the mechanical parameters can be directly determined by the recorded stress-strain curve, proper indentation models have to be exploited to infer the mechanical parameters through a reverse analysis. The measured, or more exactly, the inferred mechanical parameters of various materials including metals and polymers, have long been found to be dependent on the indentation depth. The indentation size effects in metals are widely discussed in the literature and are commonly attributed to geometrically necessary dislocations [8] . However, the size effect characteristics in polymers are more complicated than metals as the size effect in polymers has been also observed in elastic deformation [9] .
As discussed above, the variation of the gradient of a physical variable within a length scale is equivalent to the Laplacian term in the corrected formulation of the variable (see Eqs. (2) and (5)). By considering that the stress and the strain within the sample are nonlinearly distributed beneath the indenter, the complex spatial distribution of the contact pressure is supposed to influence the interpretation of the indentation data. In order to separate the possible effect from other effects, we performed micro-indentation tests on the nitrile-butadiene rubber (NBR) to ensure the deformation is elastically dominant. In order to illustrate that the effect of even a modestly nonlinear distribution of stress can be noticeable, a spherical-conical indenter was chosen of which the radius of the spherical tip is 100 µm, and the included angle of the cone is 90 • (see Fig. 5). The classical Oliver-Pharr's method is used to infer the reduced modulus from the unloading curve.
The reduced modulus inferred by the indentation model is shown in Fig. 6 and shows a noticeable dependence on the contact radius, or equivalently, on the indentation depth. For the instrumented indentation tests, the mechanical parameters such as the reduced modulus and the hardness are inferred from the measured indentation force (or the total contact force) and the contact area in terms of the theory of contact mechanics [10] . For parabolic punches, Hertz [11] assumed that the contact pressure between the indenter and the elastic half-space follows the distribution like where r is the radius of the measure pressure, and a is the contact radius; p 0 is the maximum pressure of contact between the indenter and the material sample. The relation with the geometry results in the radius of the area of contact as where R is the radius of the indenter, and E * = E 1−v 2 is the reduced modulus. The total indentation load can be readily obtained by integrating the pressure over the contact region, By combining Eq. (16) and Eq. (17), the reduced modulus can be related to the deduced contact radius and the measured indentation force through, As shown in Eq. (18), only the contact force and the contact area are relevant to the determination of the reduced modulus. Considering the contact force can be accurately measured by the instrumented indentation, the most critical parameter influencing the data analysis is the contact area. However, the area is not measured but deduced from the recorded indentation data in terms of the models of contact mechanics. As discussed in Section 2, the assumed contact pressure exerted on the material points (i.e., Eq. (15)) is supposed to be different from the externally applied pressure because of its nonlinear spatial distribution. If the correction retains only the Laplacian term, the externally applied pressure p 1 (r) can be recovered from the Hertzian distribution by solving the differential equation in the cylindrical coordinate system, The discrepancy between the distributions of the truly applied pressure and the Hertz contact pressure becomes noticeable when the length scale is comparable with or greater than the contact radius. The corrected contact radius is larger than the radius predicted by the Hertz contact theory since the corrected maximum pressure is lower than p 0 . As a consequence, the reduced modulus is overestimated because the actual contact radius is underestimated. The predicted reduced modulus is well consistent with the experimental measurement as shown in Fig. 6. By fitting the experimental measurement, the length scale of the rubber is estimated to be 80 µm, and the true reduced modulus is 10.5 MPa.

Necking of metallic materials and hyperelastic materials under tensile loading
The above two cases illustrate the effects of the nonlinearly distributed stress and strain on solving the linear elastic problems. The effect of the solution to a typical nonlinear problem, i.e., necking under uniaxial tensile deformation, is to be discussed in this section. Often observed in metallic materials and ductile polymers, necking is usually attributed to the strain-softening due to the damaged plasticity of metals or due to the unfolding of polymeric molecules [12] . Interestingly, necking is also observed in some strain-hardening hyperelastic materials under the tensile loading [3] .
In the present study, we adopt a one-dimensional model to study the deformation localization of a rod under uniaxial loading. In order to take advantage of a linear form, the conjugate pair, i.e., the first Piola-Kirchhoff stress P 1 (1st PK stress hereinafter) and the principal stretch λ 1 , are adopted to formulate the governing equation and the constitutive equation. The corrected 1st PK stress and the corrected principal stretch can be defined by simply adding the correction term, Suppose the true constitutive relation has been defined in terms of the 1st PK stress and the principal stretch, i.e., Through the Taylor series expansion, Eq. (21) can be approximated by . After incorporating Eq. (22) into the balance equation, i.e., P 1 = 0, one obtains Before the deformation gets localized, the stretch follows a homogeneous distribution and λ C 1 = λ C 1 = 0. Therefore, P 1 = P 1 =ḟ (λ 1 )λ 1 = 0. Equation (23) implies that only when f (λ) = 0, there is a inhomogeneous solution to the corrected stretch. In other words,f (λ 1 ) = 0 is a necessary condition for the deformation localization of a rod under uniaxial loading. This condition is satisfied around the inflection point of the strain-softening or the strain-hardening phase of materials (see Fig. 7). For hyperelastic materials, previous studies [14] show that the frequently adopted constitutive models such as Neo-Hookean, Mooney-Rivlin, Ogden and Gent, do not allow necking under entirely mechanical loading. This is consistent with the above finding in that these models fail to describe the mechanical behavior at large deformation where the inflection point appears. If only the Laplacian term is included in the correction, the distribution of the corrected stretch can be sought by solving Eq. (23). This is a third-order differential equation and as analyzed above, the equation has double solutions. One represents a homogeneous deformation and another represents an inhomogeneous deformation. Only the latter is of interest for the present study. For the metallic material, the stress-strain curve of a high strength steel (i.e., S700(A) [13] ) is taken to define the constitutive relation expressed by the 1st PK stress and the principal stretch. For the hyperelastic material, the Treloar uniaxial test data of a type of vulcanized rubber [14] is used to define the constitutive relation. Although the measurements may not represent the true behavior of the materials after the necking, they are good candidates to define approximate constitutive relations for semi-quantitative studies.
The predicted distributions of the engineering strain of the high strength steel and the principal stretch of the rubber are shown in Fig. 8. Due to the highly nonlinear distribution, the corrected engineering strain is less severely localized than that of the nominal engineering strain. The observation also applies to the hyperelastic material. The present preliminary semi-quantitative investigation hints that special care should be taken to determine the 'true' post-necking behavior of materials.

Comparison with nonlocal theories and gradient theories
The augmented field variables in the present study share similar mathematical formulations as that of the variables augmented by the nonlocal theories [15] . However, the basic ideas behind the two approaches are different. The major purpose of the nonlocal theories is to describe the nonlocal interaction of material points. The present study emphasizes that the nonlocal or even the global interactions originate from the nonlinear spatial distribution of the field variables. Especially, the nonlinear spatial distributions of the stress and the strain have a fundamental effect not only on the constitutive modeling of materials but also on the formulation and solution to the governing equations.
The basic ideas of the gradient theories [16][17][18] and the present study share similarities in that both methods agree that the higher order gradients of the stress, and the strain are important to describe the mechanical behavior of materials under complex loads. The classical theories define the gradient terms mainly from mathematics while the present study defines the gradient terms mainly from physics, i.e., the homogenization of materials and physical fields. After the homogenization operation, the first-order gradients disappear, and the final mathematical formulations are much simplified in the present study. In the classical theories, the length scale parameters are defined mainly to achieve a consistency of scaling analysis. While in the present work, only one length scale parameter is required and more importantly, which has a clear physical meaning.
Although the earlier developed Cosserat/micropolar theories can be classified as a specific gradient theory [19][20] , the difference between this classical theory and the present work is noticeable. The Cosserat theory is based on the assumption that micro-moments exist at each point of a continuum. The independent kinematic degrees of freedom are the displacements and the micro-rotations. In fact, the micro-rotation can be interpreted as the antisymmetric part of the first-order strain gradient [16] . Accordingly, in addition to the classical Cauchy stresses, independent couple stresses are introduced. The Cosserat theory considers only part of the gradient terms and may explain the size effect observed in the torsion of circular cylinders [21] and in the bending of plates and beams [22] . However, no size effect in tension is predicted by the theory.

Effect of solution to solid mechanics problems
Solution to the three simple but typical problems manifests two major benefits of incorporating the correction terms. The first is the improvement to the solution to the stress (and/or strain) with a drastic spatial variation. Negligence of the correction term may make it fail to explain the observed size effect. The second is the capability to predict the multiple deformation patterns often encountered in the problems involving both material and geometric nonlinearities. Mathematically, each deformation pattern corresponds to one solution to the governing equation. Therefore, a stability analysis should be inserted into the 'predict-correct' scheme in order to identify all the possible deformation patterns, and this will be investigated in our following studies.
The effect of the nonlinear spatial distribution depends not only on the degree of nonlinearity but also on another very important parameter, i.e., the length scale of the RVE. This length scale is defined with respect to averaging the physical fields. It is closely related to but may not be the same as the length scale of the microstructure of materials. As discussed above, if the nonlinearly distributed field shows a drastic spatial variation or spans the whole domain, the length scale is much different from the length scale defined with respect to the microstructure of materials. It is manifested that the length scale should be determined as a separate constitutive parameter besides the conventional mechanical parameters. A proper routine to characterize the parameter is to be investigated.

Conclusions and suggestions for future research
The present study emphasizes that a commonly observed phenomenon, i.e., the nonlinear spatial distributions of stress and strain, affects modeling the mechanical behavior of solid materials and structures. A revisit of three simple but typical problems manifests that a proper consideration of the influence is necessary for the problems of which the stress and/or the strain show drastic spatial variations. In addition to that the stress and the deformation can be more reliably predicted, the deformation patterns for nonlinear problems can be characterized by the augmented governing equations as well.
The proposed remedies seem to be promising but much work is yet to be conducted.
(i) For the highly nonlinear fields, the correction terms are preferentially expressed in integral forms, and specific formulations are to be sought before conducting analysis of practical problems.
(ii) Effective methods should be designed to determine the length scale, which plays an important role in defining the correction terms.
(iii) To obtain the true mechanical behavior of materials during the post-necking phase, proper methods should be investigated to eliminate the effect of the deformation localization.
(iv) Numerical solution technique should be developed for the solution to nontrivial problems. Compared with the classical numerical procedure, further efforts are needed to properly define the boundary conditions and to perform stability analysis in order to identify the possible deformation patterns.
Di YANG, and Yuheng CAO for their support on part of the verification.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.