Finite Strain Constitutive Modelling of Shape Memory Alloys Considering Partial Phase Transformation with Transformation-Induced Plasticity

This paper presents a unified modelling effort to describe partial phase transformation during cyclic thermo-mechanical loading in Shape Memory Alloys (SMA). To this purpose, a three-dimensional (3D) finite strain constitutive model considering TRansformation-Induced Plasticity (TRIP) is combined with a modified hardening function to enable the accurate and efficient prediction of partial transformations during cyclic thermo-mechanical loading. The capabilities of the proposed model are demonstrated by predicting the behavior of the material under pseudoelastic and actuation operation using finite element analysis. Numerical results of the modified model are presented and compared with the original model without considering the partial transformation feature as well as with uniaxial actuation experimental data. Various aspects of cyclic material behavior under partial transformation are analyzed and discussed for different SMA systems.


Introduction
Shape memory alloys (SMAs) are smart alloys possessing the unique behaviors known as Pseudo-Elasticity (PE), one-way and two-way Shape-Memory Effect (SME). After being subjected to mechanical deformation while at the high temperature austenitic phase, the alloy is able to recover the original shape through a solid-state, diffusionless phase transformation caused by the imposition of a stress (i.e., PE) and/or temperature (i.e., SME) field [1]. The physical mechanism behind SMA properties is attributed to such solid-state, diffusionless phase transformation between a high-temperature, high-symmetry austenitic phase and a low-temperature, low-symmetry martensitic phase, that takes place when the alloy is subjected to certain thermo-mechanical loading conditions [2].
Such unique behavioral characteristics make SMAs largely employed in a wide range of innovative applications [3]. Particularly, thanks to PE, good biocompatibility, and high corrosion resistance, SMAs are extensively used in the biomedical field, as stents, occluders, or surgical tools, while thanks to the SME and high output energy density, coupled with several advantages for system miniaturization (e.g., high power-to-mass ratio, maintainability, reliability, clean and silent actuation), SMAs are largely investigated as actuators in various engineering industries including, e.g., aerospace [4], civil [5], automotive [6], and recently wind energy [7][8][9][10].
When SMAs are subjected to cyclic loading, either thermal or mechanical, they most often present two important characteristics, i.e. TRansformation Induced Plasticity (TRIP) strains and partial transformation.
The first characteristic is generally exhibited during the so-called training procedure, which is typically employed to stabilize the cyclic response of SMAs before their use [11][12][13]. During this procedure, permanent microstructural changes are often introduced in the alloy, which leads to the generation of an internal stress field and large irrecoverable TRIP strains. Consequently, cyclic loading tremendously affects SMA characteristics and most importantly the performance (e.g., actuation work) and fatigue life of the SMA component. As an example, Phillips et al. [13] showed that the internal damage in a SMA actuator subjected to actuation fatigue evolves in a nonlinear manner, with a rapid damage nucleation during the training period associated to the irrecoverable strain accumulation. Moreover, the training process can be used to induce the so-called two-way SME in actuators [14], that is the capability of a SMA to change its shape reversibly during subsequent thermal cycles between its transformation temperatures without applying any external mechanical load.
The second characteristic is exhibited when SMAs are subjected to a mechanical or thermal load that is reversed before the completion of phase transformation. Reversal of transformation direction while the alloy is in a mixed phase state leads to the formation of a Partial Transformation (PT) branch. In this case, SMAs are capable of ''remembering'' their previous states associated with the transformation reversal prior to a complete phase transformation [15]. Such a feature is also known as reversal memory. Depending on the applied loading history, different PT situations may arise, including partial loops with one single reversal point, closed loops with two reversal points (see Fig. 1), and partial loops with multiple reversal points. From a microstructural point of view, such phenomenon is associated to an ordered alteration of the material microstructure, dictated by the achievement of a critical threshold energy value by each grain [16]. When a threshold is reached, a certain group of grains is transformed, and the phase of the SMA is changed. In case of a transformation reversal, alternations take place depending on this threshold. Assumed that the same energy level state is achieved after a series of loadings, the material also obtains the respective microstructure and, consequently, the same phase. Based on these assumptions, the formation of closed PT cycles and memory wipe-out phenomenon can be explained. As a result, PT largely affects the amount of generated transformation strains, the actuation frequency [15,17], and most importantly the fatigue life of the component [18]. In fact, it has been shown in ref. [12] that operation in PT conditions enhances the actuation frequency of antagonistic SMA wire actuators and reduces the developed stresses and activation temperatures, which can be translated to a possible energy efficient operation with prolonged thermo-mechanical fatigue life. Catoor et al. [19] showed that Nitinol stents experience cyclic loading from cardiac cycles and musculoskeletal movements in the mixed-phase state and their lifetime is affected by operation in PT conditions. It is clear from the above discussion that these two features are very likely to be coupled during cyclic loading and their consideration is fundamental for the successful design of SMA-based structures with increased performance and fatigue lives. Accordingly, it becomes essential an associated accurate modeling to support the design.
The second characteristic is generally accounted for by describing the behavior of conventional SMAs without considering irrecoverable strains and damage, which is sufficient for the design of components where the operating temperatures, maximum stress levels, and number of actuation cycles are relatively low. These contributions include differential equations and Duhem-Madelung type models [15,33,34], Preisach-type models [35][36][37][38][39], purely phenomenological models and state equations based on principles of thermodynamics and/or plasticity [20,[40][41][42][43][44].
To the authors' knowledge, the modeling of the complex behaviour of SMAs under the two described characteristics, accounted simultaneously in thermal and/or mechanical cycles, is still missing in the current literature.
Motivated by the described framework, the aim of the present work is to address SMA features under such two important characteristics through a unified modelling effort.
To this purpose, the proposed formulation starts from the three-dimensional phenomenological model recently published in ref. [31], that has been validated against experimental data on NiTi and NiTiHf SMA materials. The choice of dealing with the model by Xu et al. [31] is motivated by the following reasons: (i) the model is formulated in a finite strain framework, which is well suited in case of SMA actuators experiencing 30% or higher TRIP strains during their lifetime or in the presence of cracks [45]; (ii) the model describes the two-way SME exhibited by trained SMAs at load-free conditions; (iii) the model considers stress-dependent TRIP evolution under multiaxial stress state. Despite these advantages, the model by Xu et al. [31] is limited to an accurate and efficient description of only full phase transformation.
Accordingly, an efficient and easy-to-implement modified hardening function [15] is introduced in the present paper to account for PT behaviour under cyclic loading. This constitutive modeling approach has been inspired by the models presented in refs [46] and [41] and combines some of their characteristics. However, it does not require the determination of any additional material parameters rendering the calibration process easier and more straightforward. In fact, the model can be calibrated by using the same data extracted from uniaxial characterization processes often followed in the literature.
In order to evaluate model accuracy and reliability, several numerical simulations are performed, including uniaxial pseudoelastic and actuation tests and the finite element analysis of a tube-shaped actuator where the effect of geometric non-linearity is quite evident due to the large rotations involved. Results from numerical simulations of uniaxial actuation tests are validated through a comparison with few data from a conducted experimental test on a NiTiHf high temperature SMA. Various aspects of SMA behavior including TRIP and PT are analyzed and discussed.
The paper is organized as follows. Section ''Model Formulation'' presents model equations, while Sect. ''Numerical Results'' shows the results of numerical simulations. Section ''Comparison with Experiments for NiTiHf High Temperature SMA'' reports the comparison with experiments. Finally, conclusions and summary are given in Sect. ''Conclusions''.

Model Formulation
The present section focuses on the formulation of the proposed model. The extension proposed herein starts from the thermodynamically-consistent finite strain model formulation recently presented by Xu et al. [47], and its continuous development considering TRIP strain and loadfree two-way SME after cyclic loading [31]. The interested reader may refer to [31,47] for details about the original model formulation and to Appendix A for modeling preliminaries.

Kinematic Assumption
In the framework of macroscopic modeling and of finite strain continuum mechanics, an additive decomposition of the rate of deformation tensor D into an elastic part D e , a transformation part D tr , and a TRIP part D tp , is considered as follows [31]: Combing Eqs (A.12)-(A.15), the total logarithmic strain tensor h (or true strain) can be additively split into an elastic part h e , a transformation part h tr , as well as a TRIP part h tp : Here, h tr is the transformation logarithmic strain tensor accounting for the recoverable inelastic strain associated with phase transformation, while h tp is the TRIP logarithmic strain tensor describing the irrecoverable transformation-induced plastic strain.

Thermodynamic Potential
Based on the work by Xu et al. [31] and within the classical thermodynamic framework for dissipative materials, the Kirchhoff stress tensor s and temperature T are assumed as external state variables, while fn; h tr ; h tp ; bg is assumed as a set of internal state variables. Particularly, n is the martensitic volume fraction ranging between 0 (pure austenitic phase) and 1 (pure martensitic phase) and b is the internal stress tensor defining the internal stress field generated as a consequence of the training process. Accordingly, the Gibbs free energy G depends on such variables as follows: where S, a, c, s 0 , and u 0 denote the effective compliance fourth order tensor, thermal expansion second order tensor, specific heat, specific entropy at the reference state, and specific internal energy at the reference state (defined at temperature T 0 ), respectively, while f is the hardening function. Particularly, S, a, c, s 0 , and u 0 are evaluated in terms of the martensitic volume fraction n by means of rule of mixtures. The specific form adopted for f allows for an accurate and efficient description of PT, as detailed in the next subsection. Finally, it is noted that the Kirchhoff stress tensor s is related to the Cauchy or true stress tensor r by means of the relationship s ¼ Jr. Here, J is the determinant of the deformation gradient F (see Appendix A for details) and is approximately equal to 1 since phase transformation is assumed to be volume preserving [31,47].

Hardening Function for PT Loops
An hardening function is included in the Gibbs free energy (3) proposed by Xu et al. [31] to account for the polycrystalline hardening effects. In order to describe SMA behavior under PT during cyclic thermo-mechanical loading, including reversal point memory and the associated memory wipe-out, the hardening function is here modified according to the approach proposed by Karakalas et al. [15]. Such an approach considers the scaling hypothesis by Ivshin and Pence [46] on the value of the martensitic volume fraction that is used only in the expression of the hardening function and uses interpolation functions for the approximation of the martensitic volume fraction. Therefore, after the material exhibits a transformation reversal, the value of the martensitic volume fraction is stored and a set of linear interpolation functions are defined. The two reversal values are distinguished since each one has different effect on the evolution of the respective hardening expression.
Accordingly, the hardening function is updated introducing a pair of interpolation functions to scale the martensitic volume fraction, as follows: where a 1 , a 2 , and a 3 are material parameters and the four smoothing parameters n 1 , n 2 , n 3 , and n 4 are introduced to effectively treat the smooth transition characteristics at the initiation and completion during phase transformation. The modified value for n in the expression of f for the forward transformation is of the following form: where n FR and n RF are the reversal point passing, respectively, from forward to reverse transformation and from reverse to forward. For the purpose of illustration, the reader may refer to points FR and RF in Fig. 1. Similarly, the modified values for n in the expression of f for the reverse transformation is of the following form:

Constitutive Equations
Following the classic thermodynamic principles and Coleman-Noll procedure, constitutive equations can be derived from the Gibbs free energy (3). Accordingly, the following relationships between the Kirchhoff stress and logarithmic strain and the entropy and temperature are obtained: A reduced form of dissipation inequality can be derived by substituting the above constitutive equations into the Clausius-Planck inequality [31,47]: Transformation Strain Evolution Law The evolution law for transformation strain is described as follows [31]: Shap. Mem. Superelasticity (2021) 7:206-221 209 where K fwd and K rev are, respectively, the forward and reverse transformation direction tensor, defined as: Here, h trÀr and n r represent the value of transformation strain and the martensitic volume fraction at the reversal point, respectively; s eff ¼ s þ b is the effective stress tensor that can be used to achieve the two-way SME of trained SMAs under load-free conditions, s eff 0 and s eff being its deviatoric part and von Mises equivalent effective stress; H cur is an exponential function introduced to calculate the value of current transformation strain given an effective stress state [48]: where H min corresponds to an observable two-way SME strain for pre-trained SMAs or some SMAs experiencing particular production process such as extrusion and aging under stress, s crit denotes a critical stress value below which only H min is observed, k t is a curve-fitting material parameter, and H max is the maximum (or saturated) value of transformation strain, defined as: Here, f d is the accumulation of orientated martensitic volume fraction which is defined in the next section, while and H max f are the values of H max before and after the cyclic loading. The definition of H max takes into account for a degradation phenomenon that sometimes exists for the value of maximum transformation strain as a result of the accumulation of retained martensite [49].

TRIP Strain Evolution Law
The evolution law for the TRIP strain, considering the effect of stress multiaxiality, is as follows [31]: where K tp fwd and K tp rev are the forward and reverse TRIP direction tensor, respectively, defined as: Here, C p 1 and C p 2 are material parameters dictating the magnitude and evolving trend for TRIP during cyclic thermo-mechanical loading, while f d is the accumulation of oriented martensitic volume fraction, defined as follows: in which the oriented martensitic volume fraction is calculated as: Internal Stress Evolution Law The following exponential evolution law is used for the internal stress [31]: r b being a model parameter representing the maximum (or saturated) magnitude of the internal stress and k 1 being a material parameter controlling the evolution trend.

Transformation Function
The following transformation function is adopted [31]: where p is the general thermodynamic driving force conjugated to n: and the critical value Y is constructed to be stress-dependent quantity with a reference critical value Y 0 and an additional parameter D: It is worth noting that Eq. (21) can be obtained by substituting the evolution laws (10) and (14) into the reduced form of dissipation inequality (9). Accordingly: The model is completed by the Kuhn-Tucker conditions, defined as follows:

Advantages of the Proposed Model Formulation
It is worth highlighting that the proposed formulation considering both TRIP and PT description under finite strain has several advantages. From the modeling point of view, the modification of the hardening function in Eq. (4) is simple and allows for an accurate description of PT loops (as demonstrated by the numerical tests in the next sections). Moreover, it is worth highlighting that the proposed hardening function ensures a thermodynamic consistent model formulation. This means that, for a closed partial transformation cycle, the Gibbs free energy should also return to its initial value after a loading cycle (isothermal, iso-stress, or combined thermomechanical one) to preserve the thermodynamic consistency of the model, as largely discussed, e.g., in ref. [48]. In the absence of TRIP and assuming a small strain framework, our formulation recovers the modeling formulation proposed by Karakalas et al. [12]. Following the analytical definitions reported in Appendix A of this work [12] and by means of numerical integration, it can be verified that the Gibbs free energy returns to its initial value for closed partial transformation cycle in case of both isothermal and isobaric loading cycles. Therefore, the introduction of the proposed hardening function does not affect the thermodynamic consistency of the original model formulation.
From the calibration point of view, no additional model parameters are introduced. Accordingly, model parameters are all physical and can be derived from experimental data by following the procedure described in [31].
From the numerical point of view, the algorithmic treatment of model equations within a finite element framework is efficient and effective for our case. In fact, the elastic predictor/inelastic corrector scheme proposed in [31] is here applied to compute the variables at the current load step based on the variables known from the previous one. Compared to the implementation proposed in [31], the modification of the hardening function implies the introduction of two additional variables, i.e., n FR and n RF (see Sect. ''Hardening Function for PT Loops''). These variables need to be stored, for each integration point of the finite element mesh, at the beginning of the scheme if a transformation reversal is detected. The numerical advantages of the proposed formulation compared to the one proposed in [41], where the martensite fraction values and the direction of the transformation have to be checked before and after each load increment, are thus clear.

Numerical Results
To demonstrate the capabilities of the proposed model to capture SMA behavior under training and PT thermo-mechanical cycles, several boundary value problems considering different loading conditions are investigated. Particularly, we consider various loading conditions typical of pseudoelastic wires for biomedical applications or linear and torsional actuators. For the simulation, the commercial finite element software ABAQUS/Standard is used, after having implemented the algorithm, described in [31] and modified with the proposed hardening function, into a userdefined material subroutine (UMAT).

Uniaxial Pseudoelastic Tests
At first, three pseudoelastic tests on a single eight-node hexahedral element under force control and prescribed homogeneous temperature are simulated. We adopt material parameters reported in Set 1 of Table 1 [31] related to Ni 40 Ti 50 Cu 10 (at.%) SMA. Since the coefficient of thermal expansion for NiTi is calculated to obtain a very low value and because its effect is negligible compared to phase transformation, it was assumed to be equal to zero. The material properties selected for this example were used to verify that the unified model was able to predict the material behavior presented by Xu et al. [31] under complete transformation pseudoelastic cycles. It is remarked that the proposed formulation does not require additional parameters to describe minor loops.
The first pseudoelastic test involves 25 tension cycles, where the applied stress varies between 0 and 550 MPa (see Fig. 2a), and the temperature is set equal to 360 K. During cycling, PT, coupled with TRIP, takes place. Figure 2c and e show the stress-strain and strain-loading step parameter curve, respectively, where the results from the original model by Xu et al. [31] and the proposed model are reported and compared. As it can be observed, the loading is reversed before full trasformation is achieved, determining a clear different behavior during reverse transformation for the two models. A large amount of irrecoverable TRIP strain is accumulated from the first cycle to the 25th cycle. Moreover, due to the specific loading history, TRIP strain does not stabilizes and a typical ratcheting behavior is evident.
The second pseudoelastic test involves five tension cycles at a temperature of 360 K. First, a stress of 750 MPa is applied and unloaded back to 350 MPa. Then, the applied stress varies cyclically between 350 and 560 MPa (see Fig. 2b), inducing PT and TRIP. Figure 2d and f compare the stress-strain and strain-loading step parameter curves, respectively, obtained with the two models. The loading history causes full transformation during the first loading up to 750 MPa, followed by PT during subsequent cycles. Models' responses overlap for the major transformation loop, while differ for the loops caused by PT.
The third pseudoelastic test involves two non-proportional cycles at a temperature of 360 K. The biaxial loading consists in a square path, where the axial components 11 and 33 of the applied stress vary between 0 and 550 MPa (see Fig. 3a). Figure 3b compares model formulations in terms of true strain component 11 versus component 33 curve. The path induces only PT and the difference between the two model formulations is clear when reverse transformation takes place. It is also noted that the considered loading path induces a limited, but different, amount of TRIP strain in the two model predictions.

Uniaxial Actuation Tests
Different actuation tests are simulated considering a single, eight-node hexahedral finite element in order to demonstrate the response of the constitutive model. The tests are simulated by applying a tensile stress of 200 MPa at constant temperature of 438 K (above critical A f temperature) to ensure that the material is initially at the pure parent phase of austenite, followed by thermal cycling between a maximum and a minimum value under a 200 MPa constant stress. For all these tests, material parameters from Set 2 of Table 1 [31] are adopted, related to Ni 49:9 Ti 50:1 (at.%) SMA. Since the coefficient of thermal expansion for NiTi is calculated to obtain a very low value ( $ 10 À5 ) and because its effect is negligible compared to phase Note that r b is purposely being put as zero to simplify the simulation as the two-way SME is not simulated here transformation, it is assumed to be equal to zero. As for pseudoelastic tests, the material properties selected for this example were used to verify that the unified model was able to predict the material behavior presented by Xu et al. [31] under complete transformation cycles.
The first actuation test is performed by applying three thermal cycles between 303 and 438 K to ensure full transformation, followed by two thermal cycles between 303 and 373 K to cause PT (see Fig. 4a). The second actuation test is performed by applying two thermal cycles between 303 and 438 K to ensure full transformation, followed by three thermal cycles between 340 and 438 K to cause PT (see Fig. 5a). The third actuation test is performed by applying two thermal cycles between 303 and 438 K to ensure full transformation, followed by three thermal cycles between 340 and 380 K to cause PT (see Fig. 6a). Figures 4b-c, 5b-c, and 6b-c show the strain-temperature and strain-loading step parameter curves for the three tests, respectively, where the results from the original model by Xu et al. [31] and the proposed model are reported. For all the tests models' predictions are equal for major transformation branches and in agreement with the results published in [31], while differ for minor loops. Particularly, the first, second, and third test show model response for PT during heating, cooling, and both heating-cooling, respectively. In all cases, the proposed model shows a gradual and smooth PT description compared to the original model.

Tube Actuator
To demonstrate the predictive capabilities of the developed model, a more complex actuation test is simulated, i.e. a SMA torque tube to be used as rotational actuator. The torque tube has a inner radius r equal to 3 mm and a thickness t equal to 0.1Ár mm. To reduce computational cost, only a tube segment of length L = 2/3Ár is studied, as shown in Fig. 7. The segment is meshed with 640 eightnode hexahedral finite elements (denoted as C3D8 in the Abaqus library). For the simulation, we adopt material parameters from Set 2 of Table 1. The tube segment is fixed at one end, while at the other end a torque load is applied. Figure 7 shows the analyzed geometry and the applied boundary conditions. The test is simulated by applying, first, a torque load of 3 N m at constant temperature of 410 K to ensure that the material is initially at the pure parent phase of austenite; then, 50 thermal cycles between 300 and 410 K are applied to cause full transformation, followed by five thermal cycles between 340 and 380 K at constant applied torque load to ensure PT. Figure 8a and c show, respectively, the shear straintemperature and twist angle-temperature curves obtained with the proposed model. The first 50 thermal cycles are applied to train the material and induce full transformation and clearly shown that the actuator undergoes large strains and rotations. Accordingly, models' responses are completely matching for these major loops. The latter five cycles induce PT in the tube actuator. The comparison between the proposed and original model for the PT cycles in terms of shear strain and twist angle is shown in Fig. 8b and d. During PT, TRIP strain slightly increases. The comparison between the two models again shows a clear difference between the predicted minor loops. The stabilization of the TRIP strain leads to PT loops demonstrating the reversal memory property of SMAs.  The experimental apparatus consists of a MTS Insight testing machine, which is equipped with a MTS loadcell with a 30 kN load capacity. After the specimen was aligned with the grips and positioned, a high temperature Epsilon extensometer of a 1-inch gauge length was attached on the specimen. Moreover, two type K thermocouples were also attached on the specimen to measure its temperature during the experiment. The specimen along with the grips was enclosed to a properly insulated Thermcraft thermal chamber. A controller was responsible for the determination of the heating and cooling rate. Heating of the specimen was performed by heating the air inside the chamber, while cooling was performed by allowing cold liquid nitrogen vapor to enter the chamber. An electric valve opened during the cooling phase and the flow of liquid nitrogen was manually controlled through a valve regulator. Experimental measurements from the thermocouple and the extensometer were transferred into a Windows 7 workstation where they were stored.
The experimental test consisted in applying a tensile stress of 395 MPa at constant temperature of 493 K (to ensure fully austenite state), followed by thermal cycling at constant applied mechanical load. Particularly, one major thermal cycle between 493 and 373 K was performed to induce the full transformation from austenite to martensite and back, followed by seven minor cycles between 493 and 427 K to induce PT during cooling phase.  To simulate this test, we adopt material parameters from Set 3 of Table 1, calibrated on experimental data related to the major hysteresis cycle. Experimental data related to the minor loops are used to validate the model, since the behavior under complete transformation cycles has been already verified. As it can be observed in Fig. 9, both experiments and numerical simulations highlight a small amount of inelastic strain during partial cycling, that accumulates as the number of cycles increase. Yet, when the reverse transformation is fully completed, this strain seems to not affect the strain level exhibited at the pure phase in the experimental data. The proposed model formulation (see Fig. 9a) shows good agreement with experiments compared to the original model formulation (see Fig. 9b).

Conclusions
This paper presented a unified modeling framework for the simultaneous description of TRIP as well as PT characteristic in SMAs subjected to cyclic thermo-mechanical loading.  [31]. Following that, the prediction of the baseline line model for partial phase transformation cycles are compared with ones predicted by adopting the modified hardening function proposed by Karakalas et al. [15]. To demonstrate the difference between the two different approaches three minor loop cases were considered: (i) partial transformation reversal is initiated on the major heating branch and evolves till fully martensitic state is reached, (ii) partial transformation is initiated on the major cooling branch and heating continues till the material is at the fully austenitic state and (iii) starting from the major cooling branch partial transformation loops are performed without reaching either purely austenite or martensite phase. Particularly, the comparison with the original model published by Xu et al. [31] has demonstrated the successful capability of the model in accounting both for PT reversal point memory and the associated memory wipe-out in the case of PTs coupled with the generation and saturation of TRIP strains, without the need of additional model parameters. Moreover, the comparison with experimental data under uniaxial conditions demonstrates that the proposed modeling approach is promising for the description of these phenomena. It is worth highlighting that a comprehensive experimental and numerical investigation of the coupling between the two described features under both proportional and non-proportional loading is still missing in the current literature. Accordingly, future work will focus on the validation of the proposed model on an ad-hoc experimental campaign. Haghgouyan during the experimental tests is greatly appreciated. The developed user-defined material subroutine (UMAT) will become available upon request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
Funding Open access funding provided by Università degli Studi di Pavia within the CRUI-CARE Agreement.

Appendix A: Summary of Preliminaries for Model Formulation
Considering a deformable body B, the deformation process of a material point P from reference (undeformed) configuration to current (deformed) configuration can be described by using the deformation gradient second order tensor: Here, X and x denote the position vector of P, respectively, in the reference configuration at time t 0 and in the current configuration at time t, U and V are the right and left stretch tensors, respectively, and R is a proper orthogonal rotation tensor. The right and left Cauchy-Green deformation tensors are then, respectively, defined as: The velocity gradient L can be additively decomposed into a symmetric part, i.e., the rate of deformation tensor D, and an anti-symmetric part, i.e., the spin tensor W, as follows: where v ¼ dx=dt ¼ _ x is the velocity field of P and The material and spatial Hencky (logarithmic) strain tensors H and h are, respectively, defined as: where U ¼ x À X is the displacement vector and Grad is the gradient. Based on the works [50][51][52][53], it is proven that the corotational rate of the logarithmic strain h associated with the so-called logarithmic spin X log is identical to the rate of deformation tensor D: where: where R log is a second-order rotation tensor such that Based on Eq. (A.8) and (1)