The development of the numerical integration scheme consists of the following:(i) identification of primary variables in constitutive formulation, (ii) discretization and numerical integration of equations in time domain, (iii) development of Newton–Raphson iterative scheme, and (iv) development of time sub-stepping algorithm to increase or decrease time step size, depending on the incremental values of primary variables.
3.2. Discretization of Constitutive Equations
A time integration scheme is executed in sample coordinate axes through discretizing the deformation history in the time domain and subsequently numerically integrating constitutive equations for each time step. In order to define a deformation time history, the configurations of a material point are considered at time and . In this, a relation of is used, where and represent time at the start and end of the time step, respectively. Afterwards, the magnitudes of all variables are evaluated at and , and denoted with the subscripts n and , respectively. The numerical time integration scheme is based on the following assumptions:
The total deformation gradient and velocity gradient are given.
The values of variables , , , , and at time are known.
The initial values of time-independent slip and twin systems’ vectors, elasticity tensors and , initial crystallographic orientation of crystal , initial kinetic flow rule (m and ), and hardening parameters (, , , ) are used as an input.
The output of the numerical integration scheme provides updated values of variables as:
,
,
,
, and
at time
. The constitutive equations are discretized through fully implicit time integration procedure using backward Euler scheme. According to the kinematic formulation, an incremental form of Green strain tensor in terms of the symmetric part of the velocity gradient, Equation (
25) in [
53], can be written as:
Integration of Equation (
40) using backward Euler scheme provides an updated value of Green strain tensor, as follows:
where
and
are the symmetric parts of velocity gradients
and
, respectively. The incremental values of these parameters can be illustrated as:
where an updated value of right Cauchy–Green deformation tensor in third intermediate configuration can be defined as:
Furthermore, an incremental magnitude of spin tensor
is represented as:
In addition,
is an updated symmetric part of
. Here,
is an incremental form of right Cauchy–Green deformation tensor in second intermediate configuration, which can be expressed as:
The incremental value of the plastic velocity gradient in second intermediate configuration can be defined as:
where
represents an incremental change in the volume fraction of
twin system,
is an updated value of the shear strain rate of
-slip system,
is an updated value of shear strain induced by
i-twin system,
and
are the updated magnitudes of Schmid and twin orientation tensors in second intermediate configuration, and
is an incremental change in the value of the plastic deformation gradient. Substitution of
in Equation (
43) gives the following:
Considering the effects of rigid body rotation tensor
on Schmid and twin orientation tensors, Equation (
48) can also be written as:
The forward mapping approach for Schmid and twin orientation tensors can be implemented into Equation (
49) as:
where
is the initial rotation matrix that is used to transform crystal coordinates to sample coordinate systems through Euler angles, and
and
are the initial vectors representing
-slip system in reference configuration. It was stated previously that the rotation matrix (Euler angles) can be updated through a rigid body rotation tensor as follows:
. The rotation matrix can also be used to transform Schmid orientation vectors from reference to first intermediate configuration as follows:
and
. In this condition, Equation (
50) can be rewritten as:
In Equation (
51),
and
are the slip incremental values of direction and area normal vectors of
-slip system in third intermediate configuration. It is clear from Equation (
51) that the updated rotation matrix
is used to transform Schmid orientation vectors from reference to third intermediate configuration. Similarly, twin orientation vectors are given as:
Substitution of Equations (
51) and (
52) in (
49) provides
Here,
and
are the updated Schmid and twin orientation tensors in third intermediate configuration. The updated values of shear strain for
-slip system and volume fraction of twinned martensite at time
can be defined by Equations (
54) and (
55), respectively, as:
As discussed previously, the rigid body rotation tensor
updates the crystal orientation (Euler angles) matrix
, which can be used to transform the fourth order elasticity tensors for slip
and twinned
regions to the third intermediate configuration as follows:
Any one of the elasticity tensors can be obtained in terms of another by using coordinate transformation rule defined as:
Here,
is the transformation matrix that relates lattice orientations in twinned and untwinned (slipped) regions. The components of transformation matrix
can be defined through Equation (
59), given by [
26], as:
where
is the area normal vector of the twin plane, and
is the Kroneker delta. An equivalent elasticity tensor can be calculated as:
Furthermore, an updated form of equivalent second Piola–Kirchhoff stress tensor can be estimated as:
An evolution equation for rigid body rotation tensor
is integrated through the exponential map discussed by [
58] as follows:
where an updated value of the spin of lattice
can be estimated through the skew-symmetric component
of the velocity gradient
as:
By using backward mapping, the skew-symmetric component of velocity gradient
can be estimated as:
Here,
is the updated skew-symmetric component of
. Substitution of
from Equation (
63) in (
64) provides
For small elastic strain problems, the value of the right Cauchy–Green deformation tensor
is typically small. In this case, Equation (
65) can be written as:
An incremental value of the slip resistance of
-slip system can be evaluated using a backward Euler scheme as follows:
In the present model, slip resistances for all
-slip systems are considered to be similar. Therefore, Equation (
67) can be modified as follows:
In Equation (
68),
is an incremental saturation value of slip resistance, which can be calculated as:
where
is the initial value of saturation slip resistance,
is an incremental initial value of the slip system shear strain at the initial value of saturation slip resistance,
a is the material parameter,
is the material slip-hardening parameter,
(
) represents the number of non-coplanar twin systems to slip system, and
is a material parameter. In addition, the twin resistance can be estimated using Equation (
34) as:
The resistance of all twin systems are assumed to be identical. Therefore, Equation (
70) can be expressed as:
The updated driving potentials for slip and twin mode of deformations can be estimated using Equations (
29) and (
30), respectively as:
The parameters
and
are the thermal analogues for
-slip and
i-twin systems, respectively. For an isothermal process, these factors will remain constant during the deformation process.
is a dimensionless dislocation parameter, and it is assumed to be constant throughout the deformation. An equivalent modulus of rigidity
can be estimated as:
The incremental forms of crystal defect microstrain parameters for the slip and twin modes of deformations can be calculated using Equations (
35) and (
36), respectively, as:
The updated stress-like functions
and
for slip and twin modes of deformation can be calculated using Equations (
37) and (
38), respectively, as:
Furthermore, the updated Cauchy stress tensor
can be calculated using second Piola–Kirchhoff stress tensor,
, as:
3.3. Newton–Raphson Iterative Scheme
In the preceding section, a set of couple nonlinear algebraic Equations (
61), (
62), (
68), (
71) and (
55) for the primary variables
,
,
,
, and
, respectively, were developed. In this, the updated form of these primary variables are calculated. A set of primary variables can be represented by a vector
, where
The primary variables are the main constituents (directly or indirectly) of the constitutive model. An elastic-plastic response of a system mainly governs by these variables. In order to obtain an overall response of a material, a set of five equations in terms of primary variables are constructed in a residual format. A residue of the system of equations can also be represented in a vector form as
. The components of
are given through Equations (
81)–(
85) as:
In the next step, Equations (
81) and (
82) are solved using a full Newton–Raphson (N-R) method, since these two are implicit in nature; however, the rest are explicit. In the current work, a two-level iterative scheme, similar as that presented by [
56], is used to obtain the values of primary variables. In the first level of iteration, the N-R method is used to solve Equation (
81) for
by assuming the best possible values of the other primary variables
. Once the updated value of the second Piola–Kirchhoff stress tensor is obtained, the second level of the iterative procedure is performed, which includes an N-R solution of the slip resistance
from Equation (
82). This considers an updated value of
and the estimated values of
,
, and
. Finally, updated values of the twinned martensite volume fraction
, rigid body rotation tensor
, and twin resistance
are calculated from Equations (
83)–(
85), respectively.
Convergence Criterion
The convergence criterion is required to terminate the iterative loop once the solution is assumed to be sufficiently accurate. The convergence criterion for the presented two-level iterative procedure is based on the variation of
-norm for
, and
. If the
-norm of the residuals is less than an imposed tolerance, then the incremental solution of the updated primary variables is converged (fully elastic). Otherwise, the trial value must be updated iteratively (based on the calculated value) until the residual satisfies the convergence criterion, given as:
In the current Newton–Raphson iterative scheme, iterations are carried out unless and until the variation in the
-norm of the residuals for
and
satisfy the conditions in Equations (
87) and (
88), respectively, as follows:
In the current work, a time sub-stepping algorithm (TSSA) is used in the numerical integration scheme to control the time step size. The advantage of using time sub-stepping in an iterative procedure is to improve the convergence of a solution by reducing the time increment once needed. The TSSA must also be capable of increasing the time step size at a material point where convergence can easily be achieved to reduce computational time. The TSSA can reduce and increase the time step size based on the incremental variation. In this algorithm, if any of the convergence conditions, Equations (
87) and (
88), are not satisfied, then the N-R iterative scheme will call TSSA, which controls the time step size according to the incremental variation in the values of primary variables.
3.4. Time Sub-Stepping Algorithm
The proposed TSSA is based on the ratio of the maximum and desired incremental values of the primary variable. The ratio is represented by the parameter
as:
where
is the maximum value of
over all the crystals, all the integration points are in finite element mesh, and
is a desired incremental value of
in the numerical algorithm. The value of
is normally regulated by the computational performance and accuracy of the numerical integration procedure. In a fully implicit numerical scheme, the computational performance does not usually becomes a challenge, but the accuracy does. Therefore, the value of
is mainly controlled by the accuracy of the numerical solution. In accordance with the two-level iterative scheme,
is defined as a set of three primary variables
, where
In the proposed time sub-stepping algorithm, three ranges of values are defined for the parameter
: (i) If the parameter
exceeds 1.25, the estimated value of
is rejected, and the new time increment (25% smaller than the previous) will be defined. This condition makes sure that the difference between the maximum and desired value of
will not reach beyond 25%, which could produce an inaccurate solution. (ii) If
lies in the range of 0.8 to 1.25, then the estimated solution is accepted, defining the new time step as
. In this condition, the new time step size
is more or less identical to the previous
. (iii) If
is less than 0.8, then the estimated solution is assumed to be converged, and a new time step size is defined that is 25% larger than the previous. This condition ensures that the solution is well converged, so the time increment could be increased to reduce computational time. A summary of the time sub-stepping algorithm is given in
Table 1. The numerical integration scheme for elastic-plastic deformation of a crystal based on crystal plasticity formulations is summarized in
Figure 2.