NURBS plasticity: non-associated plastic flow

This paper extends the non-uniform rational basis spline (NURBS) plasticity framework of Coombs et al. (2016) and Coombs and Ghaffari Motlagh (2017) to include non-associated plastic flow. The NURBS plasticity approach allows any smooth isotropic yield envelope to be represented by a NURBS surface whilst the numerical algorithm (and code) remains unchanged. This paper provides the full theoretical and algorithmic basis of the non-associated NURBS plasticity approach and demonstrates the predictive capability of the plasticity framework using both small and large deformation problems. Wherever possible errors associated with the constitutive formulation are specified analytically and if not numerical analyses provide this information. The rate equations within the plasticity framework are integrated using an efficient and stable implicit stress update algorithm which allows for the derivation of the algorithmic consistent tangent which ensures optimum convergence of the global out of balance force residual when used in boundary value simulations. The important extension provided by this paper is that the evolution of plastic strain is decoupled from the yield surface normal. This allows the framework to model more realistic material behaviour, particularly in the case of frictional plasticity models where an associated flow rule is known to significantly overestimate volumetric dilation leading to spurious results. This paper therefore opens the door for the NURBS plasticity formulation to be used for a far wider class of material behaviour than is currently possible. c ⃝ 2018 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons. org/licenses/by/4.0/).


Introduction
Constitutive models that provide incremental relationships between stress and strain are essential for boundary value analysis of engineering problems. Within this, one of the most common classes of material behaviour is elastoplasticity where the elastic region of stress space is bounded by a yield surface. On this yield surface the material will undergo elasto-plastic material behaviour and countless yield envelopes have been proposed since the works of Tresca [1] and von Mises [2]. However, the form of the yield function impacts on the stress integration algorithm which is required to convert the rate form of the plasticity equations into an incremental form that can be used in boundary value simulations (using the finite element method for example). This issue was overcome by the non-uniform rational basis spline (NURBS) plasticity framework of Coombs et al. [3] and Coombs and Ghaffari Motlagh [4] which allowed any smooth isotropic yield envelope to be integrated using the same numerical algorithm. However, [3] and [4] were limited to the case of associated plastic flow where the form of the yield surface governs both the yielding of the material and the evolution of plastic strains. This limits the form of material behaviour that can be predicted and, in the case of frictional plasticity models, leads to a significant overestimation of volumetric dilation. This paper overcomes this limitation by decoupling the evolution of plastic strains from the yield surface normal leading to a non-associated plastic flow constitutive framework.
In this paper we do not attempt to review all of the constitutive models available in the literature, instead an interested reader is referred to the work of Yu [5] for a general review of constitutive models. In the specific area of NURBS plasticity, beyond the work of [3,4], the only other paper that the authors are aware of is that of Coelho et al. [6] who construct NURBS response surfaces in biaxial strain and stress space based on curve fitting to experimental data. This is quite different to the approach of [3,4], and that advocated in this paper, where a NURBS yield envelope is constructed and then used within a conventional plasticity formulation. This approach provides a constitutive formulation which is valid for generalised, six-component, stress and strain space.
The NURBS plasticity formulation is combined with an implicit predictor-corrector stress integration algorithm [7] to provide an incremental relationship between stress and strain. Several papers have compared different stress integration algorithms and, as with the vast array of constitutive models, in this paper we do not attempt to review them. In this case the interested reader is refereed to the works of Anandarajah [8] and Safaei et al. [9], amongst others. The reasons for adopting an implicit algorithm in this paper are twofold: (i) they rigorously enforce the consistency conditions at the updated stress state (and in the case of NURBS plasticity, throughout the process) and (ii) allow for the derivation of the algorithmic consistent tangent that ensures asymptotic quadratic convergence of the global residual when used within a boundary value simulation.
The layout of the paper is as follows, Section 2 provides the theoretical framework for hardening non-associated flow NURBS-based plasticity, including the definition of the NURBS surface and the non-associated flow rule, isotropic hardening through the movement of control points, the form of stress integration used and the technique of energy mapped stress that allows us to interpret the stress integration method as a geometric projection. Section 3 provides details on the numerical implementation including the backward Euler (bE) stress integration process and the algorithmic consistent tangent. Numerical examples are presented in Section 4 and, finally, conclusions are drawn in Section 5.
The majority of the paper is presented in tensor form using index notation, the notable exception is the numerics that are presented in matrix-vector form for ease of implementation. Due to the geometric nature of the method presented in this manuscript, the majority of the paper is presented in principal stress and strain space with the following ordering of the principal stresses with tensile stresses taken as positive. Note that although the equations are presented in principal stress space we can do this without loss of generality of the final result as the principal quantities are simply transferred back to generalised quantities at the end of the algorithm. Generalised, 6-component, stress and strain quantities are denoted using( ·).

Non-associated flow NURBS plasticity
This section provides the essential equations required for an isotropically hardening NURBS surface with nonassociated plastic flow. There is significant overlap between the theory presented here and that of Coombs et al. [3] (for perfect plasticity with associated flow) and Coombs and Ghaffari Motlagh [4] (for isotropic hardening with associated flow), however the repetition is retained for the sake of clarity and to provide a self-contained formulation. For more detailed information on the construction of NURBS-based surfaces see the work of Piegl and Wayne [10] and the paper of Coombs et al. [3] for the particular case of perfect plasticity yield envelopes.
A general NURBS surface can be expressed as where k is the physical index, C k are the control point positions and n and m are the number of control points in the local ξ and η directions. 1 The NURBS basis functions, R i, j , are given by where N i, p and N j,q are the pth and qth-degree B-spline basis functions (see [10,11], amongst others), ξ and η are the local positions within the knot vectors and w i, j are the weights associated with the control points.

NURBS-based yield envelopes
Starting from the equation for a NURBS surface (1), a NURBS-based yield envelope [3,4] can be expressed as where (S, σ ) i is the surface outward normal (that is, the partial derivative of S with respect to stress) and σ i the principal stress state. The yield surface separates stress space into two regions: an elastic region where f < 0 and an inadmissible region where f > 0. The boundary between these two regions ( f = 0) is used to identify the point of yielding and stress states on this surface will undergo elasto-plastic deformation. The outward normal to the yield envelope can be obtained through the cross product of the two local derivatives where ϵ i jk is the Levi-Civita tensor. 2 See [10] for efficient algorithms for the calculation of the derivatives of the NURBS surface with respect to the local coordinates. Fig. 1(i) shows a bi-quadratic spherical NURBS surface with the form f = σ i σ i − r 2 y = 0, where r y is the radius of the yield envelope. The control points used to define the surface are shown by the red points and different surfaces can be obtained by moving the positions of the control points and/or modifying the basis functions, R i, j . Fig. 1(ii) also shows a number of outward normals, (S, σ ) i , to the spherical yield surface.
The number of control points required to define a NURBS plasticity yield surface will depend on the form of the yield function. The number of control points needed for some widely used yield envelopes are: • von Mises: open ended cylinder, quadratic NURBS curve in the deviatoric direction to define the circular cross section using 8 control points combined with a straight line in the hydrostatic direction using a minimum of 2 control points -16 control points in total; • Tresca: prismatic regular hexagon, linear NURBS curve in the deviatoric direction to define the hexagonal cross section using 6 control points combined with a straight line in the hydrostatic direction using a minimum of 2 control points -12 control points in total; • Drucker-Prager frictional cone: quadratic NURBS curve in the deviatoric direction to define the circular cross section using 8 control points combined with a straight line in the hydrostatic direction using a minimum of 2 control points -16 control points in total; and • Mohr-Coulomb frictional cone: linear NURBS curve in the deviatoric direction to define the irregular hexagonal cross section using 6 control points combined with a straight line in the hydrostatic direction using a minimum of 2 control points -12 control points in total. However, it may be necessary to use higher order NURBS surfaces in order to be able to avoid numerical issues with the stress integration algorithm, especially in the case of yield envelopes with corners [3,4] (see Section 4.1 for more details for the Drucker-Prager (D-P) and Mohr-Coulomb (M-C) yield surfaces). It is also possible to reduce the number of control points by only defining the yield surface over a single sextant of stress space, where σ 1 ≥ σ 2 ≥ σ 3 , as is done in this manuscript and in the work of Coombs and Ghaffari Motlagh [4]. For example, the number of control points required to define a spherical yield surface in all sextants of stress space is 45 whereas if it is defined in a single sextant only 15 points are required [4]. For associated flow plasticity theory the outward normal to the yield surface also provides the flow direction controlling the evolution of plastic strains, that iṡ whereγ is the scalar plastic multiplier rate (or consistency parameter). This plastic multiplier rate must satisfy the Kuhn-Tucker-Karush consistency conditions where the yield envelope is a function of plastic strain, ε p i , to allow for hardening/softening of the surface. These conditions enforce that the material must either be on the yield surface undergoing elasto-plastic deformation ( f = 0 andγ ≥ 0) or inside the yield surface with purely elastic behaviour ( f ≤ 0 andγ = 0).

Non-associated flow
In the case of non-associated flow the evolution of plastic strains is decoupled from the spatial gradient of the yield envelope. The plastic strains evolve according tȯ where (g, σ ) i is the gradient of the plastic potential surface. Previous attempts to extend the NURBS plasticity framework to include non-associated flow have directly specified the flow direction at the control points rather than specifying the geometry of a plastic potential surface (see Coombs [12]). However, when specifying the flow direction at the control points, it is only possible to recover associated flow over the entire yield surface if there is no coupling between the local knot coordinates (ξ, η). One example is the prismatic von Mises yield surface where the deviatoric stress is ρ = ∑ 3 k=1 σ k and ρ y is the deviatoric radius of the yield surface. When defined using NURBS plasticity, this surface combines a circle in the ξ (deviatoric) direction with a line in the η (hydrostatic) direction. In this case the radius of the yield surface is not dependent on the η value. However, if we introduce a hydrostatic dependency and extend the von Mises surface to that of D-P [13], with the form the deviatoric and hydrostatic directions are coupled and we lose the ability to recover associated plastic flow over the entire yield surface. In (9) ζ = 1 √ 3 ∑ 3 k=1 σ k is the hydrostatic stress, β = tan(φ) is the opening angle of the cone, φ is the friction angle, ζ a = c √ 3 cot(φ) is the location of the cone's tensile apex and c is the material's cohesion. Before moving onto the proposed formulation, the error associated with defining the direction of plastic flow at the control points is investigated for the D-P yield surface with a friction angle of φ = π/9 (20 • ) and a cohesion of 0 Pa. The NURBS yield surface was defined over the hydrostatic range of ζ ∈ [−2, −1]Pa using the following knot vectors with a linear and quadratic bases in the hydrostatic (local η) and deviatoric (local ξ ) directions, respectively. The direction of plastic flow was defined using the approach of Coombs [12] as where (G, σ ) k are the flow directions at the control points and the NURBS basis functions R g i, j (ξ, η) are calculated in the same way as (2). The flow direction was defined using the following knot vectors Ξ g ξ = {0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4} and Ξ g η = {0, 0, 1, 1}, and weights w g ξ = {1, 1, 1, 1, 1, 1, 1, 1, 1} and w g η = {1, 1}, where the weights in the ξ direction have been modified to exactly recover associated plastic flow for a von Mises yield surface. As with the yield surface geometry, the flow direction has linear and quadratic bases in the hydrostatic (local η) and deviatoric (local ξ ) directions, respectively.
The control points positions and associated flow directions for the D-P surface were arranged as shown in Fig. 2. In all cases the flow directions at the control points had a radial deviatoric direction with the appropriate hydrostatic component for the specified friction angle such that the vector was normal to the cone in the hydrostatic direction (as shown in Fig. 2(i)). Fig. 3(i) shows the variation in the error of the non-associativity of the plastic flow direction with the deviatoric local position on the yield surface, where the error is defined as and ∥(·)∥ denotes the L2-norm of (·). The plastic flow direction was calculated using (10) and the normal to the yield surface using (4). Note that Fig. 3(i) only shows the variation of the error in the ξ direction as the error does not vary in the hydrostatic local direction, η. There is zero error (recovering associated plastic flow) when the control points  coincide with the yield surface whereas the error is maximum at the control point local positions that do not coincide with the surface. This error is due to the coupling between the ξ and η directions and increases with an increasing angle of friction, φ, as shown in Fig. 3(ii). In all cases the maximum error is located at ξ = 0.5, 1.5, 2.5 and 3.5 and reduces to zero as φ → 0.
The error in the flow direction is due to the presence of the cross product in (4) meaning that it is not possible to specify the flow direction at the control points and guarantee the recovery of associated plastic flow. That is, even when attempting to recover associated flow, in general where and (·) refers to ξ or η, as required. It is only possible to recover associated flow over the entire yield surface when the form of the yield function in the ξ and η directions are independent. For example, for the cylindrical von Mises surface the yield radius (in the ξ direction, for example) does not depend on η parametric coordinate but this is not the case for the D-P yield envelope where the yield radius varies with η.
In this paper we adopt a more conventional approach, in terms of plasticity theory, and specify the geometry of the plastic potential surface using control points. In the case of associated plastic flow the plastic potential control point positions coincide with those used to define the geometry of the yield surface. The gradient of the plastic potential surface therefore given by where g k (ξ, η) is the plastic potential surface and G k are control point coordinates that control the shape of the surface. Note that it is not necessary to have the same basis functions for the direction of plastic flow as used to describe the geometry of the yield surface however it is assumed that a single set of control points control the form of both the yield, S k , and plastic potential, g k , surfaces.

Isotropic hardening
Consistent with the work of Coombs and Ghaffari Motlagh [4], in this paper the yield surface is allowed to expand or contract isotropically according to the level of inelastic straining at a material point, such that where the superscript (·) 0 denotes the original control point coordinates and h(ε p i ) controls the expansion/contraction of the control points. For linear isotropic hardening the incremental evolution of h is given by where is the value of the hardening function from the previously converged state and α controls the hardening/softening rate and equals zero in the case of perfect plasticity. Initially the hardening parameter is taken to be unity, that is h 0 = 1.

Stress integration
As with the other NURBS plasticity papers [3,4], we use an implicit elastic predictor, plastic corrector scheme. The starting point for the stress integration algorithm is the previously converged, or initial, stress state, σ n i and hardening parameter, h n . This state is then subjected to a strain increment, ∆ε i , giving an elastic trial stress state of D e i j is the linear elastic stiffness matrix and (ε e n ) j is the elastic strain from the previously converged, or initial, state. This trial stress state should be checked against the yield criteria for the material under consideration and if it exceeds to yield envelope then it must be corrected back onto an appropriate stress state on the yield surface. This return stress state, σ r i , is given by and ∆ϵ p j = ∆γ (g, σ ) j is the plastic strain increment and ∆γ the incremental plastic multiplier. Once the incremental plastic multiplier has been obtained the updated elastic strain is given by and the updated hardening parameter, h, from (17). The remaining question is how to obtain the incremental hardening parameter, ∆γ , or the return stress state, σ r i , such that the other quantities can be determined? In this paper we adopt a closest point projection (CPP) implicit stress integration algorithm to arrive at the updated stress state. This corresponds to the minimisation [14] of with respect to the return stress σ r i whilst satisfying the Kuhn-Tucker-Karush consistency conditions (6), where C e i j = (D e i j ) −1 is the elastic compliance matrix. An important point which is often overlooked it that despite this process being referred to as a CPP, the return stress is not the closest point geometrically in standard stress space, but rather the stress that minimises (21). The only case where the CPP is the geometric closest point on the yield surface from the trial stress state is when the Poisson's ratio is zero and C e i j = E −1 δ i j , where E is the Young's modulus and δ i j the Kronecker delta tensor. 3 However, in this paper we adopt the non-associated flow version of energy-mapped ς i space (EMSS) [15,16] to reduce this CPP minimisation to the problem of finding the point on the yield envelope that the normal to the plastic potential surface passes through when intersecting with a trial point outside of the surface. The use of EMSS converts (21) essentially removing the influence of Poisson's ratio, ν, when finding the energy-mapped return stress state ς r i . The following transformation converts between conventional and energy-mapped stress space This mapping leads to a squashing and a stretching of the yield surface in the hydrostatic and deviatoric directions respectively, as shown in Fig. 4(i) for a spherical yield surface with ν = 0.2 and ν = 0.4. Once the closest point solution in EMSS has been found, the solution can be transformed back to conventional stress space. For a NURBS yield surface we only need to map the control point coordinates for the yield, (C k ), and plastic potential, (G k ), surfaces into EMSS, the rest of the NURBS information remains unchanged.

Numerical implementation
Consistent with the perfect plasticity associated plastic flow implementation of Coombs et al. [3], here we use a coarse initial subdivision algorithm to provide the initial starting point for a bE implicit stress integration process. This is to provide an initial estimate for the local positions within the knot vectors, ξ and η in (3) that act as the primary unknowns in the CPP problem (in addition to the updated hardening parameter). The rest of the section is focused on the implicit stress integration algorithm and the derivation of the algorithmic consistent tangent. However, before focusing on the stress integration algorithm we first propose a new method to determine if the material is undergoing elastic or elasto-plastic deformation.

Yield function estimate
Determining the return stress state on the NURBS surface is only required for stress states undergoing elasto-plastic deformation. However, the yield function for NURBS plasticity (3) requires the return stress state and associated outward normal, therefore in previous NURBS plasticity approaches even for elastic stress states the appropriate return state on the yield surface must be determined. In this paper we propose an alternative technique to predict if the material is undergoing elasto-plastic deformation without determining the return stress state. The distance between the trial point and a point on the NURBS surface in the direction of the yield surface outward normal, (S, σ ) i , is In the case of a convex surface, if the trial stress state, σ t i , is inside or on the yield surface then k d ≥ 0 for all points, S i (ξ, η), on the NURBS yield envelope. However, if k d < 0 for any (ξ, η) then the trial stress state is outside of the yield surface. It is not possible to check that all points on the NURBS surface. Instead a parametric grid of (ξ, η) points are checked where the locations are evenly distributed within that local knot coordinates, that is and where ∆ξ = (ξ max − ξ min )/(n d − 1), ∆η = (η max − η min )/(n d − 1) and n d is the number of points considered in each direction. These trial knot locations are shown in Fig. 4(ii) by the white-filled circles in a single sextant of principal stress space for a spherical yield surface. The trial stress state outside of the yield surface is also shown by the black-filled circle as are the vectors used in (24) for two of the trial yield surface points.

Subdivision algorithm
In this paper we adopt an implicit bE algorithm to determine the updated local position on the yield surface, (ξ, η), and new size of the surface, controlled by h. In order to start of the Newton process of the bE algorithm an initial estimate of the return position and the hardening parameter are required. Consistent with the work of Coombs et al. [3], a coarse subdivision process is used for the initial estimate of ξ and η whereas, consistent with [4], the previous value of the hardening parameter provides the starting point for h.
The trial local position on the frozen yield surface (that is, h = h n ) for the implicit algorithm is selected finding the minimum for an initial set of candidate points on the yield surface, where ς S i is the energy-mapped yield surface. See [3] for details of the subdivision algorithm used to determine the candidate points on the yield surface.

Implicit stress integration
The EMSS bE stress integration algorithm contains three unknowns which satisfy the following residuals The tangent directions to the plastic potential surface, ( ς g, ξ ) i and ( ς g, η ) i , are obtained through taking the derivative of the energy mapped plastic potential surface with respect to the local coordinates where (·) is ξ or η depending on the required derivative. The first two residuals in (27) ensure that the return path in EMSS is in the direction of the energy mapped plastic flow rule and the third that the hardening function has reached a stationary value at the update stress state, wherẽ In the above equation ∆σ is the plastic stress increment over the stress return path. The unknowns are updated through a standard Newton-Raphson process where k denotes the Newton-Raphson iteration number. The derivatives required for the Jacobian matrix, [∂r/∂ x], are provided in Appendix B. As with the algorithm for associated flow perfect plasticity, the stress return path for Newton-Raphson procedure described in this paper starts and remains in the yield envelope. This leads to inherently stable numerics as it eliminates one of the issues associated with bE stress integration, namely the form of the yield function outside of the yield surface and its impact on the efficiency and stability of the stress return process. For example, inappropriately constructed yield equations can contain local minima and/or spurious auxiliary surfaces outside of the true yield envelope which can trap a returning stress state.

Algorithmic consistent tangent
Specification of the algorithmic consistent tangent allows for asymptotic quadratic convergence of the global out of balance force residual [17]. The tangent is first constructed in principal stress space and then transformed into six-component stress space using the principal directions associated with the trial elastic strain state (see Appendix A for details). Following the approach of [18], we can linearise the constitutive model into the following form where the derivatives are determined at the updated stress state. ∆h is the incremental form of the hardening function, and in the case of linear isotropic hardening The derivative of the flow direction with respect to stress is The derivatives of the local knot coordinates with respect to stress can be obtained from the inversion of the Jacobian matrix linking the local NURBS coordinates with the principal stress directions, that is The normal to the NURBS surface, (S, σ ) i provides the third direction, orthogonal to the tangent vectors, in the transformation. The derivatives associated with the hardening function are The derivative of the yield function with respect to h is The plasticity model presented in this paper has been expressed in principal components. It is therefore necessary to detail the mapping of this model into generalised 6-component space. The equations to convert the stress and strain measures into 6-component space are given in Appendix A, the appendix also details the transformation of the algorithmic (or elasto-plastic) stiffness matrix and the shear components of the algorithmic consistent tangent stiffness matrix. A pseudo-code for the isotropically hardening non-associated flow NURBS plasticity model is given in Fig. 5 that details fully the steps required in calculating the updated stress state and hardening parameter.

Numerical simulations
This section provides material point and boundary value simulations to demonstrate the performance of the constitutive model and the numerical stress integration algorithm described in the previous sections. However, before the numerical examples are presented we first need to define the NURBS plasticity models.

Non-associated flow plasticity models
This section provides the NURBS and control point information for the two plasticity models used in the numerical analysis presented in this paper, namely the D-P and M-C yield surfaces.

Drucker-Prager
The majority of the NURBS information for the D-P yield surface has already been given in Section 2.2. However, the D-P yield surface (9) contains a tensile apex which poses an issue for the stress return algorithm presented in this paper as the derivatives of the NURBS surface are undefined at this point. Here we follow the same approach as Coombs et al. [3] and Coombs and Ghaffari Motlagh [4] and locally round the apex, as shown in Fig. 6 with ζ a = 0. The yield surface is shown in both hydrostatic versus deviatoric stress space and principal stress space for both the original and rounded surfaces. The knot vector and associated control point weights for the grey curve shown in Fig. 6(i) are This knot and weight vector replaces those defined in Section 2.2 for the η direction. The radius of the rounding curve can be obtained from the point where the true and rounded curves depart, point C in Fig. 6(i), that is where ζ C < ζ a and the arc angle is ϑ = π/2 − φ. The hydrostatic locations of points D and E can be subsequently obtained from ζ C and R, where point D lies on the intersection between the original yield curve and a line of constant hydrostatic pressure from the tensile limit of the rounded surface. It is clear from Fig. 6(i) that introducing rounding at the apex of the yield surface generates a region were stress states should be undergoing elastic behaviour or located on the true yield surface that are actually outside of the NURBS surface and therefore inadmissible. However, the maximum hydrostatic error due to the rounding is [4] error = (ζ a − ζ C ) ) , which can be controlled through setting an appropriate value of ζ C . It is also important to note that the degree of rounding shown in Fig. 6 is for visualisation purposes only and is in excess of that used in the numerical analyses. The ξ -direction knot vectors and control point weights for both the yield and the plastic potential surfaces defined in a single sextant of stress space are Ξ ξ = {0, 0, 0, 1, 1, 1} and where the η-direction information is given above. Combining these gives the thick grey and black dashed curves shown in Fig. 6, where β g = tan(ψ) is the opening angle of the plastic potential surface and ψ ∈ [0, φ] is the dilation angle.

Mohr-Coulomb
The M-C yield surface can be defined as c and φ are the cohesion and the friction angle of the material, respectively. In this paper the M-C yield surface is represented using a bi-quadratic NURBS surface with a locally rounded tensile apex using the same approach as for the D-P surface. The M-C yield surface also is only C 0 continuous at both the compression and extension meridians and this will cause issues for the stress integration algorithm which requires the first derivative of the yield surface with respect to stress. We therefore introduce local rounding into the yield surface around these two meridians. Fig. 7(i) shows a deviatoric section through the NURBS M-C yield surface (thick grey solid line) and the true M-C yield surface (thin black dashed line) in one sextant of stress space. In the figure, the deviatoric radius has been normalised by the yield radius on the compression meridian and the normalised radius on the extension meridian is given bȳ in the case of φ = 20 • ,ρ e = 0.795. Both of the meridians of the NURBS surface have been rounded by the same length, δ, where θ e and θ c are the arc angles for the extension and compression meridians, given by and with A = √ρ 2 e −ρ e + 1. The deviatoric section is constructed using seven control points as shown in Fig. 7(i) and the knot vectors for the yield surface defined in a single sextant of stress space are  Fig. 7(ii), where, as with the D-P surface, the rounding has been exaggerated for visualisation purposes. The maximum normalised deviatoric radius error at the compression or extension meridian is given by where the subscript (·) can be e or c for the extension of compression meridians, respectively. Fig. 8(i) shows the variation of the normalised error versus local knot position for φ = 0 • , 10 • , 20 • , 30 • and 40 • . In all cases the maximum error is located at the compression meridian (ξ = 3) and in the special case of zero friction angle both the compression and extension meridians have the same maximum error. In this case the M-C surface reduces to the Tresca yield surface; a regular prismatic hexagonal yield surface aligned with the hydrostatic axis. As the friction angle increases the error at the compression meridian increases whereas the error at the extension meridian reduces, this is clearly shown in Fig. 8(ii) where the error at the compression/extension meridian is reported for φ ∈ [0, 40] degrees.

Material point investigations
Before analysing some boundary problems we will first investigate the performance of the non-associated plastic flow NURBS plasticity framework at a material point (or stress-strain) level. In all cases two subdivisions (see Coombs et al. [3] for details) were applied before the bE stress return algorithm and n d was set to 5 for the yield function estimation resulting in 25 points being checked on the yield surface to determine if the stress state was undergoing elastic or elasto-plastic behaviour.

Drucker-Prager error analysis
This section analyses the errors associated with the implicit stress integration algorithm for a perfect plasticity D-P yield surface with non-associated flow. The material had a Young's modulus of 100 Pa and a Poisson's ratio of 0.2. The cohesion was set to 0.49 Pa, the friction angle to φ = π/9 (20 • ) and the final 0.1 Pa of the yield surface apex was rounded. 4 The errors associated with the stress return algorithm were evaluated using dilation angles of π/18 and π/36 (10 • and 5 • ). The stress state was initially located on the shear meridian in one of the sextants of stress space with a hydrostatic stress of ζ = 0 kPa. This point was then subjected to a stress increment that took the trial stress state outside of the yield envelope. The space of trial states explored was ρ t /ρ n ∈ [1,6], where the t and n subscripts denote the trial and starting locations.
The errors associated with the trial state are shown in Fig. 9, using the following normalised error measure where {σ NURBS } is the stress return location associated with the NURBS model and {σ e } is the exact stress return. 5 The errors associated with ψ = π/18 and ψ = π/36 are shown on the right and left of the thick black line, respectively. The starting point for ψ = π/18 is shown by the white-shaded circle whereas the red-shaded circle is the starting point for the ψ = π/36 analysis. Although errors of almost 20% are present in the model, exactly the same level of errors are observed in the D-P yield surface integrated with a conventional implicit stress integration procedure. As expected with any predictorcorrection stress integration algorithm, the error increases as the tangential proportion of the stress increment increases. Errors also increase with increasing non-associativity, with ψ = π/18 having a maximum error of 1.66 × 10 −1 whereas for the ψ = π/36 the maximum error was 1.92 × 10 −1 , again this is due to the return path having a larger tangential component relative to the yield surface normal direction.

Mohr-Coulomb error analysis
This section analyses the errors associated with the implicit stress integration algorithm for a perfect plasticity M-C yield surface with associated and non-associated flow. The material had a Young's modulus of 100 Pa and a Poisson's ratio of 0.2. The cohesion was set to 0.49 Pa and the friction angle to 20 • and the final 0.1 Pa of the yield surface apex was rounded. Fig. 10 shows the errors in the NURBS stress return algorithm for the case of associated plastic flow. The M-C implementation of Clausen et al. [19] which explicitly includes the corners at the compression Fig. 9. Stress return error analysis for D-P with non-associated flow with ψ = π/18 and ψ = π/36. The inner circle shows a deviatoric section through the yield surface and the white and red filled circles the starting points for the ψ = π/18 and ψ = π/36 error analyses, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) and shear meridians and the hydrostatic apex was taken as the reference solution for the stress integration algorithm. Three different cases were examined, namely the comparison of a rounded NURBS surface with: [19] where the strain increment was applied in a single step; (B) δ = 0.1 with [19] where the strain increment was applied in 1000 sub steps; and (C) δ = 0.01 with [19] where the strain increment was applied in 1000 sub steps.
In all cases the strain increment was applied to the NURBS plasticity model in a single step. The errors associated with the stress return algorithm, as evaluated by (43), are shown over two sextants for the three cases in Fig. 10. The thick black and red lines show a section through the yield surface for δ = 0.1 and δ = 0.01, respectively. In all cases the stress state started at the intersection of the shear meridian and the yield surface with zero hydrostatic stress, as shown by the identified stress states in the figure.
When comparing the true and rounded surfaces when subjected to the strain increment in a single step (region A in Fig. 10) the non-zero errors exist in the trial regions that return to the rounded part of the yield surface. As expected from (42), the larger errors are associated with the compression meridian and symmetry is observed in the error distribution as the starting point is irrelevant when applying the strain increment in a single step. The maximum error for any trial point was 8.02 × 10 −2 .
Subdividing the strain increment into 1000 sub steps with the same level of local rounding (region B in Fig. 10) gives a maximum error of 1.59 × 10 −1 and this maximum error is located close to the extension meridian in the adjacent sextant to the initial stress state (the σ 1 > σ 3 > σ 2 sextant in Fig. 10). Interestingly reducing the length of the meridian rounding actually increases the maximum error in the stress return, with region C having a maximum error of 1.86 × 10 −1 . This is because more points get trapped in the corner region that should actually return to the planar part of the yield surface. However, reducing the local rounding length does reduce the errors in the sextant of the initial stress state, as shown by the σ 3 > σ 1 > σ 2 sextant in Fig. 10. Table 1 gives the maximum errors for the M-C model with associated (grey shaded row) and non-associated plastic flow. As with the D-P model, the maximum error increases with increasing non-associativity and with decreasing δ. The M-C yield surface with sharp corners (δ = 0) had the largest error for the reason explained above. The maximum number of N-R iterations for any of the trial stress states on any of the NURBS rounded M-C yield surfaces to find the updated stress state was 5.

Boundary value simulations
This section presents the results from boundary value simulations to demonstrate the performance of the proposed non-associated flow NURBS plasticity framework. In all cases two subdivisions (see Coombs et al. [3] for details) were applied before the bE stress return algorithm and n d was set to 5 for the yield function estimation resulting in 25 points being checked on the yield surface to determine if the stress state was undergoing elastic or elasto-plastic behaviour.

Plane strain rigid footing
The first analysis is that of a 1 m wide plane strain rigid footing displacing into a weightless 10 m by 5 m domain with a Young's modulus of E = 1 × 10 7 kPa and a Poisson's ratio of ν = 0.48. The yielding of the material was governed by a perfect plasticity D-P yield envelope (9) with cohesion of c = 490 kPa and a friction angle of θ = π/9 (20 • ). Dilation angles of 20 • and 10 • were considered with the second case requiring a non-associated flow rule. Table 1 Maximum stress return errors for M-C with non-associated plastic flow, where n is the number of sub-steps and the letters A, B and C correspond to the same combinations of meridian rounding and sub-steps as used in Fig. 10. Fig. 11. Rigid footing: normalised pressure versus displacement response with a D-P constitutive model. The problem was analysed using a mesh comprising of 135 eight-noded bi-quadratic quadrilateral elements integrated using reduced four-point quadrature. Due to symmetry only half of the problem was modelled. The mesh is the same as that used by [3,[20][21][22][23], amongst others. A vertical displacement of 4 mm was applied to the footing over 20 equal loadsteps.
The normalised pressure versus displacement response is shown in Fig. 11 along with an inset figure showing the mesh detail around the rigid footing. B is the footing width, p is the footing pressure and v is the vertical displacement of the footing. The associate flow (AF) NURBS plasticity response (long dashed line) is compared with the result of de Souza Neto et al. [20] (discrete points) and that of a conventional bE CPP implementation of the D-P yield surface (thick grey line). Excellent agreement is seen between the three results.
The non-associated flow result with ψ = π/18 (10 • ) is also presented in Fig. 11 for both the conventional bE (thick grey line) and the NURBS (short dashed line) models. As with the associated flow results, excellent agreement is seen between the two models. Fig. 12 shows the deformed mesh around the footing for associated (top) and nonassociated (bottom) plastic flow where the mesh has been shaded by the vertical displacement with dark grey being the maximum downwards displacement. The original mesh is shown by the fine dashed line and the displacements have been exaggerated by ×20. The associated flow case exhibits excessive volumetric dilation in the region adjacent to the footing leading to unrealistic heaving of the ground surface. Reducing the dilation angle from π/9 (20 • ) to π/18 (10 • ) significantly reduces the heave leading to a more realistic deformed surface profile. This rigid footing problem was also analysed using the M-C yield function. Fig. 13 gives the normalised pressure versus displacement response for a perfect plasticity M-C model with cohesion of c = 490 kPa and a friction angle of θ = π/9 (20 • ) for both associated and non-associated (ψ = π/18) plastic flow. The same mesh and elastic material properties used for the D-P analysis were used but in this case a displacement of 3mm was imposed using 20 equal displacement-controlled loadsteps. The results using the NURBS plasticity framework (discrete points) are presented alongside the implicit implementation of Clausen et al. [19] (solid grey lines). In both cases excellent agreement is seen between the two plasticity approaches and the normalised limit pressures agree well with the analytical limit pressure of p/c = 14.84, as shown by the dashed black line.
The rigid footing problem is now analysed using a mixture of D-P and M-C yield surfaces. The same mesh and elastic material properties as used in the previous analyses were used in this case and a displacement of 3 mm was imposed using 20 equal displacement-controlled loadsteps. The different material models were randomly assigned at an element level and the distribution of the different constitutive models is shown in Fig. 14 where the white and grey elements represent the D-P and M-C models, respectively. Only the NURBS information changed in order to simulate the different materials and the underlying numerics for the stress integration was the same in all cases. Both the M-C and D-P yield surfaces had a cohesion of c = 490 kPa and a friction angle of θ = π/9 (20 • ). Fig. 15 shows the normalised pressure versus displacement response for this distribution of constitutive models with associated and non-associated (ψ = π/18) plastic flow. The responses when all of the elements used associated flow M-C and D-P models are also shown using dashed lines. The mixed constitutive model associated flow simulation has the potential to predict a response in the grey-shaded region of the figure depending on the number, and the distribution, of the elements using the different material models. As with the previous analyses, the NURBS plasticity response, shown by the discrete points, is compared with a conventional bE implementation and excellent agreement is seen for both for the associated flow and non-associated flow simulations.

Finite deformation cavity expansion
This section presents the finite deformation expansion of a cylindrical cavity under internal pressure. The problem was modelled using an updated Lagrangian finite element formulation with the following assumptions: 1. multiplicative decomposition of the deformation gradient into elastic and plastic components; 2. linear behaviour between Kirchhoff stress and elastic logarithmic strain; and 3. exponential map of the plastic flow.  The combination of these assumptions allows the infinitesimal strain format of stress integration algorithms to be maintained whilst extending the boundary value analysis to include finite deformation mechanics. See Coombs [18] for details of the large deformation formulation adopted in this paper.
The problem was modelled using a two-dimensional axisymmetric finite element code and had an inner radius of 1 m and a fixed outer radius of 2 km (to approximate an infinite domain). The inner radius, a 0 , was expanded to 5 m using 20 equal displacement-controlled loadsteps. The perfect plasticity M-C material had the following material parameters: Young's modulus of 100 MPa, Poisson's ratio of 0.3, cohesion of 70 kPa and a friction angle of π/9 (20 • ). Dilation angles of π/9 (20 • ), π/18 (10 • ) and 0 were considered. 150 fully-integrated eight-noded elements were used  to analyse the problem and the size of the elements were progressively increased by a factor 1.1 from the inner to the outer surface (the inner and outer elements had radial lengths of 0.124 mm and 182 m, respectively). The analytical solution for this problem, assuming perfect plasticity M-C constitutive behaviour, was given by Yu and Houlsby [24]. Fig. 16 shows the normalised internal pressure versus normalised radius response for the three dilation angles, where p is the internal pressure and a the current radius. The response of the NURBS plasticity model is shown by discrete points with the dilation angles of ψ = 20 • , 10 • and 0 • being represented using white, grey and black filled circles respectively. The analytical solution given by Yu and Houlsby [24] is shown by the dashed black lines. The NURBS plasticity framework shows good agreement with the analytical solution for all of the dilation angles. Table 2 gives the global Newton-Raphson (N-R) residual for loadsteps 1, 2, 3, 19 and 20 with ψ = π/18. The global tolerance was 1 × 10 −06 and the maximum number of iterations for any loadstep was 5. The data presented in the table shows quadratic (or near quadratic) convergence of the global out of balance force, demonstrating the correct implementation of the algorithmic consistent tangent for large deformation elasto-plasticity.

Conclusions
This paper has presented for the first time an extension of the NURBS plasticity framework of Coombs et al. [3] to include non-associated plastic flow. Information regarding the shape and size of the yield and plastic potential surfaces is stored at control points covering one sextant of stress space. This allows any smooth isotropic convex yield surface to be included within the plasticity formulation without modification of the numerics used for the stress integration. Within this, the implicit backward Euler stress update procedure contains three unknowns, namely the local position of the updated stress state on the yield surface and the updated hardening parameter. Therefore, introducing nonassociated flow does not increase the number of unknowns in the stress integration algorithm as compared to the associated flow NURBS hardening formulation of Coombs and Ghaffari Motlagh [4]. In addition, when solving for these unknowns the current estimate for the updated stress state remains on the yield envelope and satisfies the consistency condition throughout the process. This results in a stable stress integration algorithm as, unlike in conventional closest point projection schemes, the stress state cannot be trapped in a local minimum or converge to an axillary surface outside of the true yield envelope (see [22] for more information on the impact of the form of the yield function on the stability of the bE stress integration process).
Unlike previous attempts to extend the NURBS plasticity framework to non-associated flow, in this paper the geometry of the plastic potential surface is represented at the control points. This guarantees the ability to recover associated plastic flow for any surface. This is not possible when specifying the flow direction at the control points due to the presence of a cross product in the yield surface normal. This paper has also proposed a new method to estimate if a trial stress state in inside (elastic) or outside (elastoplastic) of the yield surface. This reduces the computational cost of the constitutive model as the closest point to the yield envelope only needs to be determined for those stress states undergoing elasto-plastic material behaviour.
The numerical examples included in this paper have quantified the errors associated with the stress integration process for both the Drucker-Prager and Mohr-Coulomb yield envelopes and where possible the maximum errors associated with yield surface smoothing have been specified analytically. The performance of the plasticity framework has also been demonstrated through a number of small strain and finite deformation finite element analyses where analytical or existing numerical results exist. Within this, the correct derivation and implementation of the algorithmic tangent has been confirmed through the convergence rate of the global out of balance force residual.
The key extension provided by this paper is that the evolution of plastic strain is decoupled from the yield surface normal. This allows the plasticity formulation to model a more diverse range of material behaviour and, in the case of frictional plasticity, predict a more realistic volumetric response, avoiding the excessive dilation seen with associated plastic flow. where (·) denotes the six-component stress and strain quantities. The transformation matrix is given by (q 1 ) 2 (q 2 ) 2 (q 3 ) 2 q 1 q 2 q 2 q 3 q 3 q 1 (q 4 ) 2 (q 5 ) 2 (q 6 ) 2 q 4 q 5 q 5 q 6 q 6 q 4 (q 7 ) 2 (q 8 ) 2 (q 9 ) 2 q 7 q 8 q 8 q 9 q 9 q 7 2q 1 q 4 2q 2 q 5 2q 3 q 6 q 1 q 5 + q 4 q 2 q 2 q 6 + q 5 q 3 q 3 q 4 + q 6 q 1 2q 4 q 7 2q 5 q 8 2q 6 q 9 q 4 q 8 + q 7 q 5 q 5 q 9 + q 8 q 6 q 6 q 7 + q 9 q 4 2q 7 q 1 2q 8 q 2 2q 9 q 3 q 7 q 2 + q 1 q 8 q 8 q 3 + q 2 q 9 q 9 q 1 + q 3 q 7 Note the relative magnitudes of the residuals in (27), with the first two residuals will be of the order of stress squared and the third close to unity. Normalising the first two residuals (and the appropriate entries in the Jacobian matrix) with respect to Young's modulus reduces the potential for ill conditioning of the Jacobian matrix.