Silvère Konlambigue
email: silvere.konlambigue@datahertz.fr

Jean-Baptiste Pothin
email: jean-baptiste.pothin@datahertz.fr

Paul Honeine

Abdelaziz Bensrhair

† S Konlambigue

Performance Evaluation of State-of-the-Art Filtering Criteria Applied to SIFT Features

Keywords: SIFT, Feature filtering, Pruning, Entropy, Inner primary ratio (IPR), Saliency, Convolutional Neural Networks (CNN)

Feature matching is an important and crucial task in computer vision. Unfortunately, the features extracted from images are usually redundant or irrelevant, leading to many ambiguities and false positives during matching. This work provides an exhaustive study of state-of-the-art feature filtering strategies investigating different criteria, such as the contrast level, inner primary ratio, entropy, saliency and median value, in the purpose of discarding irrelevant features. The evaluation on the well-known Oxford-5k dataset demonstrates the relevance of these criteria to reduce the amount of data involved in the matching step and the number of false positives, thus leading to faster and more accurate matching with almost no loss in performance. Our results are also compared to convolutional neural networks (CNNs), providing extensive experiments that corroborate recently published work.

I. Introduction

Since the work of Krizhevsky et al. in [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] with AlexNet for visual recognition in ILSVRC (ImageNet Large-Scale Visual Recognition Challenge), the current trend in computer vision has been focusing on deep learning, especially convolutional neural networks (CNNs). For these CNNbased methods to perform well, a large amount of data and time are needed to train the network (e.g. 1.2 million training images, 50,000 validation images, and 150,000 testing images with roughly 1000 images in each of 1000 categories for ILSVRC). One solution to these drawbacks is provided by the so-called transfer learning method, where one uses a pre-trained (generic) network and finetunes its weights in order to address a specific task at hand [START_REF] Liu | Scene recognition for indoor localization using a multi-sensor fusion approach[END_REF]. This procedure results in a black box framework for which one does not control the various treatments carried out or its behavior. Also, even if it requires less paul.honeine@univ-rouen.fr † abdelaziz.bensrhair@insa-rouen.fr data than when starting from scratch, a large amount of data is still needed in practice. For all these reasons, researchers have been recently re-investigating methods like Scale Invariant Feature Transform (SIFT), since it overcomes the aforementioned drawbacks. Moreover, recent studies corroborate connections between SIFT and CNN [START_REF] Zheng | SIFT meets CNN: A decade survey of instance retrieval[END_REF], with comparable performances [START_REF] Sachdeva | Performance evaluation of SIFT and Convolutional Neural Network for image retrieval[END_REF]. The SIFT method [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] is still considered a reference in computer vision algorithms due to its robustness to geometric and photometric transformations, as well as the high discriminative power of its 128-element vector descriptors. All those properties make the SIFT method to be preferred to other image matching methods [START_REF] Chatoux | Comparative study of descriptors with dense key points[END_REF]. For computer vision applications such as object recognition [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF], image retrieval [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] and many others, feature matching is a very important task. In the naive approach, namely the exhaustive search, finding matches between two images consists in comparing each feature of the first image to all features from the second looking for their nearest neighbor in the descriptor space.

Local descriptors methods such as the SIFT method, often leads to large features databases which, combined with its descriptors dimensions, poses two problems: storage and time consuming of the nearest neighbor search. One idea to speed up the matching process is to reduce descriptors size. This can be done in several ways such as PCA-SIFT [START_REF] Ke | PCA-SIFT: A more distinctive representation for local image descriptors[END_REF], Reduced-SIFT [START_REF] Ledwich | Reduced SIFT features for image retrieval and indoor localisation[END_REF] or SURF [START_REF] Bay | SURF: Speeded Up Robust Features[END_REF]; see for instance [START_REF] Wu | A comparative study of SIFT and its variants[END_REF] for a comparative study. Another way to speed up matching is to look for approximated nearest neighbors (ANN) instead of the real ones. The binarization of descriptors [START_REF] Peker | Binary SIFT: Fast image retrieval using binary quantized sift features[END_REF], [START_REF] Zhou | BSIFT: toward data-independent codebook for large scale image search[END_REF] or binary features [START_REF] Heinly | Comparative Evaluation of Binary Features[END_REF] allow also a phenomenal saving of time during the matching through the use of the Hamming distance, instead of the Euclidean distance usually used to compare real or integer value descriptors; the price to pay is usually a reduced performance compared to the latter. These descriptors also offer a huge gain in storage precisely because binary words take up less space than integer or real ones.

Another major approach, leading to reduction of storage space and matching speed-up is feature filtering. In practice, many of the found matches are false positives [START_REF] Hartmann | Predicting matchability[END_REF] and can skew the results when their number is too high, even with the use of a robust geometric fitting method such as RANSAC [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF]. Feature filtering is a very simple way of reducing the number of false positives by eliminating irrelevant features according to a certain criterion. Most prominent methods are based on entropy [START_REF] Zivkovic | On matching interest regions using local descriptors -can an information theoretic approach help?[END_REF], median [START_REF] Zhou | BSIFT: toward data-independent codebook for large scale image search[END_REF] value, contrast level [START_REF] Foo | Pruning SIFT for scalable nearduplicate image matching[END_REF] and several others. As mentioned, beyond reducing the number of false positives, feature filtering has a double advantage: reducing the number of features involved in the matching process and thus accelerating it and gaining in storage space due to this reduced number of features. These filtering criteria are relatively simple to implement, very inexpensive and can be used upstream to the aforementioned methods (reduced-size descriptors, binary descriptors and ANN) to further speedup the matching task.

In this paper, we investigate various feature filtering criteria, implement them and evaluate their performance, in order to allow the readers to make a better choice according to their needs and the compromise they are ready to accept. Experiments are conducted on the Oxford-5k dataset [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF]. This dataset is quite challenging due to the amount of features and the substantial variations in scale, viewpoint, and lighting conditions. The feature matching strategy used in the experiments is the nearest neighbor distance ratio (NNDR) [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] and performance is evaluated using the mean average precision (mAP) measure [START_REF] Nikolopoulos | High order pLSA for indexing tagged images[END_REF]. The conducted study provides an in-depth comparative analysis with CNN results and other state-of-the-art methods [START_REF] Sachdeva | Performance evaluation of SIFT and Convolutional Neural Network for image retrieval[END_REF], [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF], [START_REF] Jégou | Hamming Embedding and Weak Geometry Consistency for Large Scale Image Search -extended version[END_REF]. This paper is organized as follows: in Section II, we give a brief review of the SIFT method and the matching process. Filtering criteria and their operations are presented in Section III. Section IV presents the dataset, the experimental settings, and the performance evaluations. Section V concludes this paper.

II. SIFT-based Image Matching

The SIFT method is one of the most investigated computer vision algorithms [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF]. It operates in two steps. In the first step, called keypoints detection, a given image is processed to extract distinctive points using the Difference-of-Gaussian (DoG) detector. Looking for these points at different scales, by constructing a Gaussian scalespace, allows the SIFT method to be robust to scale changes. Note that the Gaussian scale-space computation is very time-consuming, which can be overcome using for instance extended box filters instead of commonly used Gaussian filters [START_REF] Konlambigue | Fast and accurate gaussian pyramid construction by extended box filtering[END_REF]. A keypoint is represented as a tuple (x, y, σ, θ, R) where x and y are its coordinates in the image, σ the scale, θ the orientation and R the normalized image patch around the keypoint.

In the second step, called keypoints description, 16 small 8-bin orientation histograms are computed from 4 × 4 sub-regions covering the patch R, and then concatenated to form the 128-element keypoint descriptor: d={d 0 , . . . , d 127 }. Histograms are weighted by the gradient magnitude and a Gaussian window to give less emphasis to points that are far from the keypoint position. The gradient magnitude m(x, y) and orientation φ(x, y) are computed as follows:

m(x, y) = R x (x, y) 2 + R y (x, y) 2 , ( 1
)
φ(x, y) = arctan R y (x, y) R x (x, y) , (2)
with

R x (x, y) = R(x + 1, y) -R(x -1, y), (3)
R y (x, y) = R(x, y + 1) -R(x, y -1). ( 4
)
As R is already normalized by σ and θ, the resulting descriptor is thus scale and rotation invariant.

To find matches between two images, especially in the exhaustive search approach, one has to compare each feature from the first image to all the features from the second, looking for their nearest neighbors in descriptor space. The nearest neighbors are typically evaluated using metrics such as the Euclidean distance. In the NNDR strategy, a match between two features is accepted if and only if the ratio between the distance of the first and the second nearest neighbors is lower than some threshold, namely

d nn1 d nn2 < τ, ( 5
)
where 0 < τ < 1 denotes a predefined threshold.

III. feature filtering Criteria In this section, some criteria and their technical details are presented in order to eliminate irrelevant SIFT features. See Fig. 1 for an end-to-end image retrieval system overview.

A. Contrast-based method (CP)

In the SIFT method [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], weak contrast points, points with low DoG values, are filtered-out with respect to a threshold value. This filtering strategy improves the robustness of the detected features. In [START_REF] Foo | Pruning SIFT for scalable nearduplicate image matching[END_REF], Foo et al. took advantage of this contrast-based filtering strategy and introduced a new pruning scheme. Instead of setting up a threshold value, they suggested to select the top N most significant keypoints ranked by their contrast values. Note that this pruning scheme only affects images that have more than N features.

B. IPR-based method

Treen et al. introduced in [START_REF] Treen | Efficient SIFT matching from keypoint descriptor properties[END_REF] a new filtering criterion based on the inner primary ratio (IPR) value formulated for a 128-element descriptor as:

IPR = d 2 40 + d 2 48 + d 2 72 + d 2 80 127 i=0 d 2 i . ( 6
)
High IPR values indicate that the four inner elements (i.e., in the numerator) capitalize most of the descriptor energy, while low IPR values indicate that the energy is distributed on each bin of the descriptor. Only keeping low IPR descriptors reduces ambiguities and false matches from the matching. For example in [START_REF] Treen | Efficient SIFT matching from keypoint descriptor properties[END_REF], the IPR threshold value is set to IPR th = 0.235 after extensive experiments on a synthetic benchmark.

C. Entropy-based method

The entropy is a metric from information theory that quantifies the amount of information contained in a message. In [START_REF] Zivkovic | On matching interest regions using local descriptors -can an information theoretic approach help?[END_REF], Zivkovic et al. evaluated the amount of information hold in a 128-element vector descriptor using the definition:

H 1 (d) = - 127 i=0 p i (d) log 2 p i (d), (7)
with

p i (d) = d i 127 k=0 d k . ( 8
)
On the other hand, Dong et al. considered in [START_REF] Dong | High-confidence near-duplicate image detection[END_REF] another approach for estimating the descriptor entropy. In their definition, as the descriptor is integer value d k ∈ {0, . . . , 255}, the entropy is evaluated in the form:

H 2 (d) = - 255 i=0 pi (d) log 2 pi (d), (9)
with pi (d) = |{k | d k = i}| 128 , ( 10
)
where | • | represents the number of elements in the set.

D. Saliency-based method

In a field of view, points that stand out from their surround, namely the salient points, visually captivate attention (focus-of-attention) in the scene. Given an input color image, the method proposed by Itti et al. in [START_REF] Itti | A model of saliency-based visual attention for rapid scene analysis[END_REF] computes a saliency map, namely scalar quantities at each location in the image. This map is computed by combining feature maps from a multi-scale feature extraction method, the center-surround, on intensity, color channels and orientation images derived from the input image; see Fig. 2, obtained using the code provided by Harel et al. in [START_REF] Harel | Graph-based visual saliency[END_REF].

In order to keep only the most relevant features, one has to select features with the highest saliency values. The method was used in [START_REF] Sergieh | Improving SURF image matching using supervised learning[END_REF] on the SURF features in order to filter out irrelevant features by classifying them as significant (i.e., highly salient) or insignificant (i.e., less salient). The saliency threshold value for this binary classifier varies from 0.2 to 0.3. In [START_REF] Zhao | Fast SIFT scene matching algorithm based on saliency detection and frequency segmentation for downward-viewing images[END_REF], [START_REF] Awad | Saliency filtering of SIFT detectors: Application to cbir[END_REF], the saliencybased feature filtering was used on SIFT features. More specifically, instead of defining a global threshold for all images as in [START_REF] Sergieh | Improving SURF image matching using supervised learning[END_REF], the threshold is made depend on each image and set to 3 times the average salience of the image in [START_REF] Zhao | Fast SIFT scene matching algorithm based on saliency detection and frequency segmentation for downward-viewing images[END_REF]. In [START_REF] Awad | Saliency filtering of SIFT detectors: Application to cbir[END_REF], no threshold was set but features were pruned in order to select the most relevant, such that the ratio between the kept features to their initial number meets the user-defined ratio, varying from 10% to 40%.

E. Filtering by descriptors median value

In the SIFT method [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], features points located along the edges are detected based on their principal curvatures and eliminated because they are unstable and prone to false matches. Despite this strategy, Zhou et in [START_REF] Zhou | BSIFT: toward data-independent codebook for large scale image search[END_REF] that some of these points are still present and degrade the matching accuracy. Similar to Treen et al. [START_REF] Treen | Efficient SIFT matching from keypoint descriptor properties[END_REF], Zhou et al. also noted that, for these keypoints, most of the coefficients of their descriptors bins are of low magnitude leading to low median values. One is then able to reduce irrelevant points by filtering out low median value features.

In [START_REF] Zhou | BSIFT: toward data-independent codebook for large scale image search[END_REF], the threshold value is set to M th = 5.5 for integer values descriptors in the range [0 , 512].

IV. Experiments and Results

A. Dataset

To effectively evaluate performance, experiments must be conducted on a large and quite challenging dataset. For this purpose, we evaluate each of the filtering criteria on the Oxford-5k dataset [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] implementing an image retrieval system. This dataset contains 5062 high resolution images distributed in 11 different landmarks, denoted L1 to L11, and indexed by 55 query images, 5 images per landmark with associated region of interest (ROI). Each query has four ground-truth labels:

• Good -the ROI is clearly visible • Ok -more than 25% of the ROI is visible • Junk -less than 25% of the ROI is visible (including high occlusion and high distortion) • Absent -the ROI is not present (namely, images that are not in the first three classes) In performance assessment, according to the protocol of the Oxford-5k dataset, Good and Ok labels are considered true positives, while Absent labels are false positives. Junk labels are considered neutral and do not affect the performances. Considering the SIFT method, these 5062 database images yield about 17 millions features. Note that our experiments implement the Multiple Match Removal (MMR) method [START_REF] Mikolajczyk | Indexing based on scale invariant interest points[END_REF].

B. Performance measure

To evaluate the retrieval performance on the Oxford 5k dataset, the mean average precision (mAP) is used, as described in the following.

In the matching stage, each query is compared to all of the database images. A descending ranked list is then returned according to the level of similarity defined as:

sim(I q , I s ) = m q,s m q , ( 11
)
where m q,s is the number of features that satisfy (5), between query image I q and database image I s , and m q is the number of the query features. Note that sim(•, •) ∈ [0, 1] with sim(•, •) ≈ 1 (resp. ≈ 0) means images are very similar (resp. dissimilar).

From the ranked list of images, one evaluates the average precision (AP) for query q as follow [START_REF] Nikolopoulos | High order pLSA for indexing tagged images[END_REF]:

AP(q) = N r=1 P (r) × rel(r) number of relevant images , ( 12
)
with r the retrieved image rank, N the number of retrieved images, rel(r) an indicator equaling 1 if the image at rank r is relevant relatively to q and 0 otherwise. P (r) is the precision at cut-off r defined as the ratio of the number of relevant retrieved images at rank r, over r.

As the AP is calculated for each query q, the mAP is obtained by averaging over all queries:

mAP = Nq q=1 AP(q) N q × 100%, (13)
where N q is the total number of query images.

To compare the methods we evaluate the change in mAP, relatively to the baseline mAP, against the percent- age of filtered data. The reduction level ρ over storage in the database is defined as:

ρ = 1 - s n s N s × 100%, (14)
where n s , resp. N s , is the number of remaining features, resp. the total number of features before the filtering, from the database image I s .

C. Performance evaluation

The results of the experiments are given in TABLE I. In the experiments, the baseline is the original SIFT. Its mAP value on the Oxford-5k dataset is mAP = 65.91%, which is underlined in TABLE I. In this table, the values in columns L1-L11 represent the mAP values according to each landmark, and the columns "mAP" and ρ are respectively the averages over all landmarks and the reduction rates defined in [START_REF] Zhou | BSIFT: toward data-independent codebook for large scale image search[END_REF]. The first six lines recall the state of the art; see [START_REF] Zheng | SIFT meets CNN: A decade survey of instance retrieval[END_REF] for a complete comparison of SIFT and CNN based methods. Note that we match and even improve these performances without using any codebook.

In the literature, the reduction rates often put forward are around 12%. Fig. 3 shows the loss in mAP according to the percentage of reduction for each of the aforementioned methods. The right figure is a simple zoom out of the left one. From this figure, three observations can be made. First, none of the filtering criteria improves the mAP. Each of them results in a loss of performance depending on the rate of reduction. Secondly, for reduction rates of less than about 15%, most methods seem equivalent to each other, except the contrast-based method (CP) which appears to be the best one. Finally, at higher rates, the entropy and median provide the lowest drop in mAP. This is expected because it is well understood that a weak entropy descriptor must have a low median. The IPR-based method is less good but joins those of the entropy when reduction factor is important. Note that for reduction rates greater than 40%, the saliency-based approach is worst than random filtering.

In TABLE I we give details of mAP for thresholds that do not exceed 12% threshold. Also, in our experiments, an IPR threshold value IPR th = 0.235 allows a data reduction of about 12% with only 1% loss of mAP and therefore corroborates the recommendations made in [START_REF] Treen | Efficient SIFT matching from keypoint descriptor properties[END_REF]. It is worth noting the difference in thresholds for both types of entropy. Although it is the same descriptor entropy that is evaluated in both approaches, the difference in the thresholds values comes from the nature (integer or real) of the descriptor, as described in ( 7) and ( 9) in Section III-C. Through experiments, we have observed that one should pay attention to this detail as it can have a highly negative impact on the performance; see Fig. 4.

V. Conclusion

In this paper we have reviewed a number of criteria proposed in the state of the art that can be used for SIFT feature filtering. The obtained results match and even improve these performances of state-of-the-art methods, including deep neural networks with CNN. Although some criteria seem equivalent in terms of mAP performance, according to a certain rate of reduction, two important aspects need to be taken into account when choosing a criterion: the ease of implementation and the cost in computing time. On these bases, we advise the criteria based on the entropy, for the following reasons: it is directly computed from the descriptor, practical when one uses an already implemented version of a computer vision algorithm, and also because a filtering according to this criterion is based on the richness of the descriptor.

Fig. 1 .

 1 Fig. 1. The SIFT-based image matching schema with the features filtering block as a post-processing step.

Fig. 2 .

 2 Fig. 2. Example of a saliency map: left -original image, rightsaliency map (overlay), obtained using Itti et al. algorithm [26]. The saliency values increase here from blue to red.

Fig. 3 .

 3 Fig. 3. Performance comparison of filtering criteria. The right figure is a zoom in for ρ ≤ 25%.

Fig. 4 .

 4 Fig. 4. Histogram distribution of entropy H 1 and H 2 extracted from the Oxford-5k dataset.

SIFT+BoVW[500k]+td-idf [32] -----------61.30 -SIFT+BoVW[20k]+X [33] ---------

TABLE I

The mAP (%) performance evaluation on the Oxford-5k dataset. The values for every landmark (L1 to L11) are given when available, as well as their average (column mAP ) and the reduction level ρ (%) over the database. The first five lines are taken from the respective references. BoVW


M B Rubin
email: mbrubin@tx.technion.ac.il

S Forest

Analysis of material instability of a smooth elastic-inelastic transition model

Keywords: elastic-inelastic transition, limit load, localization, instabilities, smooth model

Analysis of material instability of a smooth elastic-inelastic transition model

Miles B. Rubin, Samuel

1. Introduction [START_REF] Lubliner | A new model of generalized plasticity and its numerical implementation[END_REF] generalized previous work by the first author and developed a rate-independent overstress model with a smooth elastic-inelastic transition and a finite elastic range for which the inelastic deformation rate depends linearly on the stress rate. Their model included both kinematic and isotropic hardening. This model was generalized for large deformations in [START_REF] Panoskaltsis | On large deformation generalized plasticity[END_REF]] but the function h used in this generalization Preprint submitted to International Journal of Solids and Structures January 25, 2020 limited attention to sharp elastic-inelastic transitions. [START_REF] Einav | The unification of hypo-plastic and elasto-plastic theories[END_REF] generalized previous hypo-plastic and hyper-plastic models which predict rateindependent smooth stress-strain curves with no finite elastic range and for which the inelastic deformation rate depends linearly on the total strain rate. It is also noted that [START_REF] Kolymbas | Bifurcation analysis for sand samples with a non-linear constitutive equation[END_REF] considered a constitutive equation for large deformations of sand which is modeled as a rate-independent material with no finite elastic range and with nonlinear dependence on the rate of deformation tensor. [START_REF] Hollenstein | Modeling a smooth elastic-inelastic transition with a strongly objective numerical integrator needing no iteration[END_REF][START_REF] Hollenstein | Erratum to: Modeling a smooth elastic-inelastic transition with a strongly objective numerical integrator needing no iteration[END_REF] developed a large deformation model which exhibits a smooth elastic-inelastic transition using overstress with a yield surface. This model can be considered as a unifying generalization to large deformations of the models in [START_REF] Perzyna | The constitutive equations for rate sensitive plastic materials[END_REF]; [START_REF] Lubliner | A new model of generalized plasticity and its numerical implementation[END_REF]; [START_REF] Panoskaltsis | On large deformation generalized plasticity[END_REF]; [START_REF] Einav | The unification of hypo-plastic and elasto-plastic theories[END_REF]] since it models both rate-independent and ratedependent materials exhibiting: smooth elastic-inelastic transitions with finite elastic ranges (like for metals) and smooth stress-strain curves with no finite elastic ranges (like for soils). In this paper, this model will be referred to as a smooth model for short.

The smooth model developed in [START_REF] Hollenstein | Modeling a smooth elastic-inelastic transition with a strongly objective numerical integrator needing no iteration[END_REF]] is an Eulerian formulation of constitutive equations based on the work of [START_REF] Eckart | The thermodynamics of irreversible processes. iv. the theory of elasticity and anelasticity[END_REF] who was the first to propose an evolution equation directly for an elastic deformation measure which includes an inelastic deformation rate. In particular, this Eulerian formulation does not need definitions for total deformation or inelastic deformation measures. The magnitude of the inelastic deformation rate in the evolution equation ( 9) is controlled by a scalar function Γ and its direction is determined by a symmetric tensor A p that is defined by an elastic distortional deformation tensor and it independent of the rate of deformation. In the general form of the smooth model, the scalar Γ is specified by Γ = (a 0 + b 0 g ) + (a 1 + b 1 g ) ε ,

where a 0 , b 0 , a 1 , b 1 are non-negative constants, g is a yield function that is independent of rate, the Macaullay brackets g are defined by

g = max(0, g) , ( 2
)
and ε is the equivalent total distortional deformation rate defined in (7).

The main physical features of this large deformation overstress smooth model can be discussed with reference to the function Γ. First, it is observed that there are no loading and unloading conditions and the model predicts smooth response for all values of the constants a 0 , b 0 , a 1 , b 1 . If a 0 = b 0 = 0, the model predicts rate-independent response, otherwise it predicts rate-dependent response. If b 0 = b 1 = 0, the model has no finite elastic range, which is typical of soils, whereas if a 0 = a 1 = 0, the model has a finite elastic range, which is typical of metals. Moreover, for this latter case the model predicts vanishing inelastic deformation rate when g ≤ 0 and it predicts a smooth elastic-inelastic transition as g becomes positive. This smooth response is a more physical representation of the elastic-inelastic transition in metals than the sharp transition predicted by the standard yield surface model. Also, in this paper, attention is limited to the simplest rate-independent case with a finite elastic range by specifying a 0 = b 0 = a 1 = 0 .

(3) [START_REF] Hill | Some basic principles in the mechanics of solids without a natural time[END_REF] analyzed constitutive equations for solids without a natural time for which the non-symmetric Piola-Kirchhoff stress rate Π is related to the rate Ḟ of the deformation gradient by a fourth order tensor L, such that

Π = L • Ḟ . (4)
This equation describes rate-independent material response when L is a homogeneous function of order zero in Ḟ, which can be nonlinear in F. In general, the elastic response in this model is hypoelastic since a strain energy function may not exist. However, suitable restrictions can be placed on L to ensure that the elastic response in this model is hyperelastic. [START_REF] Hill | Some basic principles in the mechanics of solids without a natural time[END_REF] also considered a special case when Π is determined by the derivative of a potential function with respect to Ḟ. In particular, [START_REF] Hill | Some basic principles in the mechanics of solids without a natural time[END_REF] analyzed general aspects of uniqueness and stability in this class of materials. [START_REF] Rice | The localization of deformation[END_REF] used a piecewise linear form of the rate-type constitutive equation (4) to study localization of plastic deformation. The localization mode was determined by the condition that the jump in the rate of the traction vector applied on a material surface be stationary to a jump in the rate of deformation given by

Ḟ = Ḟ0 + g ⊗ N , ( 5
)
where Ḟ0 is the uniform rate of deformation far from the localization, N is the unit normal to the localized material surface in its reference configuration, g is a vector that needs to be determined by the analysis and a⊗b is the tensor product between two vectors a, b. Rice explained that nonlinear dependence of L in (4) on Ḟ0 can be handled by considering an unlimited number of models based on different L.

Since inelastic distortional deformation rate in the smooth model depends nonlinearly on the current rate of deformation tensor D, the rate of traction applied to the material surface also depends nonlinearly on D. Consequently, application of Rice's criterion to the smooth model requires special attention. It will be shown that if vanishing of the current rate of the traction vector is used as an instability criterion for the smooth model that this criterion yields the same conditions for instability as those for instability of a velocity perturbation from a homogeneous static state with Ḟ0 = 0.

With regard to rate-independent constitutive equations for metals with a yield surface and a finite elastic range it is noted that [START_REF] Petryk | Material instability and strain-rate discontinuities in incrementally nonlinear continua[END_REF] included inertia in a dynamic analysis of instability in an incrementally nonlinear continua with a constitutive equation of the type (4) based on a velocity potential. An important conclusion of that work is that "the derived condition for the instability of uniform quasi-static straining does not coincide with that for instability of equilibrium and both differ in general from the familiar condition of strong ellipticity loss formulated for the tangent moduli". [START_REF] Hutchinson | Shear band formation in plane strain[END_REF] analyzed shear band formation in plane strain for kinematic hardening plasticity and for a plasticity theory based on yield surface corner development. Also, [START_REF] Petryk | Post-critical plastic deformation in incrementally nonlinear materials[END_REF] studied post-critical plastic deformation using the concept of a representative nonuniform solution in a homogeneous rate-independent material characterized by an incrementally nonlinear corner theory of plasticity. Specifically, indeterminacy of a post-critical representative solution was removed by eliminating unstable solution paths. In this regard, it was shown in [START_REF] Forest | A rate-independent crystal plasticity model with a smooth elastic-plastic transition and no slip indeterminacy[END_REF]) that the smooth model eliminates indeterminacy in rate independent crystal plasticity because the model does not require loading and unloading conditions.

With regard to rate-independent constitutive equations for geological materials with no finite elastic range it is noted that [START_REF] Kolymbas | Bifurcation analysis for sand samples with a non-linear constitutive equation[END_REF] analyzed bifurcation in a constitutive equation for sand with nonlinear dependence on the rate of deformation tensor. [START_REF] Chambon | Localization criteria for non-linear constitutive equations of geomaterials[END_REF] also conducted an assumed mode bifurcation analysis of a shear band using a nonlinear constitutive equation of the type (4) that was characterized by ν linear relations, as suggested by [START_REF] Rice | The localization of deformation[END_REF].

Additional analysis of questions of loss of uniqueness and stability of boundary value problems can be found in books like [START_REF] Nguyen | Bifurcation and stability in dissipative systems[END_REF][START_REF] Petryk | Material instabilities in elastic and plastic solids[END_REF][START_REF] Nguyen | Stability and Nonlinear Solid Mechanics[END_REF][START_REF] Besson | Non-Linear Mechanics of Materials[END_REF][START_REF] Bigoni | Nonlinear Solid Mechanics: Bifurcation Theory and Material Instability[END_REF]. The explicit results provided in these books heavily rely on the structure of the rate-independent elastoplastic constitutive equations, and especially on the existence of a multibranch tangent operator relating the strain rate to the stress rate.

The main objective of the present work is to consider material instability of the rate-independent smooth model. An outline of this paper is as follows. Section 2 presents the basic equations of the nonlinear Eulerian formulation of inelasticity. Section 3 presents an analysis of the limit load and Section 4 analyzes material instability by perturbing a homogeneous static deformation. Section 5 presents another analysis of material instability and Section 6 develops the small deformation equations. Section 7 discusses analytical results for example problems. Section 8 presents details of the numerical implementation of the small deformation equations and compares the finite element results with the analytical results for the example problems discussed in Section 7. Finally, Section 9 presents a discussion of the main results. Details of the tensorial notation used in this paper can be found in [ [START_REF] Rubin | Cosserat theories: shells, rods and points[END_REF], Ch. 2].

Basic equations of the Eulerian formulation of inelasticity

This section summarizes the main equations of the nonlinear Eulerian large deformation formulation of inelasticity. Both the standard rate-independent model and a smooth elastic-inelastic transition model developed in [START_REF] Hollenstein | Modeling a smooth elastic-inelastic transition with a strongly objective numerical integrator needing no iteration[END_REF] will be considered.

A material point is located by x in the present configuration and the velocity v is given by v

= ẋ , (6)
where a superposed (˙) denotes material differentiation with respect to time t. The velocity gradient L, rate of deformation D, deviatoric part D of D and the effective total distortional strain rate ε are defined by

L = ∂v/∂x , D = 1 2 (L + L T ) , L = 1 3 (D • I) I + L , D = D - 1 3 (D • I) I , ε = 2 3 D • D = ( 2D 3 ε ) • D , ( 7
)
where I is the second order unit tensor and A • B = tr(AB T ) is the inner product between two second order tensors (A, B).

Motivated by the work in [START_REF] Eckart | The thermodynamics of irreversible processes. iv. the theory of elasticity and anelasticity[END_REF][START_REF] Leonov | Nonequilibrium thermodynamics and rheology of viscoelastic polymer media[END_REF][START_REF] Flory | Thermodynamic relations for high elastic materials[END_REF][START_REF] Rubin | Calculation of hyperelastic response of finitely deformed elastic-viscoplastic materials[END_REF], the Eulerian formulation proposes evolution equations directly for elastic deformation measures. Specifically, in the absence of inelastic dilatational rate, the elastic dilatation is equal to the total dilatation J, which is determined by the evolution equation

J = JD • I , (8)
and the elastic distortional deformation is a second order, symmetric, positivedefinite, unimodular tensor B e , which is determined by the evolution equation

Ḃ e = LB e + B e L T - 2 3 (D • I)B e -ΓA p , A p = B e -( 3 B -1 e • I ) I , Γ ≥ 0 . (9)
In this equation, Γ controls the magnitude of inelastic deformation rate and A p controls its direction. When Γ vanishes the response is elastic and B e is the unimodular part of the left Cauchy-Green deformation tensor so the model reproduces general isotropic elastic response. Also, when loading vanishes with L = 0 and Γ is positive, inelastic deformation rate causes B e to approach the identity tensor, which causes deviatoric stress to vanish. Next, introducing the deviatoric distortional elastic strain g e by

g e = 1 2 B e , B e = B e - 1 3 (B e • I) I , (10)
the yield function g can be specified by

g = 1 - κ γ e , γ e = 3 2 g e • g e , (11)
where κ is a hardening variable that is determined by an evolution equation.

9 I] • D -Γ , (14)
Here, the strain energy Σ per unit mass is specified for a compressible neo-Hookean material in the form

ρ 0 Σ = µ[kU (J) + 1 2 (α 1 -3)] , α 1 = B e • I , k = 2(1 + ν) 3(1 -2ν) , (15)
where (kµ, µ) are the constant zero-stress bulk and shear moduli, respectively, ν is Poisson's ratio, U controls the energy of dilatation, which satisfies the restrictions

U (J) ≥ 0 , U (1) = 0 , dU dJ (1) = 0 , d 2 U dJ 2 > 0 , d 2 U dJ 2 (1) = 1 , ( 16
)
α 1 is pure measure of elastic distortion and ρ 0 is the constant stress-free reference mass density. In the purely mechanical theory, the integral of the rate of material dissipation D is defined by the rate of work done on the body minus the rates of kinetic and strain energies, which in local form is given by

D = T • D -ρ Σ ≥ 0 , (17)
and the Cauchy stress T is determined by the hyperelastic form

T = -p I + T , p = p(J) = -µk dU dJ , T = J -1 µB e , (18)
where use has been made of the conservation of mass to determine the current mass density ρ ρ = J -1 ρ 0 .

Using these constitutive equations, the rate of material dissipation requires

D = 1 2 J -1 µΓ[B e • I - 9 B -1 e • I ] ≥ 0 . ( 20
)
which can be shown to be satisfied for all states of the material [START_REF] Rubin | Calculation of hyperelastic response of finitely deformed elastic-viscoplastic materials[END_REF]. The restriction (17) for the purely mechanical theory is identical to the restriction in the thermomechanical theory [START_REF] Rubin | Calculation of hyperelastic response of finitely deformed elastic-viscoplastic materials[END_REF] that the rate of material dissipation ρθξ ≥ 0 is non-negative when attention is restricted to isothermal response (i.e., θ = constant) and the strain energy function Σ is assumed to be equal to the isothermal Helmholtz free energy. Furthermore, it is noted that the evolution equations for the elastic deformation measures (J, B e ) determine elastic deformations from the zero-stress state of the material given by

T = 0 for J = 1 , B e = I . (21)
In addition, the evolution equations ( 8) and ( 9) are Eulerian in the sense that they are insensitive to arbitrariness of: a reference configuration; an intermediate stress-free configuration; a total deformation measure; and a plastic deformation measure [START_REF] Rubin | Removal of unphysical arbitrariness in constitutive equations for elastically anisotropic nonlinear elastic-viscoplastic solids[END_REF].

Determination of the stress rate

With the help of ( 8), ( 9), ( 10) and ( 18) it can be shown that

J µ Ṫ = (D • I)(kJ 2 d 2 U dJ 2 I - 5 3 B e ) + LB e + B e L T - 2 3 (B e • D) I -2Γg e . ( 22
)
If Γ vanishes then the response is elastic and if Γ is a homogeneous function of order one in the rate D then the response is rate-independent. Otherwise, the response is rate-dependent. Moreover, from (11) it follows that

∂g ∂g e = ( 3κ 2γ 2 e ) g e . (23)
Consequently, for this simple model, the rate of inelasticity in ( 22) is normal to the yield function g in (11).

The standard rate-independent model

For the standard rate-independent model, the derivative of the yield function is given by ġ = ĝ -Γg ,

ĝ = 3κ 4γ 3 e [ 1 3 (B e • I)B e + B e B e - 8γ 2 e 9 I] • D , g = κ γ e + H . (24)
Then, the loading conditions and the consistency condition are given by Γ = 0 for (g < 0) and (g = 0 with ĝ ≤ 0) ,

Γ = 1 1 + H ( 3 4κ 2 )[ 1 3 (B e • I)B e + B e B e - 8κ 2 9 I] • D for (g = 0 with ĝ > 0) .
(25)

The smooth transition model

For the general smooth-transition model [START_REF] Hollenstein | Modeling a smooth elastic-inelastic transition with a strongly objective numerical integrator needing no iteration[END_REF]] the function for Γ is specified by (1) and in this paper, attention is limited to the simplest rate-independent case for metals with a finite elastic range by specifying (3) so that Γ is given by

Γ = b 1 ε g . ( 26
)
Then, with the help of ( 7), ( 22) and ( 26) the stress rate can be expressed in the form

J µ Ṫ = LB e + B e L T + (kJ 2 d 2 U dJ 2 I ⊗ I - 5 3 B e ⊗ I) • D - 2 3 (I ⊗ B e ) • D -2b 1 g (g e ⊗ 2D 3 ε ) • D , (27)
which shows that the tangent stiffness is not symmetric and that it depends nonlinearly on D.

Using either ( 25) or ( 26), the evolution equations ( 8), ( 9) and ( 13) predict rate-independent response since they are homogeneous of order one in rate. Moreover, they predict elastic response when the rate of inelasticity Γ vanishes.

Analysis of the limit load

To study the limit load, consider the normal component df n of an element of force df applied on a material surface with unit normal n and elemental material area da defined by

df = Tnda , df n = n • df = T • (n ⊗ n)da . ( 28
)
Using the results that for a material surface

d(nda) dt = [(D • I) I -L T ]nda , ṅ = [(D • n ⊗ n) I -L T ]n , (29)
it follows that a limit load occurs when the material derivative of df n vanishes

ḋf n = A • (n ⊗ n)da = 0 , A = Ṫ -LT -TL T + [ 4 3 (D • I) + (D • n ⊗ n)]T . ( 30
)
In general, the elastic distortional state of the material B e can be expressed in the form

B e = 1 3 α 1 I + B e , (31)
where the invariant α 1 is defined in (15). It was shown in [START_REF] Rubin | Calculation of hyperelastic response of finitely deformed elastic-viscoplastic materials[END_REF] that for any value of the deviatoric tensor B e the invariant α 1 can be determined by the solution of the cubic equation which requires B e to be unimodular det( 1 3

α 1 I + B e ) = 1 , (32)
with a unique real solution with (α 1 ≥ 3) and (α 1 = 3 for B e = 0). Moreover, the evolution equation for B e is obtained by taking the deviatoric part of ( 9) to obtain

Ḃ e = LB e + B e L T - 2 3 (B e • D) I - 2 3 (D • I)B e -ΓB e , (33)
To analyze the limit load it is convenient to consider the special case of uniform proportional loading with the velocity field v and the uniform velocity gradient L given by

v = Dx , L = D = 1 3 (D • I) I + D . ( 34
)
More specifically, the deviatoric rate of deformation D and the elastic distortional deformation B e are specified by

D = 3 2 ε N , B e = 2 2 3 γ e N , (35)
with N being a unit tensor defined by the constant Lode angle β and the constant orthonormal triad p i , such that

N = 2 3 [cos( π 6 + β)p 1 ⊗ p 1 + sin(β)p 2 ⊗ p 2 -cos( π 6 -β)p 3 ⊗ p 3 ] , N • N = 1 , - π 6 ≤ β ≤ π 6 . ( 36
)
In these expressions, the eigenvalues of D are defined by ε and the eigenvalues of B e are defined by γ e . Also, the maximum eigenvalues of D and B e are associated with p 1 and the smallest eigenvalues are associated with p 3 . Moreover, for β < 0 it is possible to think of extension in the p 1 direction with compression in the transverse p 2p 3 plane and for β > 0 it is possible to think of compression in the p 3 direction with extension in the transverse p 1p 2 plane.

For proportional loading the constitutive equation ( 18) yields

T = -p I + J -1 µB e , Ṫ = -J dp dJ (D • I) I + J -1 µ[-(D • I)B e + Ḃ e ] . (37)
Thus, for uniaxial tension in the p 1 direction with (n = p 1 ) and (β = -π/6) it follows that

N = 2 3 [p 1 ⊗ p 1 - 1 2 (p 2 ⊗ p 2 + p 3 ⊗ p 3 )] , σ = T • (p 1 ⊗ p 1 ) = -p + J -1 µB e • (p 1 ⊗ p 1 ) , T • (p 2 ⊗ p 2 ) = -p + J -1 µB e • (p 2 ⊗ p 2 ) = 0 , (38)
and

Ṫ • (p 2 ⊗ p 2 ) = -J dp dJ (D • I) + J -1 µ[-B e + Ḃ e ] • (p 2 ⊗ p 2 ) = 0 . (39)
Moreover, for the simple case with U in (15) and p in (18) are specified by

U = J -1 -ln(J) , p = kµ( 1 J -1) , (40)
it can be shown using the symbolic program Maple that ( 38) and ( 39) yield

σ = 2µγ e 1 + 2γe 3k , J = 1 + 2γ e 3k , D • I = (α 1 + 2γ e ) ε -2Γγ e 2γ e + 3k . (41)
Also, it can be shown using (30) that the limit load requires

ḋf n = kµ[(9kα 1 + 4α 1 γ e -4γ 2 e ) ε -2Γγ e (9k + 4γ e )]da (3k + 2γ e ) 2 = 0 , (42)
which can be solved to obtain

2Γγ e = α 1 [1 - 4γ 2 e α 1 (9k + 4γ e ) ] ε . ( 43
)
In addition, the solution of (32) for uniaxial tension in the p 1 direction is given by

α 1 = A + 4γ 2 e A , A = 1 2 108 -64γ 3 e + 12 81 -96γ 3 e 1/3
.

(44)

The standard model

For the standard rate-independent model with (γ e = κ) during loading, (25) yields

2Γκ = α 1 ( 1 + 2κ α 1 1 + H ) ε , (45)
so that ( 43) requires

1 + 2κ α 1 1 + H = 1 - 4κ 2 α 1 (9k + 4κ) , (46)
where α 1 (κ) is determined by substituting (γ e = κ) in (44).

The smooth model

For the rate-independent smooth model, Γ is given by ( 26) so for loading (43) requires

γ e -κ = α 1 2b 1 [1 - 4κ 2 α 1 (9k + 4κ) ] , (47)
where α 1 (γ e ) is given by ( 44). Although H does not appear explicitly in this expression, the values of (γ e , κ) are coupled by integration of the evolution equations, which depend on H.

Analysis of localization by perturbing a homogeneous deformation

This section considers small perturbations superimposed on a homogeneous deformation. In all cases, inertia and external body forces are neglected for simplicity. Specifically, consider a deformation characterized by uniform homogeneous fields for which ∂J/∂x = 0 , ∂B e /∂x = 0 , ∂T/∂x = 0 , ∂κ/∂x = 0 , L = 0 . (48)

In particular, this state satisfies equilibrium

divT = ∂T/∂x • I = ∂T/∂X • F -T = 0 , ( 49
)
where X is the location of a material point in a fixed reference configuration and F is the deformation gradient from the reference to the present configuration. Moreover, using the fact that

d dt (F -1 ) = -F -1 L , (50)
the material derivative of (49) yields

d dt (divT) = ∂ Ṫ/∂X•F -T -∂T/∂X•L T F -T = ∂ Ṫ/∂x•I-∂T/∂x•L T . (51)
Consequently, instantaneously the equation for maintaining equilibrium about this homogeneous state reduces to

d dt (divT) = ∂ Ṫ/∂x • I = 0 . (52)
Next, using ( 22) it follows that for the uniform state (48), the instantaneous equilibrium equation ( 52) requires

[kJ 2 d 2 U dJ 2 I - 5 3 B e - 2 3 (B e • I) I] ∂(D • I)/∂x + (∂L/∂x) • B e + B e ( LT ∂L/∂x • I) - 2 3 (B e • ∂D/∂x) -B e ∂Γ/∂x = 0 . ( 53
)
To study potential instability modes, consider a body which has been homogeneously deformed and is at rest, and consider a perturbation with a superposed velocity and associated velocity gradient given by

v = ( γ w ) exp(ξ)m , L = γ exp(ξ)(m⊗n) , D = 1 2 γ exp(ξ)(m⊗n+n⊗m) ξ = wn • x , m • m = 1 , n • n = 1 , γ > 0 , w > 0 , (54)
where the constant unit vector m defines the direction of velocity, the constant unit vector n defines the direction of its dependence on x, γ is a positive constant controlling the magnitude of L and w is a positive constant wave number. Here, attention is limited to a body which is either finite in the direction n or is semi-infinite with ξ ≤ 0 to ensure that the perturbation remains bounded. Using these expressions, it follows that

∂L/∂x = γw exp(ξ)(m⊗n⊗n) , ∂D/∂x = 1 2 γw exp(ξ)(m⊗n+n⊗m)⊗n , ε = 1 √ 3 γ exp(ξ) 1 + 1 3 (m • n) 2 ∂ ε/∂x = 1 √ 3 γw exp(ξ) 1 + 1 3 (m • n) 2 n . ( 55
)
Since the expression w exp(ξ) appears in all terms in (53) it can be factored out of the resulting equation. Also, for inelastic loading the values of Γ developed below are non-zero and for elastic response they are zero.

The standard model

For the standard rate-independent model with Γ specified by ( 25) it can be shown that

Γ = γ exp(ξ)Γ , ∂Γ/∂x = γw exp(ξ)Γn , Γ = 1 1 + H ( 3 4κ 2 )[ 1 3 (B e • I)(B e • m ⊗ n) + (B e B e • m ⊗ n) - 8κ 2 9 (m • n)] . (56)

The smooth model

For the rate-independent smooth model with Γ specified by (26) it can be shown that

Γ = γ exp(ξ)Γ , ∂Γ/∂x = γw exp(ξ)Γ n , Γ = 1 √ 3 b 1 g 1 + 1 3 (m • n) 2 .
(57)

The perturbation vector

For either of the specifications ( 25) or ( 26), equation ( 53) yields a perturbation vector a of the form

a = (m • n)[kJ 2 d 2 U dJ 2 n + 1 9 (B e • I) n - 2 3 B e n] - 2 3 (B e • m ⊗ n) n + 1 3 (B e • I) m + (B e • n ⊗ n) m -ΓB e n = 0 (58)
with Γ given by (56) for the specification (25) and by (57) for the specification (26).

4.3.1. For m = αn with α = ±1 For the case when m is parallel to n, it can be shown that (58) requires

a = α[kJ 2 d 2 U dJ 2 + 4 9 (B e • I) + 1 3 (B e • n ⊗ n)] n -( 2 3 α + Γ)B e n] = 0 . (59)
Next, let s be an arbitrary unit vector normal to n

s • n = 0 , (60)
and rewrite (59) in the component form

a • n = α[kJ 2 d 2 U dJ 2 + 4 9 (B e • I)] -( 1 3 α + Γ)(B e • n ⊗ n)] , (61a)
a • s = -( 2 3 α + Γ)(B e • s ⊗ n) . (61b) 4.3.2. For m × n = 0
For this case is it convenient to define the unit vector

s = m × n |m × n| , (62)
and to express the components of a in the forms

a • n = (m • n)[kJ 2 d 2 U dJ 2 + 4 9 (B e • I) + 1 3 (B e • n ⊗ n)] - 2 3 (B e • m ⊗ n) -Γ(B e • n ⊗ n) = 0 , (63a)
a • m = (m • n)[{kJ 2 d 2 U dJ 2 + 1 9 (B e • I)}(m • n) - 4 3 (B e • m ⊗ n)] + 1 3 (B e • I) + (B e • n ⊗ n) -Γ(B e • m ⊗ n) = 0 , (63b)
a • s = -[ 2 3 (m • n) + Γ](B e • s ⊗ n) = 0 . ( 63c
)
Given a state of the material defined by (J, B e , κ), localization can occur if the scalar equations ( 61) are satisfied by some combination of (n, s, α) or if the scalar equations ( 63) are satisfied by some combination (m, n).

Another analysis of localization

Rice (1976) analyzed localization in plastic deformation by considering a body that experiences homogeneous plastic deformation with uniform velocity gradient L. Using the evolution equation ( 29) for the unit normal n to a material surface and the expression for the traction vector t applied to this material surface

t = Tn , (64)
it follows that ṫ = Ṫn + T ṅ .

(65) [START_REF] Rice | The localization of deformation[END_REF] considered the increment of the rate ∆ ṫ due to an increment ∆L of the form

∆L = g ⊗ n , ( 66
)
where g is a vector to be determined. Rice's localization condition proposes that the material experiences localized plastic flow if ∆ ṫ is stationary

∆ ṫ = 0 . (67)
As a special case, consider a homogeneous uniform material which has been homogeneously deformed, is currently in a state at the elastic-inelastic boundary and is at rest. Then, instantaneous application of the velocity gradient

L = g ⊗ n , (68)
yields a jump in the rate ṫ so that Rice's criterion for this case requires ṫ = 0 for L = g ⊗ n .

(69)

Moreover, using (29) it follows that for this velocity gradient the unit normal m remains constant ṅ = 0 for L = g ⊗ n .

(70) Also, Rice's assumed a piecewise-linear constitutive equation of the form (4). Within the context of the smooth model, the inelastic distortional deformation rate depends nonlinearly on the total deformation rate D. [START_REF] Rice | The localization of deformation[END_REF] explained that nonlinear dependence of L in (4) on Ḟ0 can be handled by considering an unlimited number of models based on different L. This latter approach could include the analysis of the smooth model. Instead, here a homogeneously deformed homogeneous unifom material is considered which is at rest in a state at the elastic-inelastic boundary and is instantaneously subjected to a velocity gradient of the form

L = γm ⊗ n , γ > 0 , m • m = 1 , n • n = 1 , (71)
with the rate of the traction due to (71) given by ( 65). For this velocity gradient it follows from (29) that n remains constant (70). Motivated by the special case (69) of Rice's criterion [START_REF] Rice | The localization of deformation[END_REF] it is assumed that the condition for material instability requires the rate of the traction vector to be stationary ṫ = Ṫn = 0 for L = γm ⊗ n .

(72)

Next, using (71) it can be shown that

D = 1 2 γ(m ⊗ n + n ⊗ m) , D = γ[ 1 2 (m ⊗ n + n ⊗ m) - 1 3 (m • n) I] , ε = 1 √ 3 γ 1 + 1 3 (m • n) 2 . (73)

The standard model

For the standard rate-independent model with Γ specified by (25), L specified by (71) and with use of (73), it can be shown that

Γ = Γ γ , (74)
with Γ given by (56).

The smooth model

For the rate-independent smooth model with Γ specified by ( 26), use can be made of (73) to show that Γ has the form (74) with Γ specified by (57).

Equivalence of the traction rate and perturbation vector criteria

Using (10), ( 54) and ( 56) or (57), equation ( 22) can be written in the form

J µ γ Ṫ = exp(ξ)[(m • n)(kJ 2 d 2 U dJ 2 I - 5 3 B e ) + (m ⊗ n)B e + B e (n ⊗ m) - 2 3 (B e • m ⊗ n) I -ΓB e ] . (75)
Then, for uniform deformation (52) yields

J µ γ div( Ṫ) = w exp(ξ)[(m • n)(kJ 2 d 2 U dJ 2 I - 5 3 B e ) + (m ⊗ n)B e + B e (n ⊗ m) - 2 3 (B e • m ⊗ n) I -ΓB e ] • n , (76)
which with the help of ( 58) can be expressed as

[ J µ γw exp(ξ) ]div( Ṫ) = a . ( 77
)
Consequently, for the velocity field ( 54), the traction-rate criterion (72) is equivalent to the perturbation vector criterion (52).

Small deformation equations

For small deformations, it is convenient to introduce the total strain ε, its deviatoric part ε and the total volumetric stain ε v , such that

ε = 1 3 ε v I + ε , ε • I = 0 . (78)
Also, the equivalent total distortional strain rate ε in ( 7) is defined by

ε = 2 3 ε • ε . ( 79
)
Since there is no inelastic volumetric rate of deformation in the smooth model, the volumetric elastic strain ε ev equals the volumetric total strain ε v and the elastic strain ε e can be expressed as

ε ev = ε v , ε e = 1 3 ε v I + ε e , ε e • I = 0 , (80)
where ε e is the deviatoric part of ε e . The evolution equation for elastic distortional deformation rate ( 9) is approximated by

ε e = ε -Γε e , (81)
the equivalent elastic distortional deformation γ e in ( 11) is specified by

γ e = 3 2 ε e • ε e , (82)
and the evolution equation for the equivalent plastic strain in ε p (12) becomes

εp = 2 3 D p • D p • = 2 3 Γγ e , D p = Γε e . (83)
In addition, the yield function g and evolution equation for hardening κ are specified by ( 11) and ( 13), respectively. The strain energy function (15) and rate of material dissipation ( 17) are given by

ρ 0 Σ = 1 2 µ(kε 2 v + ε e • ε e ) , U = 1 2 ε 2 v , D = T • ε -ρ 0 Σ ≥ 0 . ( 84
)
Then, the constitutive equation for stress becomes

T = µ(kε v I + 2ε e ) , (85)
and the rate of material dissipation reduces to

D = Γµε e • ε e ≥ 0 , (86)
which is automatically satisfied. Moreover, it follows from ( 25), ( 26), ( 78), ( 79), ( 81) and ( 85) that

1 µ Ṫ = k εv I + 2( ε -Γε e ) , (87)
where for the standard model Γ = 0 for (g < 0) and (g = 0 with ĝ ≤ 0) ,

Γ = 1 1 + H ( 3 2κ 2 )[ε e - 4κ 2 9 I] • ε for (g = 0 with ĝ > 0) , (88)
and for the smooth model

Γ = b 1 ε g = b 1 2 3 ε • ε g . ( 89
)
6.1. Monotonic proportional loading For monotonic proportional loading, the total distortional strain rate ε and the elastic distortional strain ε e can be expressed in the forms

ε = ε 3 2 N , ε e = γ e 2 3 N , ε ≥ 0 , ε > 0 , ( 90
)
where the total distortional strain ε is non-negative, its rate ε is positive and N is a constant unit deviatoric tensor characterizing the direction of loading defined in (36). Using these expressions, the evolution equation ( 81) for the elastic distortional strain ε e can be rewritten in the form

γe = 3 2 ε -Γγ e , (91)
which is solved using the initial conditions

ε(0) = 0 , κ(0) = κ 0 . (92)
Specifically, using the evolution equation ( 13), this equation integrates to obtain

γ e = 3 2 ε - κ -κ 0 H . (93)
6.1.1. The standard model Using the approximations

J = 1 + ε v , B e = I + 2ε e , D = ε , (94)
and neglecting higher order terms in small quantities in (25) it follows that Γ is specified by

Γ = 0 for γ e < κ , Γ = 3 ε 2κ(1 + H) for γ e = κ , (95)
so the evolution equations ( 91) for γ e and (13) for κ are given by

dγ e dε = 3 2 , dκ dε = 0 for ε < 2 3 κ 0 , dγ e dε = 3 2 ( H 1 + H ) , dκ dε = 3 2 ( H 1 + H ) for ε ≥ 2 3 κ 0 , (96)
which integrate to yield

γ e = Min( 3 2 ε, κ 0 )+ 3 2 ( H 1 + H ) ε- 2 3 κ 0 , κ = κ 0 + 3 2 ( H 1 + H ) ε- 2 3 κ 0 . (97)
6.1.2. The smooth model Using ( 11) and ( 26), the evolution equations ( 91) for γ e and (13) for κ are given by

dγ e dε = 3 2 , dκ dε = 0 for ε < 2 3 κ 0 , dγ e dε = 3 2 -b 1 (γ e -κ) , dκ dε = Hb 1 (γ e -κ) for ε ≥ 2 3 κ 0 , (98)
which integrate to yield

γ e = 3 2 ε - 3 ε -2 3 κ 0 2(1 + H) + 3 2b 1 (1 + H) 2 [1 -exp{-b 1 (1 + H) ε - 2 3 κ 0 }] , κ = κ 0 + H 3 2 ε -γ e .
(99)

Small deformation limit load for uniaxial tension

For small values of γ e it can be shown that

α 1 = 3 . (100)
6.2.1. Limit load for the standard model Using (100) and neglecting κ relative to unity, it follows that the condition (46) for the limit load in the standard model requires

H = 0 . (101)
Then, using (97) the values (γ eL , κ L , ε L ) of (γ e , κ, ε) for the limit load are given by

γ eL = κ L = κ 0 , ε L = 2 3 κ 0 > 0 , (102)
which corresponds to the onset of inelastic deformation rate.

Limit load for the smooth model

Using (100) it follows that for small deformations the limit load associated with (47) for the smooth model requires

γ e -κ = 3 2b 1 . (103)
Moreover, using (99) it can be shown that

γ e -κ = 3 2b 1 (1 + H) [1 -exp{-b 1 (1 + H) ε - 2 3 κ 0 }] , (104)
so the values (γ eL , κ L , ε L ) of (γ e , κ, ε) for inelastic loading at the limit load are given by

γ eL = κ 0 + 3 2b 1 (1 + H) [1- H ln(-H) (1 + H) ] , κ L = κ 0 - 3H 2b 1 (1 + H) [1+ ln(-H) 1 + H ] , ε L = 2 3 κ 0 - ln(-H) b 1 (1 + H) for ε L - 2 3 κ 0 ≥ 0 . ( 105
)
In contrast with the result (102) of the standard model that predicts a limit load at the onset of inelastic deformation rate, the smooth model predicts a limit load for non-positive values of H in the range

-1 ≤ H ≤ 0 , (106)
with a finite value of ε L only for (H < 0). In this regard, it is noted that even for a constant value of κ the smooth model exhibits apparent hardening. Also, it can be shown the value ε L in (105) is the same as the value of ε which causes a peak value of the function γ e (ε) in (99).

To compare the solutions for different values of (H, b 1 ) it is convenient to solve the first of (105) for the value κ 0 which causes the limit load to occur at a specified value γ eL of elastic distortional deformation. Table 1 records the values (κ 0 , κ L , ε L ) which cause the same limit load for two values of H and different values of b 1 . Figure 1 shows the loading curves of γ e versus ε for the values recorded in Table 1 with symbols indicating the values at the limit load. From these results it can be seen that increasing b 1 and decreasing H both cause ε L to decrease.

Small deformation perturbation vector criterion

Using (84) and the approximations (94), the perturbation vector criteria (63) based on the mode characterized by ( 54) and ( 55) reduce to

a • n = (m • n)[k + 4 3 + 2 3 (ε e • n ⊗ n) - 4 3 (ε e • m ⊗ n)] -2Γ(ε e • n ⊗ n) = 0 , (107a)
a • m = (m • n) 2 (k + 1 3 ) - 8 3 (m • n)(ε e • m ⊗ n) + 1 + 2(ε e • n ⊗ n) -2Γ(ε e • m ⊗ n) = 0 , (107b)
a • s = 2 3 (m • n) + Γ](ε e • s ⊗ n) = 0 , (107c)
where for the standard rate independent model, Γ in ( 56) is approximated by

Γ = 1 1 + H ( 3 2κ 2 )(ε e • m ⊗ n) , (108)
and for the smooth transition model use is made of (11) to rewrite Γ in (57) in the form

Γ = b 1 √ 3 1 - κ γ e 1 + 1 3 (m • n) 2 . ( 109
)
6.3.1. Discussion of the solution procedure Noting that the perturbation vector criterion (107) are independent of the volumetric strain ε v and the constitutive equations do not depend on the Lode angle β, it follows that the stress state for material instability is independent of the pressure. Thus, with regard to material instability, states of compression in the p 3 direction correspond to states of tension in the p 1 direction for the same magnitude |β| of the Lode angle. Consequently, the full range of material instability for the models discussed in this paper can be explored by limiting attention to the range

- π 6 ≤ β ≤ 0 . ( 110
)
For monotonic loading with a specified value of β, the state of the material is parameterized by the value of ε and is determined by the values of (κ 0 , H) in the solution (97) for the standard model and by the values of (b 1 , κ 0 , H) in the solution (99) for the smooth model. For the standard model it is assumed that material instability initiates at the onset of yield when (ε = 2 3 κ 0 , γ e = κ 0 ). For this model, material instability first occurs at the largest negative value of H which admits real solutions of the criterion (107). In contrast, for the smooth model the value of ε at the onset of material instability is not known and must be determined. Specifically, the onset of material instability is characterized by the smallest values of ε in the solution (99) which admits real values solutions of the criterion (107). As for the standard model, the existence of material instability depends on the value of H with the largest negative value of H determining the strongest material that will admit material instability.

Analysis of solutions of the perturbation vector criterion (107)

In general, ε e can be expressed in the form (90). Assuming that (ε e • s ⊗ n = 0), the equation (107c) requires

Γ = - 2 3 (m • n) , (111)
which ensures that (m • n = 0) for nonzero rate of inelasticity. Thus, (107a) reduces to Next, consider the general unit vector

k + 4 3 + 2(ε e • n ⊗ n) - 4 3 (ε e • m ⊗ n) = 0 . ( 112
n = cos(φ n )[-sin(θ n )p 3 + cos(θ n )p 1 ] + sin(φ n )p 2 , (114)
which is defined by the angles (θ n , φ n ). The vectors (m, n, s) form a righthanded triad of unit vectors with s being orthogonal to the plane of (m, n) and with m not necessarily being orthogonal to n. Consequently, it follows from ( 113) that s and m can be defined by

s = b |b| , b = (ε e n)×n , m = cos(θ m -θ n )(n×s)+sin(θ m -θ n )n , ( 115
)
where θ m is another angle to be specified.

Solution for uniaxial extension in the p 1 direction

For uniaxial extension in the p 1 direction the Lode angle is specified by β = -π/6. Due to symmetry in the p 2p 3 plane it follows that there is no physical dependence on the angle φ n so this angle can be set equal to zero. This causes (m, n) in ( 54) to lie in the p 1p 3 plane of maximum shear with (m, n, s) specified by (see Fig. 2)

m = cos(θ m )p 3 + sin(θ m )p 1 , n = -sin(θ n )p 3 + cos(θ n )p 1 , s = p 2 for θ n + π 2 > θ m . ( 116
)
For this case, (107c) is satisfied automatically and it follows from (107a) that Γ is given by

Γ = 1 2(ε e • n ⊗ n) [ 4 3 (ε e • m ⊗ n) -(m • n){k + 4 3 + 2 3 (ε e • n ⊗ n)}] . (117)
Then, substituting ( 117) into (107b) yields the condition

F (θ m , θ n , γ e , β, κ) = a • n = 0 , (118)
with (γ e , β, κ) characterizing the state of the material at localization. For a proper choice of ε v this solution can characterize uniaxial stress in the p 1 direction and for another choice of ε v it can characterize plane strain in the p 1p 3 plane with extension in the p 1 direction. It will be shown in Section 8 that the failure mode for unaxial stress is necking which is a limit load that requires a three-dimensional analysis and is not a material instability mode characterized by the velocity field (54) and the perturbation vector equations ( 107). However, the failure mode for plane strain extension is a material instability mode characterized by (107).

Solution for simple shear in a plane containing p 1

For simple shear in a plane containing p 1 , the Lode angle is specified by β = 0 • . Due to symmetry in the p 2p 3 plane it follows that maximum shear occurs in all planes containing p 1 and there is no physical dependence on the angle φ n so this angle can be set equal to zero. This causes (m, n) in (54) to lie in the p 1p 3 plane of maximum shear with (m, n, s) specified by ( 116). It will be shown in Section 8 that this is a material instability mode characterized by the perturbation vector equations (107).

Examples

For all the examples discussed in this section attention is limited to the small deformation theory with Poisson's ratio specified by

ν = 1 3 . ( 119
)
Also, to study the influence of changes in the material parameters (b 1 , H) in the smooth model a reference material was specified with the inelastic material parameters In particular, the material response due to these changes is normalized by determining the value of κ 0 that causes the peak value γ eL of γ e associated with the limit load analysis in subsection 6.2.2 to be the same as that predicted by ( 120)

b 1 = 1000 , κ 0 = 0.008 , H = -0.3 . ( 120
)
γ eL = 0.0090372 . ( 121
)
The procedure described in section 5 was used to determine the state and values (ε f , γ ef , κ f ) of (ε, κ) and the values of (θ m , θ n ) at the onset of material instability for each set of values (κ 0 , b 1 , H). Also, the value of ε L associated with the limit load was determined for each set of material parameters.

7.1. Material instability results for extension (β = -π/6) For extension in the p 1 direction the Lode angle is specified by

β = - π 6 , (122)
and the results are recorded in Table 2. The onsets of material instability (determined by the values ε f ) are marked with symbols. Also, the value H = -0.15 is close to the maximum negative value of H which admits real values of the angles (θ m , θ n ). It is interesting to note that all 9 materials listed in Table 2 have nearly identical values for the angles (θ m ≈ 50 • , θ n ≈ 44 • ) for material instability.

For comparison, the value Γ in (108) was used to solve the equations ( 117) and ( 118) for the standard model to determine the maximum negative value of H which admits real solutions for (θ m , θ n ) at the onset of material instability. Specifically, the value of κ 0 was specified to be equal to γ eL in (121) to obtain

κ 0 = γ eL = 0.0090372 , H = -0.22418 , θ m = 48.400 • , θ n = 42.024 • , (123)
which are values close to the values of (θ m , θ n ) in Table 2 predicted by the smooth model.

It is important to emphasize that the smooth model predicts material instability for a range of values of H and in particular it predicts material instability for larger negative values of H than the standard model, which means it predicts material instability for less softening.

Also, the values recorded in Tables 1 and2 for (H = -0.15) indicate that the strains (ε L , ε f ) predicted by the the smooth model for the limit load and for material instability, respectively, increase with decreasing values of b 1 for the same peak load γ eL . In addition, it is noted that the values of ε f are larger than those of ε L which indicates that material instability occurs after the limit load had occured. For simple shear in the p 1p 3 plane of maximum shear the Lode angle is specified by

β = 0 • , (124)
and the results are recorded in Table 2. The onsets of material instability are marked with symbols. Also, the value H = 0 is the maximum non-negative value of H that admits real values of the angles (θ m , θ n ). Again, it is interesting to note that all 9 materials listed in Table 2 have nearly identical values for the angles (θ m ≈ 50 • , θ n ≈ 44 • ) at the onset of material instability. Also, it is observed from Fig. 5 and Fig. 6 that material instability occurs very close the the peak load with the points at the onset of material instability indicated for H = 0 ocurring before ε equals infinity due to the numerical convergence criterion.

For comparison, the value Γ in (108) was used to solve the equations ( 117) and ( 118) for the standard model to determine the maximum negative value of H that admits real solutions for (θ m , θ n ) at localization. Specifically, the value of κ 0 was specified to be equal to γ eL in (121) to obtain

κ 0 = γ eL = 0.0090372 , H = 0 , θ m = 45.078 • , θ n = 44.957 • , (125)
which are values close to the values of (θ m , θ n ) in Table 2 predicted by the smooth model.

For this problem, it is again noted that the smooth model predicts material instability for a range of values of H and the maximum value of H = 0 is the same for both the smooth and standard models.

Finite element validations

The aim of the finite element simulations provided in this section is to validate the analytical predictions for instability and localization. Specifically, these simulations validate that the analysis predicts the relevant instability or localization models. 1 for the corresponding values of κ 0 ).

Numerical implementation

The small strain version of the rate-independent smooth model discussed in Section 5 has been implemented in the implicit finite element code Zset (Z-set package, 2013). The numerical integration of the constitutive equations is performed using an explicit second order Runge-Kutta method with automatic time stepping and an implicit Newton method involving the consistent tangent matrix of the model. In the case of convergence problems which frequently occur when simulating strain localization phenomena, the Newton integration can be switched to the Runge-Kutta method at each Gauss point of the finite element mesh.

The simulations presented in this work are performed for a softening material with modulus H = -0.15 and for three values of parameter b 1 = 500, 1000, 1500, as considered in the analytical results of Tables 1 to 3. For each value of b 1 a different value of κ 0 is chosen such that the maximum stress is the same for all tensile curves, according to the values indicated in Tables 1 to 3. To illustrate the overstress effect in the smooth model, tensile loading and unloading curves are presented in Fig. 7. It is apparent that lower values of b 1 lead to an increase of the overstress, as expected. The nonlinear part of the unloading branch is clearly visible and corresponds to the decrease and vanishing of the overstress. After reloading, the apparent yield stress is smaller than the value before the onset of unloading.

Finite element simulations of strain localization in rectangular plates are performed in the next section for various loading conditions. The rectangular plate is characterized by a ratio of length to width equal to 2.5. The plate is modeled with a geometrical imperfection which is an invisible V-notch corresponding to the translation in the direction 1 (along the plate's width direction) of one boundary node by 2.5 × 10 -4 times the plate's width. This defect induces a negligible reduction in the cross-section by this same amount. The presence of this defect ensures triggering of a localization band starting always at the same location for simulations involving a softening material.

The localization simulations in this section are the first ones based on the smooth elastoplastic model proposed by [START_REF] Hollenstein | Modeling a smooth elastic-inelastic transition with a strongly objective numerical integrator needing no iteration[END_REF]. Finite element simulations of localization phenomena in the absence of internal length (no regularization) are still useful to detect the localization point and also the orientation of the localization bands. Even if the band thickness is mesh-dependent, the prediction of the orientation of the localization band is quite robust as demonstrated for an anisotropic elastoplastic media in [START_REF] Forest | Strain localization in single crystals: Effect of boundaries and interfaces[END_REF]. Both the predicted localization point and the orientation can therefore be compared to the analytical predictions. In the absence of regularization, there is no point in using finer meshes than those used in the Fig. 8 to 11.

The detection of instability modes in finite element simulations (as well as in experiments) can be based on several criteria. The spurious meshdependence observed in simulations is regarded as an indicator of localization occurence. A reliable quantitative indicator for the onset of localization, at least for some suitable loading conditions, is the occurrence of unloading events outside the localization band [START_REF] Petryk | Material instabilities in elastic and plastic solids[END_REF].

In the following simulations, a localization band is said to develop when the total strain increases inside the band at the same time that it decreases in the surrounding material. To detect the formation of such bands an indicator of unloading is defined as the sign of the mechanical power of T • ε of the deviatoric stress T . This indicator is computed at each time step of the simulations and used to determine the first unloading event and the time step for which a full localization band has formed. 1 for the corresponding values of κ 0 ).

Plane strain extension of a rectangular plate

The formation of a localization band during plane strain extension is shown in Fig. 8. The strain is concentrated into a one to two GP (Gauss point) thick band in which it continues increasing whereas unloading is observed in the red zone outside of the band. Before the onset of material instability, the strain state of the plate is almost homogeneous. The equivalent elastic strain vs. axial strain curves at a GP close to the imperfection are exhibited in Fig. 9. The states of the GP at the onset of unloading in the mesh and for the fully developed band are indicated with symbols. They provide an interval of detection of material instability in the FE simulations. Each simulated curve stops when convergence can no longer be achieved in the finite element simulation due to strong localization.

Figure 9 shows that the theoretical predictions somewhat overestimate the total strain at the onset of material instability. The difference between the theoretical and finite element predictions decreases with increasing values of b 1 , (i.e. for decreasing overstress values). In the limit of large b 1 , the classical results obtained for the standard model are retrieved. This suggests that the existence of the overstress postpones the onset of material instability 2 for the other simulation parameters).

compared to the reference standard case. Using a finer finite element mesh leads to a thinner band and slightly earlier occurrence of the onset of material instability (not presented here). Irrespective of the mesh size the orientation of the band is always 45 • with respect to the tensile axis. This value is close to the theoretical value of the angle θ n recorded in Table 2.

Simple shear of an infinite plate

For the simulation of simple shear, a square plate with a central defect is considered. It is subjected to periodic boundary conditions in order to avoid boundary effects in shear. A 2D mesh is used under plane strain conditions. The material of the central element has a value κ 0 which is 0.015% smaller than the surrounding elements. The defect is nearly indiscernable but is sufficient to trigger localization starting from this location. The mean total shear strain applied to the plate increases monotonically and the strain field remains quasi-homogeneous until material instability occurs causing the two orthogonal shear bands shown in Fig. 10. The shear bands are exactly one element thick and are oriented at 45 • from the principal shear axes, in exact agreement with the predicted value of θ n in Table 3. The curves in Fig. 11 show that the theoretical predictions of the elastic strain γ ef at the onset of material instability are in excellent agreement with the numerical simulations. However, the theoretical predictions of the total strain ε f at the onset of material instability are a bit smaller than those of the numerical simulations, with this difference becoming smaller with increasing b 1 . The theoretical predictions are taken from Table 3. It is well-known that shear loading is prone to localization, which is predicted to take place at the maximum of the shear curve.

Uniaxial tension of a rectangular plate

The case of uniaxial tension of a rectangular plate must be treated using a fully 3D model. It is well-known that the 2D plane stress model is physically unsound because it ignores the compatibility conditions with respect to the out-of-plane direction. As a result, the regularizing effect of the out-of-plane necking of the thin sheet on shear banding can be observed only in the 3D case [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]. Plane stress conditions artificially promote strain localization. The 3D plate of Fig. 12 has the same in-plane dimensions and the same geometrical imperfection as in section 8.2. Four quadratic elements through the thickness are considered for the computational study of localization in a thin plate in tension. It has been checked that the loading curve up to the onset of localization is only slightly dependent on the number of elements through the thickness and the in-plane mesh size. However, refinement of the in-plane mesh size significantly influences post-failure response as a result of shear band formation. The formation of a localization band inclined at about 54 • from the tensile axis is shown in Fig. 12. This band orientation is classically observed during thin sheet deformation for standard plasticity models using the von Mises criterion. It is also found here for the smooth model. In addition it is noted that the localization mode characterized by the velocity field (54), the perturbation vector equations ( 107) and the values of θ m , θ n in Table 2 is not consistent with the necking mode predicted by the numerical simulation.

Due to the combination of in-plane shear banding and out-of-plane necking modes, the limit load criterion derived in section 6.2.2 is better suited to predict the onset of localization in uniaxial tension. This coincides with the well-known Considère criterion classically used in tension to detect the onset of necking in tensile specimens. The tensile curves obtained for three values of the b 1 parameters and H = -0.15 are shown in Fig. 13. The predicted limit points are indicated on the curves according to Table 1, together with numerically detected localization points. It is found that band formation in the FE simulation occurs only slightly after the predicted limit loads. This demonstrates the consistency of the limit load analysis with the theoretical failure criterion developed for the smooth model.

Conclusions

The rate-independent form of the smooth transition elastic-inelastic model developed in [START_REF] Hollenstein | Modeling a smooth elastic-inelastic transition with a strongly objective numerical integrator needing no iteration[END_REF]] differs from the standard rateindependent model in that it depends nonlinearly on the direction of the total deformation rate. Consequently, analysis of modes models of material instability for the smooth model need special attention. The nonlinear equations were used to develop an expression for the maximum normal force that can be applied to a material surface as a condition for the limit load. Also, a perturbation of an arbitrary homogeneous nonlinear state of stress was used to develop an expression for a perturbation vector that must vanish to maintain equilibrium. In addition, a special case of the material instability criterion developed in [START_REF] Rice | The localization of deformation[END_REF], which requires the rate of traction applied to a material surface to be stationary for an applied shearing rate was developed and shown to yield the same condition that the perturbation vector vanishes.

The small deformation equations for the smooth model were used to obtain analytical results for the limit load and perturbation vector conditions for example problems. The smooth model exhibits an overstress that increases with decreasing values of the parameter b 1 and it exhibits softening, which tends to cause localization, with more rapid softening occurring for smaller values of the parameter H. The evolution equations of the small deformation smooth model were solved for monotonic proportional loading to obtain expressions for stress-strain curves that allow for analytical solutions of the conditions for material instability. To study the influence of changes in b 1 , H the initial value κ 0 of the hardening parameter was determined so that the limit load occurs at the same specified value γ eL of elastic strain for each pair b 1 , H. The values κ L of the hardening variable and ε L of the total strain at the limit point for each pair b 1 , H are recorded in Table 1.

The main results of this Table are that the limit strain ε L increases with decreasing values of b 1 and increasing values of H.

In addition, the small deformation equations were programmed as a constitutive model in a finite element code to obtain numerical solutions for these example problems. Numerical simulations of example problems with H = -0.15 were considered to examine the influence of the parameter b 1 . Qualitatively, the theoretical prediction of the dependence of localization on b 1 is the same as that in the numerical simulations for all examples. The numerical and theoretical values γ ef of the elastic strain and ε f of the total strain at the onset of material instability were determined. These values of γ ef are in excellent agreement with the theoretical values and the differences between the numerical and theoretical values of ε f decrease with increasing values of b 1 .

Plane strain extension

The theoretical and numerical predictions of the angle θ n between the normal to the shearing surface and the direction of extension are in very good agreement. However, the theoretical values of the total failure strain ε f at the onset of material instability are larger than those of the numerical simulations.

Simple shear

The theoretical and numerical predictions of the angle θ n are in exact agreement. However, the theoretical values of the total failure strain ε f at the onset of material instability are a bit smaller than those of the numerical simulations.

Uniaxial tension of a 3D thin plate

The numerical simulation of a thin plate requires a 3D formulation to capture the out-of-plane necking mode which is better predicted by the limit load criterion than the perturbation vector criterion based on a shearing mode of material instability. The numerical simulations predict the observed value of θ n . Also, the numerical and theoretical values of the elastic strain ε ef at the limit load are in excellent agreement.

One of the main conclusions of this work is that the smooth model predicts a tendency for the total strain at the onset of localization to be delayed by increase in the overstress.

Figure 1 :

 1 Figure 1: Unaxial tension: Loading curves predicted by the smooth model for different values of H showing the influence of b 1 . The limit load points are marked with symbols.

Figure 2 :

 2 Fig. 1 Sketch of velocity direction m and the normal n to the failure plane

Figure 3 :

 3 Figure 3: Loading curves for extension in the p 1 direction with (β = -π/6) showing the influence of b 1 for different values of H. The onsets of material instability are marked with symbols.

  Figure 3 presents zoomed in figures near the peak of the γ eε curves showing the influence of changes in b 1 for three values of H. Figure 4 presents zoomed in figures near the peak of the γ eε curves showing the influence of changes in H for three values of b 1 .

Figure 4 :

 4 Figure 4: Loading curves for extension in the p 1 direction with (β = -π/6) showing the influence of H for different values of b 1 . The onsets of material instability are marked with symbols.

Figure 5 :

 5 Figure 5: Loading curves for simple shear (β = 0) in the p 1p 3 plane of maximum shear showing the influence of H for different values of b 1 . The onsets of material instability are marked with symbols.

  Figure 5 presents zoomed in figures near the peak of the γ eε curves showing the influence of changes in b 1 for three values of H. Figure 6 presents zoomed in figures near the peak of the γ eε curves showing the influence of changes in H for three values of b 1 .

Figure 6 :

 6 Figure 6: Loading curves for simple shear (β = 0) in the p 1p 3 plane of maximum shear showing the influence of b 1 . The onsets of material instability are marked with symbols.

Figure 7 :

 7 Figure 7: Loading and unloading curves for the smooth model in uniaxial tension for H = -0.15 and three values of parameter b 1 (seeTable 1 for the corresponding values of κ 0 ).

Figure 8 :

 8 Figure8: Simulation of localization under plane strain extension loading conditions: Finite element mesh and contour map of the ε-field (left) and of the unloading indicator (red: unloading, blue: continued loading) at the last simulation step (right). The results are obtained for H = -0.15 and b 1 = 500 (see Table1for the corresponding values of κ 0 ).

Figure 9 :

 9 Figure 9: Simulation of localization under plane strain tension loading conditions: response curves at the integration point (GP) close to the imperfection. The points identifying the localization interval in the numerical simulation are denoted by filled symbols and the theoretical localization points (ε f , γ ef ) from Table 2 are denoted by open symbols. The results are obtained for H = -0.15 and three values of b 1 .

Figure 10 :

 10 Figure 10: Simulation of shear localization: Finite element mesh and contour map of the ε-field (left) and of the unloading indicator (red: unloading, blue: continued loading) at the last simulation step (right). The results are obtained for H = -0.15 and b 1 = 500 (see Table2for the other simulation parameters).

Figure 11 :

 11 Figure 11: Simulation of shear localization: response curves at the integration point (GP) close to the imperfection. The points identifying the localization interval in the numerical simulation are denoted by filled symbols and the theoretical localization points (ε f , γ ef ) from Table 2 are denoted by open symbols. The results are obtained for H = -0.15 and three values of b 1 .

Figure 12 :

 12 Figure 12: Simulation of strain localization in a 3D thin plate in uniaxial tension with 4 quadratic elements in the thickness: The finite element mesh and contour map of the ε-field are shown together with a detailed picture of the imperfection zone. The results are shown for H = -0.15 and b 1 = 1500.

Figure 13 :

 13 Figure13: Simulation of localization in uniaxial tension in a 3D thin plate using four quadratic elements in the thickness: response curves at the integration point (GP) close to the imperfection. The points identifying the onset of material instability in the numerical simulation are denoted by filled symbols and the theoretical limit load points (ε L , γ eL ) from Table1 aredenoted by open symbols. The results are obtained for H = -0.15 and three values of b 1 .

Table1: Small deformation theory: List of the parameters for the limit load predicted by the smooth model for uniaxial stress.

Acknowledgements

This research was partially supported by MB Rubin's Gerard Swope Chair in Mechanics. MB Rubin would also like to acknowledge Mines ParisTech for hosting him in November 2017 during part of his sabbatical leave from Technion.

Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

500 0.66894e-2 0.80818e-2 0.13530e-1 0.83527e-2 0.48981e-2 49.550 44.389 1000 0.78633e-2 0.74741e-2 0.97824e-2 0.86938e-2 0.69663e-2 49.552 44.386 1500 0.82546e-2 0.69910e-2 0.85310e-2 0.88080e-2 0.76563e-2 49.557 44.393 H = -0.3 500 0.69629e-2 0.80818e-2 0.93304e-2 0.88848e-2 0.54297e-2 49.559 44.395 1000 0.80000e-2 0.70533e-2 0.76777e-2 0.89609e-2 0.72333e-2 49.573 44.422 1500 0.83458e-2 0.67105e-2 0.71268e-2 0.89864e-2 0.78347e-2 49.565 44.406 H = -0.9 500 0.74845e-2 0.70969e-2 0.74370e-2 0.89983e-2 0.55431e-2 49.573 44.422 1000 0.82609e-2 0.65606e-2 0.67310e-2 0.90178e-2 0.72901e-2 49.684 44.643 1500 0.85196e-2 0.63822e-2 0.64956e-2 0.90242e-2 0.78723e-2 49.746 44.768

Table 2: Small deformation plane strain uniaxial extension in the p 1 direction: List of the parameters and the failure values ε f , γ ef , κ f , θ m , θ n at the first point of localization. The value of the total strain ε L associated with the maximum value γ eL = 0.0090372 for the limit load criterion is also included.

500 0.60372e-2 0.23830e-1 0.90370e-2 0.60372e-2 45.225 45.251 1000 0.75372e-2 0.14928e-1 0.90371e-2 0.75372e-2 45.223 45.249 1500 0.80372e-2 0.11960e-1 0.90371e-2 0.80372e-2 45.229 45.259 H = -0.15 500 0.66894e-2 0.89228e-2 0.90372e-2 0.60373e-2 45.180 45.160 1000 0.78633e-2 0.74738e-2 0.90372e-2 0.75373e-2 45.180 45.160 1500 0.82546e-2 0.69908e-2 0.90372e-2 0.80372e-2 45.273 45.347 H = -0.9 500 0.74845e-2 0.70968e-2 0.90372e-2 0.60373e-2 45.186 45.173 1000 0.82609e-2 0.65608e-2 0.90372e-2 0.75373e-2 45.267 45.334 1500 0.85196e-2 0.63822e-2 0.90372e-2 0.80372e-2 45.431 45.663

Table 3: Small deformation simple shear: List of the parameters and the failure values ε f , γ ef , κ f , θ m , θ n at the first point of localization.