A rate-independent crystal plasticity algorithm based on the interior point method

.


Introduction
Constitutive models based on crystal plasticity are widely used in Finite Element analysis for investigating the affects of microstructural parameters on the macroscopic mechanical behavior.Many different implementations of time integration algorithms based on crystal plasticity models exist which might be broadly classified into two main categories, namely visco-plastic and rate-independent models.
Crystal plasticity, as defined by the Schmid law, is based on the understanding that when a certain critical shear stress is reached on a slip plane and a slip direction, i.e. a slip system, dislocations become mobilized and start to move in that specific direction resulting in plastic deformation [1].Due to the symmetry of the crystal lattice however the number of possible slip systems exceed the number of independent plastic strain rate tensor components [2,3].This causes numerical problems for multi-surface crystal plasticity due to linear dependence of the system of equations resulting from the flow rule [4][5][6].
Motivated by observations that there is a relation between the exerted shear stress and the rate of slip, visco-plastic crystal plasticity algorithms utilize the ratio of the resolved shear stress to the critical resolved shear stress to obtain a flow rule [7,8].This approach inherently regularizes the system of equations and slip may be allowed to occur on every slip system for all stress states [9] or only for the active set [4].In [4,5] the correspondence of the algorithmic relation between the physically motivated visco-plastic modeling of crystal plasticity and penalty regularization of multi-surface plasticity is demonstrated.As the penalty factor increases the solution becomes less and less rate dependent sacrificing numerical stability due to ill-conditioning.There exists rate independent solutions without resorting to multi-surface plasticity such as [10,11].
In the pure rate-independent approach the flow rule is determined by the consistency condition which requires that for slip to occur the resolved shear stress must be equal to the critical resolved shear stress on a slip system and it must remain equal as E-mail address: e.s.perdahcioglu@utwente.nl.https://doi.org/10.1016/j.cma.2023.116533Received 19 June 2023; Received in revised form 1 October 2023; Accepted 9 October 2023 long as plastic deformation continues [3].This implies that when on a certain slip system the critical resolved shear stress is not reached yet the slip rate must vanish.Algorithmically this can be treated as an optimization problem where the plastic dissipation is maximized [5,12].
Rate-independent models must be able to handle linearly dependent equations that result from the underlying physics of multisurface plasticity.Numerical problems might occur due to two types of singularity.When multiple slip systems are exactly at the same state (i.e.identical hardening behavior) and have exactly the same resolved shear stress there is no unique solution as to how slip will be partitioned.The other possibility, which might be more commonly encountered, is related to the Taylor ambiguity mentioned above.When more than 5 slip systems become simultaneously active the set of slip rates on these systems loose their uniqueness as different combinations will lead to the same plastic strain rate tensor.This will also result in singularity as any of the possible solutions will lead to a theoretically admissible result.When this happens numerical techniques can be utilized to choose one of the admissible solutions.Singular value decomposition [3,4,13] is one of the possible numerical approaches and in [4] as well as [6] a more detailed treatment of the problem can be found.In [5] an augmented Lagrangian approach is presented which exploits the regularization behavior of the penalty method in obtaining a rate-independent solution.
Rate-independent algorithms that utilize the concept of the active set of slip systems might also suffer from the fact that during iterations the active set might need to be changed.For instance when the closest point projection algorithm is used [13][14][15] a slip system that was inactive in the trial set might become active after convergence.It is equally possible that an active slip system becomes inactive as the stress state and its magnitude changes during iterations.In this case the increment must be repeated with an updated active set until convergence is obtained and there remains no violations of slip conditions.
An alternative algorithm is based on the interior point method which also considers crystal plasticity as an optimization problem.A variation of this, infeasible primal-dual interior point method, has been previously shown to yield a robust algorithm for crystal plasticity for small strain formulations [12].The robustness of the interior point method follows as a result of two important properties which are very relevant in crystal plasticity.First is that all the iterations are carried out in the feasible domain, meaning that there is no violation of the yield functions that might lead to divergence due to changes in the active set.Secondly the algorithm is inherently regularized against linear dependence of the constraints.
In the current study the interior point method is adopted for large deformation crystal plasticity implementation.Together with the closest point projection algorithm the algorithms are derived based on the maximization of plastic dissipation within a backward Euler time integration scheme.The kinematics of crystal plasticity is treated in rate form for straightforward implementation to existing Finite Element software and algorithmic tangents for both algorithms are derived.The motivation of this study therefore is to investigate the suitability of the interior point method to crystal plasticity in a large deformation setting with the potential of having a robust rate independent algorithm.To this end, aspects of the implementation such as the handling of Taylor ambiguity and large rotations are discussed.The influence of the interior point parameter on the accuracy of the results and the efficiency of the simulations is investigated.Finally the performance of the proposed algorithm is demonstrated with two Finite Element simulations of a volume element where each grain uses the developed constitutive model.

Crystal plasticity formulation
The crystal plasticity framework that the developed algorithms are based upon can be described as a combination of kinematics, flow-rule and hardening.Kinematics defines the decomposition of the total deformation into elastic and plastic parts.Below, the rate form of the decomposition that is used in the algorithm is described starting from the multiplicative decomposition of the total deformation gradient.The flow-rule that describes the evolution of plastic deformation in relation with the stress is derived using the principal of maximum plastic dissipation.Finally the hardening of the material due to slip activity that dictates the relation between the critically resolved shear stress as a function of the slip activity on all slip systems is described.

Kinematics
In crystal plasticity it is commonly accepted that plastic slip occurs on slip systems related to the lattice structure of crystals, i.e.Schmid rule.These slip systems are defined as a combination of unit vectors of slip direction ( () ) and slip plane normal ( () ) for a given slip system .Slip activity that is associated with dislocation motion will result in a simple shear type of deformation that leaves the lattice unrotated.There can exist also elastic deformation that results in a change of stress.The total deformation that is a combination of elastic and plastic parts can be written commonly in the following multiplicative form where  p describes the plastic deformation and Fe describes the elastic counterpart.The elastic deformation relates an intermediate configuration that is defined by the plastic deformation to the current configuration.Any rigid rotation of the material point will also be described as a part of the elastic deformation gradient.Starting from this decomposition the total velocity gradient  can be obtained The plastic velocity gradient tensor is given by the sum of plastic slip rates ( ̇() ) on slip systems defined in the current configuration as [16,17] where  is the number of slip systems and The velocity gradient given in Eq. ( 2) can be decomposed first into a sum of deformation rate,  and spin tensor  (symmetric and skew-symmetric parts, respectively) which in turn can be decomposed as elastic and plastic as Considering Eq. ( 3) it can be seen that for Lp to be isochoric ŝ and m must remain orthogonal, i.e. tr( Lp ) = 0. Following the common assumption that the elastic stretches are small [3,[18][19][20], the slip vectors,  and , are transformed to the deformed configuration using only the rotation resulting from the polar decomposition of Fe , see Eq. (A.8).Since this transformation preserves orthogonality of these vectors it can be concluded that plastic deformation as defined in Eq. (3) preserves volume.Another consequence of this is that during deformation the initially unit vectors  and  remain unit after transformation, ŝ and m, i.e.

∀𝛼 ∶ |𝐬
The elastic part of the deformation rate dictates the change of stress in the material point while the plastic part results in dissipation due to slip activity, as described in the following section.The plastic part of the spin tensor is responsible for the rotation of the lattice with respect to the material which gives rise to texture evolution in polycrystals.
A more detailed treatment of the above relations is given in Appendix A.

Flow rule
In the rate independent formulation, slip occurs only on the slip systems where the resolved shear stress ( () r ) of the slip system  is equal to its slip resistance ( ()  f ).Hence, we can define  () for each slip system as: where  represents the vector with accumulated slip amounts for all slip systems.The systems for which the equality in Eq. ( 7) holds are called active set of slip systems.In this study, the slip rates ̇() for all systems are defined to be positive which necessitates consideration of positive and negative slip directions as separate slip systems [4].
The Cauchy stress is determined by an elastic constitutive definition that might stem from a hyperelastic as well as a hypoelastic formulation provided that it is written in rate form.The rate formulation is preferred since it does not require numerically approximating the incremental deformation gradients.The rate form also leads to a straightforward procedure in determining the algorithmic tangent.In this respect the presented algorithm does not necessitate a particular choice of elastic behavior definition.
A hypoelastic choice is demonstrated due to the fact that in most plasticity applications the elastic stretches remain small and therefore nonlinear effects are negligible whereas rotations are significantly large and of concern.In most commercial Finite Element software the elastic behavior in the global frame is defined by the relation between a corotational rate of Cauchy stress tensor and rate of elastic deformation.Which can be generalized as: where Ĉe is the elasticity tensor of the lattice (in case of orthotropic elasticity) rotated to the current configuration.The resolved stress on the slip system  is calculated by the projection where P() introduced as the symmetric Schmid tensor which geometrically projects the applied stress on the slip system with the slip direction ŝ() and slip plane normal m() .Principle of maximum dissipation has been previously used to derive the flow rule of elastic-plastic materials as well as crystal plasticity (see for instance [21][22][23]).For the derivation below the strict dissipation approach (as in [21]) is followed for simplicity where no defect energy associated with plastic deformation is considered.On the other hand the equivalence between the strictly dissipative and energetic hardening potentials has been demonstrated in [24].
Dealing with large deformation crystal plasticity necessitates choosing the correct stress measure that is power conjugate with the plastic deformation.In [5,23], for instance, the Mandel stress has been used in the intermediate configuration as the power The constrained optimization problem defined above can be reformulated in its dual form using the Lagrangian as where  () are the Lagrange multipliers associated with every slip system.Optimality of the Lagrangian requires its derivative with respect to  to vanish resulting in with the Karush-Kuhn-Tucker (KKT) conditions Eqs. ( 12) and ( 3) together with the last KKT condition assure that For a slip system () to be considered active either the condition  () = 0 or equivalently  () > 0 must be satisfied.

Hardening
In the following, the main mechanism for work hardening is considered to be the impediment of dislocation motion by increase in the forest dislocation density.For each slip system, a Taylor type hardening law [2] is employed where a physically motivated interaction matrix is considered between dislocation density and slip resistance Here  0 is a constant term representing lattice friction in the absence of immobile dislocations,  is the shear modulus,  is the Burgers vector length,  () is the total dislocation density of slip system  and  () is the interaction matrix.
The coefficients of matrix   are usually determined by discrete dislocation dynamics simulations taking into account the geometric relationship between the slip systems.Each coefficient characterizes the strengthening of slip system  due to increase of dislocation density on .Six independent interactions of type: self, coplanar, collinear, orthogonal, glissile and sessile, are sufficient to determine all coefficients [25][26][27].In [27] a straightforward procedure for identifying interaction types between slip systems is given, which is summarized in Appendix D. Table 1 shows the interaction strengths determined by Kubin et al. [26] via discrete dislocation dynamics simulations for an fcc structure.
The evolution of the dislocation density on a slip system  () is governed by the slip rate of that slip system only.As in [28,29] here a conventional approach is followed where the evolution of  () occurs according to the linear ordinary differential equation The terms  ∞ and  ∞ are empirical constants that control the values for saturation of statistically stored dislocation density and the rate of saturation [29].
During the iterative stress update procedure a linearization of the yield function of slip systems with respect to plastic multipliers on each slip systems is required.For the above-mentioned flow stress definition this is obtained as follows: E.S. Perdahcıoğlu

Stress update algorithms
Having defined the flow and hardening rules the next step is to formulate a time integration scheme where the stress at the end of a time step is found which satisfies the flow rule.For this purpose a backward-Euler approach is followed and the incremental problem is solved using an optimization algorithm.
It is important to note that the presented algorithms conserve the rate form of equations presented in Section 2.1 and an approximate time integration to find the elastic and plastic deformation gradients are not required.Only an approximation to the incremental rotation is carried out (as described in Appendix B) which results in an orthogonal rotation matrix.This removes the necessity of adding constraints to ensure volume preserving plastic deformation.
Two rate-independent algorithms are implemented as described below.

Closest point projection (CPP)
Optimality of the Lagrangian given in Eq. ( 11) requires it to be stationary and consequently its partial derivatives with respect to  and  to vanish.Rewriting Eq. ( 12) in terms of the total deformation rate and the stress also considering that it is not satisfied at the start of the increment gives: or alternatively: where  ≡   and  ≡ ̇  and Eq. ( 14) has been implicitly taken into account.
Together with the residual vector on the yield functions for all active slip systems,  ∈ A: the result is a system of nonlinear equations which can be solved using the Newton-Raphson algorithm to determine the plastic multipliers as well as the stress increment together.Using Eq. ( 19) instead of ( 18) is equivalent but advantageous for the scaling of the residual vector for in terms of stress magnitudes.
For efficient implementation of the solution process all tensorial variables and operations are converted to matrix-vector form using Mandel notation and from here on in the text denoted by an underline.The details of the vectorization of variables are discussed in Appendix C.
Linearizing the system of equations, the iterative variables are related to the residuals as: where The first row of equations in (21) can be used to solve for  as a function of  as: which can be substituted in the second row leading to the condensed final set of equations to be solved: During iterations the iterative values of Lagrange multipliers, , might lead to the condition that: for one or more slip systems which results in a non-physical value for the plastic multiplier.Therefore the following correction is applied to the iterative solution which scales the update such that all  are positive: Here  is a multiplier that will scale the iterative solution vector.It is taken as the minimum value of − () ∕ () among the set of slip systems for which the update resulted in a negative .When the solution vector is scaled it can be seen that for that slip system the iterative solution will be exactly the negative of the incremental slip and therefore the updated slip will vanish while all other incremental slip amounts will remain positive.
When converged the system of equations given in (21) have residuals that are equal to zero vectors.This fact and Eq. ( 18) can be exploited to find the algorithmic relation between the incremental deformation rate and the stress: It is important to note that to be consistent with the backward-Euler approach the P tensor at every iteration is updated using the solved plastic multipliers.This requires computing after every solution the incremental rotation that is given in Appendix B.
The difficulty with the above algorithm is to determine the active set of slip systems, i.e. A.Here we refer to [4] which gives a robust algorithm to initialize and update the set of active slip systems during the solution algorithm, which is summarized as follows.An initial trial set of active slip systems is chosen, which can be either the stored active set in the state variables that converged in the previous increment or that consisting of all systems that violate their yield conditions for the elastic trial stress.Once the solution converges two conditions may arise: or, In the case ( 28) is true then the most violating slip system is added to the active set and in the case of (29) the slip system that minimized condition ( 26) is subtracted from the active set.Following these, all variables are reset to their initial values and iterations are repeated.This algorithm however can take a few iterations to determine the active set that will satisfy both conditions.
Furthermore, in the case of linearly dependent slip systems which may occur when the resolved shear stress on two or more slip systems are identical and hardening is isotropic the Jacobian matrix used in the iterative solution given in Eq. ( 24) becomes singular.In this case either SVD or the perturbation method must be used to resolve the indeterminacy.For more information we refer to [4] where different algorithms are investigated in detail for accuracy and efficiency of the computations.
Additionally and in more relevance to the physics of plasticity, as pointed out by several authors [5,10], the selection of the active set has a direct influence on the solution.It might be possible to obtain multiple converged solutions satisfying all the above mentioned criteria, ( 28) and (29).

Interior point method (IP)
In the case of solving for the maximization of plastic dissipation using the interior point method the Lagrangian given in Eq. ( 11) must be rewritten.One of the main differences of the interior point method compared to closest point projection is that there is no need for determining which slip systems are active.Another important difference is that iterations for updating the stress take place in the elastic domain.These two properties make it more likely to reach a global maximum of the plastic dissipation.
Different formulations of the dual form of the optimization problem for the interior point method can be found [12,30].One of the simplest of these can be written as where  (will be referred to as IP parameter) is chosen as a very small number.There are a few important points that can be observed related to this equation.First being that the constraints, i.e. yield functions, are now not exactly satisfied since any  cannot be exactly zero at the optimum.And consequently the KKT conditions are now different than those given in Eq. ( 13): This can be reformulated equivalently using slack variables, , as: with The slack variables are only introduced for algorithmic convenience with the purpose of turning the inequality constraints on the yield functions into equality constraints, thereby also recovering the Lagrange multipliers in the formulation directly.Note therefore that the formulations (32), (33) are equivalent to (30), (31).
The last KKT condition in Eq. ( 33) differs from the consistency equation that needs to be satisfied for consistent plastic deformation.Due to the fact that the product of the yield function value and the Lagrange multiplier does not vanish, choosing the plastic multiplier as equal to the Lagrange multiplier as given in Eq. ( 14) will have two consequences.The first complication will arise from the fact that even when the value of the yield function is a large negative number plastic deformation, although very small in magnitude, will take place.This is a consequence of  not vanishing for non-zero .When  is chosen as a very small E.S. Perdahcıoğlu number, e.g. ≈ 10 −8 , the amount of erroneous slip will be at least a few orders of magnitude smaller than the amount of physical slip.The second complication is that when deformation stops but incrementation in the Finite Element software (or similar) goes on, there will be creep-like plastic deformation.This is a consequence some slip systems being sufficiently close to the yield point to be considered active.It is important to note that however although theoretically these violate the physics of plastic flow when the  is chosen very small the relevance of these complications reduce to that of computational accuracy.
As in the previous section, the stress and the plastic multipliers that maximize the dissipation are found iteratively by searching for the values that ensure that the Lagrangian is stationary: The residual vector on the yield functions for all slip systems becomes: L  ()  =   () =  () +  () (35) and additionally for the slack variables: Once again the result is a system of nonlinear equations the unknowns being the plastic multipliers, the stress increment and the vector of slack variables.
Linearizing the system of equations, the iterative variables are related to the residuals as: where w and Δ indicate matrices with s and s at the diagonal.The last row of equations in (37) can be used to derive the relation between  and : which when substituted in the second row yields: Further substitution of Eq. ( 23) in the second row the final set of equations to be solved becomes: The same scaling update given in Eq. ( 26) is used to ensure that all incremental plastic multipliers and slack variables are positive.Similar to Eq. ( 27) the algorithmic tangent for the interior point method can be derived:

Link to penalty and augmented Lagrangian methods
Treating the stress update as a constrained optimization problem has led the way in development of various algorithms whose goal is to approximate the solution to the exact primal formulation.Among these in the crystal plasticity literature most commonly the penalty approach is encountered which can be seen as an alternate definition for the viscoplastic flow rule.In this approach, generally, the dual form of the optimization problem is formulated as where φ() is an alternative but equivalent form of the yield functions for each slip system,  and  are the penalty parameters.Various choices for these parameters are possible, for instance the choice of  = 1 and  ≫ 1 recovers the common viscoplastic formulations whereas  ≥ 1 and  = 1 can be seen as the linear penalty formulation.
Considering that at the optimum the gradient of the Lagrangian vanishes, the plastic multiplier for each slip system is recovered which can be generalized into the common viscoplastic form [4] Eq. ( 24) for CPP or Eq. ( 40) for IP  =   − Ĉe P   + =   ,  =    Eq. ( 26) Eq. ( 26) IP only end 3) Eq. ( 19) for CPP and Eq. ( 34) for IP   , Eq. ( 20) for CPP and Eq. ( 35) for IP  w , Eq. ( 36 Similar forms of this equation can be found in literature in terms of different flow or creep laws (e.g.[7]).Penalty approach enforces the constraints by penalizing their violation as opposed to the interior point method where every iteration is enforced to stay in the feasible domain.Another difference is the utilization of a logarithmic barrier which enables stricter enforcement of the constraints with stable convergence behavior .The augmented Lagrangian approach can be considered as an extension on the penalty method and being so it exploits its positive properties such as the lack of the need for the active set and regularization against Taylor ambiguity [5].Furthermore using an iterative procedure by shifting the penalized function towards the feasible domain at each iteration it increases the accuracy of the method for a certain choice of penalty parameters.In the augmented Lagrangian approach the iterative update of the plastic slip can be written as where the rightmost term can be substituted with any type of viscoplastic flow rule.Generally for the same level of accuracy a much smaller penalty factor can be used thereby stabilizing the method against numerical problems.Although it is clear that the augmented Lagrangian method is a significant improvement on the penalty method the underlying flow rule of both methods are very similar.

Algorithm
The algorithm is summarized schematically in Table 2.It can be implemented as a standalone routine that returns the updated stress for a given deformation increments as well as in a user defined subroutine of any permitting (commercial) FE software.The incremental spin tensor must be supplied in order to correctly update the orientation of the slip systems with respect to the current configuration.On the other hand if it is not available it can be supplied as a zero tensor which implies that the material coordinates are rotation free.This however will still result in a rotation of the slip systems with respect to material configuration.Therefore, even for small strain analysis, the algorithm will result in texture evolution.
Furthermore, algorithmic tangents for both stress update algorithms, CPP and IP are formulated.These result in a stable convergence behavior in FE analysis.It might occur however that while using CPP, during global equilibrium iterations successive trial displacements might trigger different slip systems.Local stress update iterations in this case might converge and return the algorithmic tangent for the current active set.If however the active set during global (FE) iterations changes the behavior might become non-differentiable and therefore convergence issues might arise.This is an unavoidable consequence of using the active set approach and to the knowledge of the author no solution exists at the stress update level to mitigate this behavior.On the other hand a low-level information exchange between the FE solver and the stress update algorithm, in the case of for instance an in-house FE code, might be exploited to stabilize convergence despite changes to the active set.

Example problems
The algorithms formulated above have been implemented as a standalone FORTRAN module that can be interfaced with Abaqus UMAT, user subroutine.The aim of these examples is to demonstrate the stability and feasibility of the developed algorithms to be used as material models in a full field crystal plasticity Finite Element simulation.

Linear dependence of slip systems
In order to address the relevant problems with rate independent crystal plasticity implementation the simple shear deformation of a single crystal proves to be a good demonstrator.A non-hardening fcc type crystal lattice is chosen which is deformed under shear.For this example the crystal rotations in the algorithm are turned off and a shear strain of 0.02 is prescribed in 1000 increments, while the spin tensor  is also set to .The lattice axes are not parallel to the global axes and are rotated with Euler angles being (36 • /18 • /90 • ) according to Bunge convention (see Appendix D for details).Material properties used for this analysis are shown in Table 3.
In the CPP algorithm mentioned above this results in a singular Jacobian matrix when more than 5 slip systems become active (Taylor ambiguity) due to the linear dependence of slip systems.Therefore the amount of slip is distributed as a result of the generalized inverse approach mentioned above and described in [4].The IP algorithm however is impervious to this condition and is able to reach a solution as shown below.For validation, the input data is chosen as given in [4] and the activation as well as the progression of slip activity on different slip systems are shown in Fig. 1.
The results of the CPP algorithm are identical to that given in [4] however a large discrepancy between the IP algorithm can be observed after a shear strain of 0.6 × 10 −2 when the number of active slip systems suddenly jump to 6.The discrepancy continues until the end of the deformation.In order to give a better understanding of this situation the individual stress components and normalized plastic dissipation that is computed at every time increment are shown in Figs. 2 and 3, respectively.It can be seen clearly that the stress response as well as the dissipation obtained with two different algorithms, which have different slip systems active, are identical throughout the deformation.This can be understood by considering that both algorithms aim to find the stress and the plastic multipliers that maximize the plastic dissipation.Due to the linear dependence of the plastic multipliers different combinations lead to the same dissipation.

Crystal rotation
Another important aspect of the implementation is the correct update of the crystal rotation.During simple shear deformation the material axes will rotate due to the associated material spin whereas the crystal lattice will rotate only due to the non-inelastic spin causing different slip systems will be activated and inactivated during the deformation.The results of both algorithms are shown in Fig. 4 and can be compared against [5].This problem is also a good demonstrator of the consequences of choosing a hypoelastic material behavior (corotational) instead of a hyperelastic one.
The material is chosen to be the same non-hardening material with an fcc crystal structure, the parameters used for the simulation are given in Table 4.For this example the algorithm is driven with a prescribed total deformation gradient with  12 = 4.0 in different number of increments.The backward Euler approach ensures that the rotation of the crystal is enforced to be the one at the end of the increment.As shown in the figure, the algorithm stably converges with an increasing number of increments and that the crystal rotation is captured accurately.In Fig. 4 the rotation of the crystal is shown using a stereographic projection of the crystal axes [1,0,0],[0,1,0] and [0,0,1] (using [31]) throughout the simulation which show a good agreement with [5] and [13].The agreement  of the results with the hyperelastic implementation in [5] can be explained by the fact that within the stress limits obtained in this material the non-linear elastic regime is not reached and with the corotational formulation the stress rotations are properly captured.
In the current implementation the update of slip system orientations are done every iteration and do satisfy the system of equations at the end of the increment, in line with the Backward Euler approach.As can be seen however in the algorithmic tangents developed for the algorithms the effect of the change of slip systems on the linearization is not taken into account explicitly.This causes the convergence rate to be less than quadratic.Another possibility is to use a staggered rotation approach where the orientation of the slip systems are kept constant during an increment and only updated at the end of the increment.This will result in an increase in the global convergence rate at the FE level.In order to quantify the inaccuracy this method would create the same  simulation has been carried out with staggered rotation for 1000 increments and compared with the case where rotation is updated implicitly.The results are shown in Fig. 5.
It can be seen that the difference in the results are very small and as expected the staggered approach lags slightly behind in capturing the effect of crystal rotation.The difference however depends on the increment size and this approach should be attempted carefully.On the other hand the example shows a total shear of 4 that is reached in 1000 increments which can be considered as a large increment size.It corresponds to a shear strain of 0.4 in 100 increments, which is more of a common range of strain in actual simulations.

IP parameter
Iterative solutions of the interior point algorithm do not violate the constraints but also the active constraints are not exactly satisfied.The  parameter introduced in Section 3.2 determines how close the solutions can get to the constraint; the smaller the  the smaller is the value of the constraint.As can be expected this has an effect on the computed plastic multipliers since they are equal to the Lagrange multipliers in the algorithm.A larger  results in a smoother transition from elastic to plastic domains for every slip system.This is advantageous in terms of a FE solution algorithm due to the fact that the transition point becomes differentiable leading to more stable convergence.On the other hand when the parameter is large it can be expected that the transition occurs  earlier and accordingly at a certain stress level the plastic strain becomes larger.When the material shows hardening this might also lead to additional hardening due to an increase in the amount of plastic deformation.
In order to show these effects the set of material parameters are chosen to reflect a hardening fcc type Aluminum-like material, as shown in Table 5.The coefficients for the hardening interaction between slip systems are those given in Table 1.The parameter  is varied between 10 −3 and 10 −7 and a uniaxial tensile deformation is prescribed in 200 increments where the strain  11 is prescribed while all other strain components are solved using the condition that the stress components except  11 vanish.The results are plotted in Fig. 6.
It can be seen that for  values smaller than 10 −3 the variation in the results are negligible.For  = 10 −3 the yielding starts earlier and therefore for the same stress level there is more plastic strain which results in more hardening in later stages in accordance with the expected behavior.

Poly-crystal FE simulation
The above algorithm is implemented in commercial FE software Abaqus as a UMAT user defined material subroutine.Abaqus provides the previous stress tensor rotated to the current configuration and also the incremental deformation rate tensor.It however does not provide the incremental spin tensor and therefore for being able to update the crystal orientations it must be computed using the provided deformation gradient tensor at the start and the end of increment.The details of this computation is given in Appendix E. The geometries used as volume elements in an FE simulation are generated using the open source Voronoi tessellation software Voro++ [32].It has been slightly modified to output the values in double precision which is needed for compatibility with other software and also with the capability of performing Lloyd iterations for a more regularized grain structure similar to the implementation in [33].Once a periodic tessellation is obtained the hull is found considering the periodicity of the cell surfaces.
The generated Voronoi cells are then input to the open source mesh generation software GMSH [34].The previously identified periodic hull is required as input in order to obtain a periodic mesh.The GMSH version available at the time of utilization is only capable of creating a tetrahedral mesh with a periodicity constraint.The generated mesh data is then used to create an input file to Abaqus where periodic boundary conditions are applied as multi point constraints between the corresponding nodes at each periodic surface.The simulations are carried out using Abaqus software version 2022.
The simulation consists of a tessellation of 100 grains with random orientations generated using the scipy python library (function scipy.stats.special_ortho_group.rvs()).The mesh is made of 127346 C3D4 (linear tetrahedral) elements.
A uniaxial deformation is prescribed in the  direction of a total of  11 = 0.2 engineering strain with all remaining surfaces being traction free.The mesh and the deformed structure are shown in Fig. 7 and the obtained homogenized engineering stress-engineering strain response is shown in Fig. 8.The simulation is repeated three times, with  = 10 −4 in 1000 increments, with  = 10 −5 in 2000 increments and with  = 10 −6 in 3000 increments.The number of increments are chosen to result in the fastest computation which resulted in wall-clock duration of 72, 142 and 245 min, respectively, on an Intel Xeon workstation using 8 CPUs running at 2.90 GHz.
Fig. 6 shows that there is a slight difference in the overall stress-strain behavior caused by the large  value that results in numerical slip activity as discussed in Section 4.3.On the other hand there is a significant gain in using a higher  due to the fact that convergence behavior is positively affected by a smooth transition to plasticity for slip systems.

Simulation of neck formation
In order to demonstrate the stability and efficiency of the interior point algorithm in a large deformation finite element setting the necking of a strip under plane strain tension is studied.A similar study is carried out in [5] to show the efficiency of the augmented Lagrangian algorithm.The strip is 8 mm in width and 60 mm in height with an imperfection imposed in the center as a 90% reduction in width to promote localization.It is discretized using 400 quadratic quadrilateral elements with reduced integration (CPE8R in Abaqus software).The material is taken as a single crystal with an orientation (57 • /28 • /0 • ) and parameters given in Table 6 which are tuned to promote early localization.A total elongation of 5 mm is prescribed at the top edge and the simulation required 1267 increments to conclude.The distribution of the total plastic slip as elongation proceeds is shown in Fig. 9 showing the extent of localization in terms of plastic slip as well as the displacement of the mesh in the middle section.The force-displacement response of the strip, given in Fig. 10, shows that localization starts as early as about 1.45 mm in elongation.

Conclusion
A novel large deformation algorithm for rate-independent crystal plasticity is presented.The starting point being the maximization of plastic dissipation, the interior point method is utilized in deriving the flow rule and the solution procedure.
The algorithm is impervious to the ambiguity that occurs when linear dependence of slip systems is encountered and provides a smooth transition to yielding that result in a robust and efficient determination of slip rates together with the stress.
The algorithm is developed in rate form and in the current configuration which makes it very suitable for large deformation implementation in finite element software.The handling of large rotations and the efficiency of the computation are demonstrated with numerical examples.

Declaration of competing interest
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.

Fig. 5 .
Fig. 5. (left) Shear stress vs shear strain in a single crystal with orientation (0 • /0 • /0 • ) under prescribed deformation  12 = 4.0 computed using the IP algorithm with the implicit and staggered rotation approaches in 1000 increments.(right) Trace of the positive ⟨100⟩ as well as [111] and [ 11 1] directions computed with the staggered approach shown on the stereographic projection.

Fig. 6 .
Fig. 6.Effect of the IP algorithm parameter, , on the overall stress vs strain behavior during uniaxial tension of an fcc single crystal with the orientation (5 • /11 • /17 • ).

Fig. 7 .
Fig. 7. Periodic volume element comprising 100 randomly oriented fcc crystals.(left) Original mesh with colors showing different material definitions.(right) Deformed mesh after  11 = 0.2 with the contour showing the total slip amount.

Fig. 8 .
Fig. 8. Engineering stress vs engineering strain response of the poly-crystal for three different  values.

Fig. 9 .
Fig. 9. 2D plane strain tension simulation of a single fcc crystal oriented in (57 • /28 • /0 • ).Figures show the progression of the deformed mesh superposed with the amount of total slip for elongations of 1 to 5 mm, left to right, respectively.

Table 1
[26]ficients for different interaction types for face centered cubic lattice[26]. p .Equivalently one can use power conjugates Cauchy stress and Lp at the current configuration.This leads to the formulation of maximum dissipation in the current configuration

Table 4
Crystal plasticity parameters used in Section 4.2.

Table 6
Crystal plasticity parameters used in Section 4.5.