Vectorized MATLAB Implementation of the Incremental Minimization Principle for Rate-Independent Dissipative Solids Using FEM: A Constitutive Model of Shape Memory Alloys

: The incremental energy minimization principle provides a compact variational formulation for evolutionary boundary problems based on constitutive models of rate-independent dissipative solids. In this work, we develop and implement a versatile computational tool for the resolution of these problems via the ﬁnite element method (FEM). The implementation is coded in the MATLAB programming language and beneﬁts from vector operations, allowing all local energy contributions to be evaluated over all degrees of freedom at once. The monolithic solution scheme combined with gradient-based optimization methods is applied to the inherently nonlinear, non-smooth convex minimization problem. An advanced constitutive model for shape memory alloys, which features a strongly coupled rate-independent dissipation function and several constraints on internal variables, is implemented as a benchmark example. Numerical simulations demonstrate the capabilities of the computational tool, which is suited for the rapid development and testing of advanced constitutive laws of rate-independent dissipative solids.


Introduction
New experimental techniques enable more thorough investigation of the complex response of materials to mechanical loading, which opens space for the development of more elaborate material models and physical simulations. On the macroscopic (continuum thermodynamics) level of modeling, such development involves deducing more complex constitutive laws which complement the fundamental balance laws and side conditions (boundary and initial) so that the response of a material body in time to external stimuli can be determined via solving evolutionary boundary value problems. For the development of complex constitutive laws characterizing the materials and material systems, various thermodynamic frameworks have been developed in the literature [1][2][3][4][5]. They allow for the formulation of a wide range of models in a very concise and consistent way. The incremental energy minimization approach can be considered a compact variational formulation of the evolutionary boundary value problem for models of rate-independent dissipative solids which were derived within such frameworks [6][7][8][9]. Let us note that rate-independent processes are invariant under a change in time scale [10].
To obtain a time discretization of an evolutionary boundary value problem which is suitable for implementation in finite-element software, we introduce a partition of the time interval from time 0 to T in the form 0 = t 0 ≤ t 1 ≤ · · · ≤ t N = T . In the time-discrete setting, the response of the system at time t n+1 can then be determined by solving the incremental (energy) minimization problem inf α J τ (α, α n ) (1) with α representing the set of all thermodynamic descriptors of the system (e.g., thermodynamic variables vital for the system's evolution). The superscript τ denotes the time-discretized counterparts of time-continuous functionals of the corresponding weak formulation, and the subscript n denotes the values of the previous time step. The (Lagrangian) functional J τ usually combines three physically motivated terms: where the incremental energy functional is E τ , the incremental dissipation functional is D τ , and the incremental external work functional is P τ . Let us note that the minimization can be subject to some additional internal constraints (kinematic and physically-based). For (hyper)elastic materials, D τ disappears, α reduces to the displacement field (u), and E τ involves only gradient(s) of u (cf. [11]). Through different mathematical forms, the dissipation functional allows one to describe different types of path-dependent material responses. Rate-independent (also termed activated) dissipative processes are characterized by a positively one-homogeneous dissipation (i.e., such a unction d(α) satisfies the relation d(sα) = sd(α) for all positive scalars s) and hence, the dissipation functional D τ is inherently non-smooth. Maybe the most common example is the conventional model of (isothermal) von Mieses plasticity, with displacement and inelastic strain as two thermodynamic descriptors appearing in α (e.g., [4]).
The rate-independent problem in Equation (1) can be resolved via different numerical approaches. Often, the particular form of Equation (2) allows splitting the minimization problem into a nested form and solving a sequence of "structural" nonlinear minimization problems and "material" nonlinear minimization problems [5]. The form of the "structural" problem formally corresponds to the principle of minimum potential energy, and the solution procedure is analogous to resolving an elastic body problem. The process of resolving the "material" problem is often termed a state update procedure. Specific numerical techniques have been intensively developed even for the most common constitutive laws [12]. Many finite-element suites allow straightforward implementation of the nested form-they effectively utilize the alternative minimization approach-and require a material tangent stiffness operator for the interconnection between the two problems (see, for example, the examples for smart materials in [13][14][15]).
On the other hand, the variational structure of Equation (1) allows the direct application of optimization methods borrowed from mathematical programming. Such an approach may be particularly attractive for models of materials with multiple, strongly coupled dissipative processes, such as coupled plastic flow, phase transformation, and micro-damage in so-called TRIP steels [16], coupling of phase transformation and plasticity in NiTi shape memory alloys [17], or the coupling of damage and plasticity in sensitive clays [18]. Indeed, in such (rate-independent) systems, the evolution cannot be split into mutually independent processes and treated separately, such as by conventional active set search strategies [18]. Instead, only a single "global yield function" (borrowing the conventional plasticity terminology) driving the evolution can be derived, and the conventional numerical treatment may then become elaborate and conditional upon the particular mathematical form of the model [17].
Recently, Moskovka and Valdman [11] developed and successfully verified the fully vectorized implementation of hyperelastic constitutive models into the finite element method (FEM) through mathematical optimization. Their approach relies on a concurrent (loop-free) evaluation of (both linear and gradient) energy contributions to the global energy functional and efficient evaluation of energy gradients (needed in optimization) by employing the concept of nodal patches. The applied vectorization has already been shown to be efficient and flexible in developing versatile FEM implementations [19,20] of nodal and edge basis functions and assemblies of stiffness matrices (cf. [21]).
Inspired by their work, the present contribution aims to further extend that concept for the energy functional of the type in Equation (2) (i.e., with a non-negligible dissipative contribution). As a representative example, we employ the previously developed constitutive model for the phase transformation and reorientation processes in shape memory alloys [22], which belongs to the class of models with strongly coupled rate-independent dissipation processes. This requires several extensions and modifications of the computational framework in [11]: (1) The set of thermodynamic descriptors of the system is enriched by internal variables (and temperature) fields. (2) Both the displacement field and fields of internal variables enter the minimization process. (3) Since the time evolution of the system becomes path-dependent (i.e., the evolution strongly depends on the particular loading history of the system), the values of the descriptors from the previous time step, α n , must be stored. (4) Finally, due to the activated nature (rate-independence) of the dissipation processes, the objective function becomes non-smooth (although it remains convex). Moreover, concerning common loading modes of shape memory alloys, temperature parametrization and the Neumann boundary condition must be implemented.
The main aim of this work, however, goes beyond the mere implementation of the particular constitutive model. Our activity targets the development of a versatile tool for the development and testing of new constitutive laws for materials with multiple dissipative processes. Such efforts often require an iterative approach consisting of repetitive tuning of several contributions to energy or dissipation functionals and comparing the basic experimental datasets with computational simulations performed with the model. The numerical implementation is usually the most laborious part of the process, especially when the resolution procedure requires tailored coding (e.g., related to active set search strategies). In this sense, an in-house built FEM-based computational tool-which combines versatile variational approaches, universal optimization methods, and a user-friendly programming environment-may often provide a reasonable alternative to "heavy-duty" commercial engineering computation platforms. Indeed, straightforward and rapid adaptation of the code to a new constitutive law outbalances the limited computational efficiency (due to the employment of uncustomized optimization algorithms) and user comfort. The present work constitutes a first step for recasting this concept into a viable paradigm. This paper is organized as follows. The constitutive model is briefly introduced in Section 2, where the mathematical formulation of the incremental energy minimization problem is also provided. Section 3 deals with the numerical implementation, particularly the optimization approach based on finite element discretization. Computational examples and conclusions are presented in Sections 4 and 5, respectively.

Constitutive Model of Shape Memory Alloys and Incremental Energy Problem Formulation
We briefly summarize the core constitutive model for shape memory alloys originally introduced in [22,23] and further refined in [24][25][26]. Let us note that more experimental details on this important class of functional materials can be found in [27], for example, and a great variety of constitutive models tackling the peculiarities of their thermomechanical response can be found in the literature (e.g., the recent works of [28][29][30][31][32][33]).
The common additive decomposition of the total strain ε to the elastic part ε el and the transformation-related part are employed in the small strain realm [28,30]: Scalar ξ denotes the volume fraction of martensite, and ε tr is the mean transformation strain tensor of martensite, which is symmetric and traceless; in other words, it holds that ε tr 11 + ε tr 22 + ε tr 33 = 0 (4) for the diagonal components ε tr ii , and ε tr ij = ε tr ji for the non-diagonal components. The natural constraint on ξ is complemented by a directional constraint on the transformation strain: where the material function · takes a specific form (see [22]) so that the tension-compression asymmetry and material anisotropy are recovered: where I 2 (x) = 2 3 x and I 3 (x) = 4 det(x) I 2 (x) 3 , with x being a tensor. Two particular functions constitute the core of the model. The first one is the energy function (thermodynamic potential) f T , consisting of the elastic part (the first two terms), the chemical part (the third and fourth terms), and the reorientation-hardening part (the last term): where tr(x) and dev(x) stand for the trace and the deviatoric parts of a tensor x, respectively. The second one is the dissipation function d T , in which both transformation and reorientation of martensite are considered: Let us note that the superscript T emphasizes that the temperature T is considered only as a prescribed (albeit variable) parameter in quasistatic loading regimes, which is the focus of this work. The implicit dependence of descriptor fields on space coordinates is omitted for brevity. Function q(T) in Equation (7) is purely temperature-dependent, and hence its particular form is not relevant for minimization. The table in Section 4 provides a brief description of all model parameters (see [22] for further details).
In the time-discrete setting, the rate of the constitutive state is considered to be constant in the time increment τ. Hence, the one-homogeneous dissipation function d T can be reformulated into an algorithmic expression as follows (see [23] for details): Now, we consider a uniform time discretization with τ = T /N and t 0 = 0 so that t n+1 = (n + 1)τ, n ∈ {0, . . . , N − 1}. We can specify the incremental functionals for a body from shape memory alloys as and construct the complete functional according to Equation (2): Here, Ω ⊂ R 3 is the geometric representation of the physical body with (assumed Lipschitz) boundary ∂Ω, and Γ N ⊂ ∂Ω represents its subset where the Neumann boundary condition is applied. The terms F vol n+1 , F surf n+1 , and T n+1 represent the corresponding time discretization of the prescribed volumetric forces, surface forces, and temperature, respectively.
The reformulated incremental energy minimization problem in Equation (1) determining the (discrete) time evolution at t n+1 then reads as on the proper functional spaces also incorporating the Dirichlet boundary condition. The reader is referred to [23] for details regarding the rigorous mathematical justification, treatment, and solution existence results of this formulation, relying on the "energetic solution concept" deeply elaborated upon in [9]. Let us note that due to the one-homogeneous nature of the dissipation function, time plays a role of a mere parameter in the response of the material, as expected for rate-independent constitutive laws.

Numerical Implementation
Within the incremental minimization approach, the time evolution of the material system is fully determined by resolving a sequence of time-incremental minimization problems given by Equation (14). We combine the finite element method (FEM) with classical optimization methods in order to solve each problem of this sequence efficiently (cf. [8]).

Minimization Strategy
In contrast to the common nested procedure, where the minimization problem in Equation (14) is split into two subproblems and resolved via alternating minimization [5], we stick to the monolithic approach (i.e., minimizing all control variables at once). Thus, we avoid the procedure of transferring the results from one subproblem to the other, involving the construction of the material Jacobian matrix (material tangent modulus) and its incorporation into the structural Jacobian matrix.
As shown in [23], each single-minimization problem is convex, the objective functional is naturally non-smooth (i.e., not continuously differentiable due to asymmetry with respect to the sign ofξ and the Frobenius norms in Equation (8)), and the design variables ξ and ε tr are constrained by the inequalities in Equation (5). Hence, after finite element discretization, we face a non-linear, non-smooth constrained convex minimization problem, which is tractable by local optimization methods.
Inspired by [22], and with no real impact on the physical background of the model, we further restrict the internal variables as follows (cf. Equation (5)): Indeed, from the thermodynamics point of view, the system is always a phase mixture at transformation-relevant temperatures, even though the volume fraction of the unfavorable phase may be macroscopically negligible.
We then impose a non-linear smooth transform between these bounded sets and suitable unbounded ones (see the details in Appendix A for details) so that respective minimization problems can be reformulated as with unconstrained variables ξ, ε tr . In this way, it is possible to employ lightweight unconstrained minimization methods instead of constrained ones. Let us note that to circumvent the common troubles linked with setting the correct values of the internal variables at the initial point of calculation ξ 0 , ε tr 0 (cf. [34]), a specific "zero increment" computation is performed first (see details in Appendix B).

Details on Implementation in MATLAB
The complete computational tool extends the FEM codes of [11] related to a 3D minimization of hyperelastic energies. The complete code is freely available at (accessed on 21 October 2022) https://www.mathworks.com/matlabcentral/fileexchange/119538 for downloading and testing. The constitutive core of the code of [11] was modified because the determination of the material response now additionally depends on the fields of internal variables and temperature. Not only do the internal variable fields ξ and ε tr enter the minimization process as the subject of minimization (together with the displacement field), but the distribution of them in Ω in the previous time step also appears in non-linear terms in D τ SMA . Such a dependence on the values in the previous time step is inevitable in dissipative models and imposes their "path dependence" (cf. linear dependence of the displacement field in the previous time step in P τ SMA , which effectively does not enter the minimization process). Hence, internal variables must be stored within the code, and their values from the previous step must be made available for the minimization procedure. In contrast, no particular treatment is needed for the temperature field, since the distribution of the temperature is prescribed. Thus, it is a priori known within each time step and each element of triangulation. Let us note that the uniform temperature distribution within the body is assumed in the present implementation; the necessary modifications toward a full thermally coupled model were elaborated upon in [26].
The computational domain Ω ⊂ R 3 is approximated by its triangulation T into closed tetrahedral elements in the sense of Ciarlet [35]. The displacement field u is then approximated in components in the space of globally continuous and piecewise linear nodal basis functions P 1 (T) defined in T, which is the lowest-order polynomial choice possible. Consequently, the small strain tensor ε(u) is approximated as a piecewise constant function P 0 (T). Thanks to Equation (4), five independent components of the symmetric traceless tensor ε tr and one scalar ξ add up to six independent design variables, which are also approximated in P 0 (T).
The code heavily utilizes vectorized evaluations of the gradients ∇u for the prescribed triangulations T inherited from the previous work [11]. All necessary gradients of basis functions in P 1 (T) are precomputed and stored in a structure-type data object "mesh" together with the geometrical properties of the triangulation. An additional structure "patches" is utilized for the evaluation of the gradient ∇ u J τ SMA . It contains information about the nodal patches (i.e., about sets of tetrahedral elements adjacent to every node of the triangulation, including gradients of the corresponding basis functions defined over them). Since in each computation increment the values of the internal variables from the previous increment are needed, an additional structure "varstruct" is constructed to store them. It also contains the prescribed value of the temperature and (density of) the surface forces, which can evolve. The values of the internal variables, stress, strain, and other parameters resulting from the computation are stored in the structure "output".
As an extension of [11], the surface traction term Γ N F n+1 surf (u − u n ) dS is included. We implement the simplest case F n+1 surf = g n+1 surf ν, where the vector ν denotes the normalized outer normal to the boundary part Γ N , and g n+1 surf is a given (density of) constant scalar surface force per unit area. Then, the vectorized implementation in the spirit of [36] utilizing information about the faces of Γ N and their outer normals reads as follows: Together, they fix 2 · 3 + 9 · 2 + 12 = 36 components of the displacement vector. Therefore, the remaining 27 · 3 − 36 = 45 components of the displacement vector are searched together with 6 · 48 = 288 components of internal variables. Thus, the total size of the minimization problem in Equation (14) is equal to 45 + 288 = 333. Figure 1 (right) illustrates the corresponding Hessian sparsity pattern. The blue dots mark the nonzero entries of the square block related only to 45 displacement degrees of freedom, and the black dots form the lines of the nonzero entries in the square block related to internal variables only. The clusters of red dots denote the nonzero entries in the remaining two blocks; they indicate the coupling between the components of displacement and the internal variables.
For the minimization, we employ the trust region optimization method available in the MATLAB Optimization Toolbox. As demonstrated in [37], providing additional information in terms of the gradients and sparsity pattern of the Hessian matrix (i.e., positions (indices) of its nonzero entries) substantially accelerates the computations. Numerical determination of the gradients (with a difference scheme) needed in optimization requires relatively little additional effort, and thus, it can be particularly advantageous for the development and refinement of constitutive laws of new materials with complex microstructural processes, where the quest for a particular form of the constitutive law is the principal endeavor.
The finite element discretization directly affects the construction of the gradient and the sparsity pattern. Whereas internal variables are approximated in P 0 (T), there are hence no inter-elemental links, and P 1 (T) is used for components of the displacement field. Thus, a nodal perturbation directly modifies the computations for all adjacent tetrahedral elements. The vectorized evaluation of the gradient components was developed and described in detail in [11], and the full vectorization of the constitutive core implemented in this work are thus vital for the efficient integration of gradient-evaluating minimization methods.

Computational Examples
In this section, we present three computational simulations to demonstrate the capabilities of our computational tool. All computations were performed on a MacBook Air (M1 processor, 2020) with 16 GB of memory. Table 1 presents the values of the material parameters inspired by [38] and used in all computational examples in this section, where D is set to the fourth-order identity tensor in Equation (6).
The first two subsections validate the constitutive law in simple loading modes common for shape memory alloys. The simulations were performed on a simple computational cubic domain of a unit length of 1 cm, which was subjected to variations in temperature, and mechanical loading was applied on the top face. The displacements on the bottom and two side faces were fixed in directions parallel with the corresponding face normals (to avoid rigid body motion). The domain was divided into eight equal cuboids, each consisting of six equal tetrahedral elements. Due to the invoked uniform loading modes, the results within all the elements were equal within the numerical error.
The example in the last subsection deals with a more complex (inhomogenous) mechanical state of the material body to demonstrate the capabilities of the complete computational tool.

Example 1: Uniaxial Tension-Compression Tests at Different Temperatures
First, we present the results of uniaxial tests performed via Dirichlet-type boundary conditions. Figure 2 presents the results of the simulation of the uniaxial tension-compression tests at five different constant temperatures, namely −50 • C, −25 • C, 0 • C, 25 • C, and 50 • C. Each simulation was divided into four steps (distinguished by "relative time" in the plots): (1) the temperature of the body was changed from room temperature (25 • C, or above A f ) to the desired one under no change in displacement, (2) the corresponding displacement of nodes on the top face of the cube was increased up to 0.06, (3) the displacement was decreased to −0.05, and (4) the displacement was again increased to zero. Except for the first step, the temperature was held constant (and uniform).
The results in Figure 2a show the evolution typical for the shape memory alloys: a pseudoplastic type of response at low temperatures (−50 • C and −25 • C, with coinciding graphs), superelastic type of response at high temperatures (25 • C and 50 • C), and an intermediate response in between (0 • C). This interpretation is backed by the plot of the evolution of the volume fraction of martensite within the steps in Figure 2c. For the two lowest temperatures, the almost full martensitic state was reached after cooling (end of step 1, i.e., relative time equals 1), whereas for the rest, transformation to martensite was induced after loading was applied (step 2). Since 0 • C was below the temperature needed for the initiation of the reverse transformation to austenite, most of the martensite remained stable even during the unloading stages of the test (in steps 3 and 4). This contrasted with the two higher temperatures, where austenite was recovered during unloading. Figure 2d shows the evolution of the relevant component of the transformation strain tensor, where the difference between the minimum and maximum reached values imposed via Equation (5), and Equation (6) resulted in the pronounced tension-compression asymmetry in Figure 2a. For instance, the critical stress for the transition to martensite was higher, and the plateau strain was lower in compression compared with the tension in the superelastic regime. The complex processes linked with energy storage and dissipation gave rise to the patterns of J τ SMA -functional minimum in Figure 2b, which are entirely distinct from the quadratic ones known in Hookean elasticity.

Example 2: One-Way Shape Memory Effect and Martensite Stabilization Effect
To demonstrate the ability to employ the Neumann type of boundary condition, we simulated the so-called one-way shape memory effect, an important phenomenon utilized in applications [27]. The cube was cooled down to −25 • C first and then loaded to 600 MPa (via imposing the Neumann condition on the top cube face), unloaded via complete removal of the stress, and finally heated to room temperature (25 • C). The results in terms of the stress-strain-temperature space and the evolution of the volume fraction of martensite with the temperature are depicted in Figure 3a,b, respectively. After unloading (at −25 • C), the material remained martensitic, with the strain reaching almost 4%. Only when it was heated enough did it transform back to austenite and return to the original configuration (shape). Hence, the strain was released, as observed in the experiments in [27].
For comparison, Figure 3b also includes the situation, in which the interlude of loading and unloading was omitted (i.e., the material was just cooled down and heated up without straining). It is worth noting that the temperature interval for the reverse transition in the one-way shape memory effect (approximately between 10 • C and 20 • C) significantly shifted upward with respect to the stress-free thermal cycle (0 • C and 10 • C, cf. M s and M f in Table 1). This temperature shift was due to the so-called "martensite (mechanical) stabilization", and it has been experimentally well-documented [39][40][41]. In the constitutive model, this is due to the particular form of the dissipation function in Equation (8), which also recognizes the contribution linked with the reverse transformation of reoriented martensite.  stress-free one-way SME (b) Volume fraction of martensite over temperature corresponding to (a) (red line), and in a stress-free thermal cycle (blue line).

Example 3: Bending of a Shape Memory Alloy Beam
The performance of the complete FEM computational tool is now demonstrated in the example of the bending of a beam. Beams manufactured from shape memory alloys have been intensively investigated due to their superb actuation and energy harvesting capabilities [42][43][44]. Although analytical models can be very useful for quick concept validation and preliminary design studies [45][46][47], full-scale FEM simulations are often indispensable in detailed studies [42,44,48,49], since they are not constrained by assumptions on a particular geometry, material response, boundary conditions, etc. The simulated problem mimics a simple beam with a rectangular cross-section (with its width being half of the height) loaded at both ends and supported in the middle so that bending is invoked. Thanks to the expected structural symmetry of the situation with respect to two (out of three) symmetry planes of the beam, it was enough to model only one-fourth of the structure and hence fix the displacement of nodes adjacent to those planes of symmetry in a direction normal to the planes. The particular geometric model of the symmetry segment (length of 0.4 cm, height of 0.1 cm, and width of 0.025 cm) together with applied triangulation is depicted in Figure 4. The structured mesh with 20 × 1 × 10 identical small cuboids each containing 6 tetrahedral elements was used to capture primarily the details in the x-z symmetry plane. Nodes on the plane x = 0 cm and y = 0.025 cm were constrained according to the symmetry requirements, nodes satisfying x = 0 cm and z = 0 cm were also constrained in z displacement (avoiding rigid body motion), and nodes satisfying x = 0.4 cm and z = 0 cm ("bottom end edge") were linearly prescribed uniform displacement in a direction z so that they reached position z = −0.06 cm in 6 increments. A uniform temperature of 25 • C was considered during the whole simulation, and hence superelastic behavior was induced. Let us remind the reader that the current work adopted small strain and isothermality assumptions; a discussion of the influence of such premises on the results of computational simulations of shape memory beams can be found in [49]. The deformed configurations overlaid with contour plots showing the distribution of the volume fraction of martensite in each computational increment are shown in Figure 5. A continuous transformation was documented, starting from both the upper and lower sides of the beam and penetrating toward the central part containing the neutral plane. Let us note that validation of the SMA model via comparison with the experimentally determined distribution of phases has been performed recently [24,50].  Figure 6a shows the distribution of the diagonal component of the stress tensor in the x direction in the reference configuration in the last increment, (i.e., corresponding to Equation (6) in Figure 5). Due to the tension-compression asymmetry, its maximum tension was approximately 2/3 of the absolute maximum in compression, and the neutral plane shifted toward the compressed part of the material (cf. analytical computations in [46,51]). For comparison, the same computational simulation was performed at −25 • C (i.e., in the pseudoplastic regime of the material), and the distribution of the diagonal component of the stress tensor in the x direction is shown in Figure 6b. Both plots use the same color coding, which allows for direct comparison of the numerical values. Clearly, thanks to the lower net stiffness of martensite (which was reoriented during loading), the stress maxima in this were lower than in the superelastic case in Figure 6a. Such strong temperature dependence of the structural stiffness is employed in many applications [52].

Conclusions and Future Outlook
This work presents a successful implementation of the incremental energy minimization approach for rate-independent dissipative solids within a finite element framework. A model of shape memory alloys rigorously formulated in [23] served as a benchmark example of such a type of constitutive law. It features a strongly coupled, rate-independent dissipation function. The solution of an evolutionary boundary value problem of a solid body is constructed as a sequence of solutions of adjacent time-incremental boundary value problems in the spirit of the energetic solution concept for rate-independent materials [9]. Each of these fixed-time problems is considered to be a genuine minimization task (Equation (1)) discretized in space via the finite element method and resolved by conventional optimization tools. In this respect, the chosen monolithic approach represents a consistent alternative to the staggering minimization strategy employed in [23] and follows a thus far unexplored pathway offered by the energetic solution concept.
The vectorized implementation into MATLAB also presents the first step in a longterm activity focused on building a lightweight FEM-based computational tool tailored for rapid development, preliminary testing, and validation of new constitutive laws of dissipative materials with complex microstructural processes. Here, the motivation comes from material science and engineering, where the rapid progress in the manufacturing and characterization of new materials is often not fully exploited by the constitutive models. Among other factors, this is because the process of validation of new constitutive laws utilizing tailored implementation into computationally-aided engineering software suites is usually laborious and time-consuming. Employing a versatile FEM-based tool provides a viable alternative, especially in the initial stage of constitutive law development. Of course, when the form of the constitutive law is finally established, its tailored implementation in such suites (including developing its own constitutive subroutines) is inevitable for the resolution of more complex boundary problems (e.g., with geometry requiring adaptive remeshing or involving additional engineering phenomena as contacts). Following this idea, we adopted and substantially modified the vectorized code from [11] so that it could be used for both hyperelastic and rate-independent constitutive laws. The modifications include the introduction of internal variables and their history into the minimization process, adding a solution parameter (representing the temperature field), or handling the Neumann boundary conditions. The capabilities of the whole computational tool were tested on a problem of bending a shape memory alloy beam. The version presented in this work can be downloaded for free from a dedicated website.
The benefits of the applied approach come hand in hand with some limitations. The first one is related to the variational formulation. Although the class of constitutive models complying with Equation (2) is wide, a special mathematical and numerical treatment is needed for some types of models (e.g., the concept of bipotential [53] for materials with non-associative laws or the homogenization techniques for heterogeneous materials [54]) (see also the mathematical background in [9]). Second, the dimension of the minimization problem in the monolithic formulation is naturally higher than in the case when the problem is split in the nested formulation, which poses a challenge for computational resources. Third, applying universal minimization methods is often less efficient than employing customized algorithms and may lead to prohibitive computational time consumption for more complex evolutionary boundary problems. Whereas the first issue may require completely new numerical development, relatively milder modifications such as parallelization or algorithms of non-smooth optimization can help to at least partially remove a burden from the other two constraints.
Hence, one direction of planned future refinements of the numerical implementation resides in improving the computational speed. The time needed for the computation of gradient inputs to the optimization routines could be reduced, such as by including the methods of automatic differentiation (cf. [55]). Alternatively, one can even abandon the intrinsic MATLAB optimization tools and employ some progressive optimization methods tailored to the inherent structure of the rate-independent constitutive laws as, for example, fast iterative shrinkage-thresholding algorithms adapted for non-smooth functions [56]. Reformulation of the models within the finite strain (large deformation) theory poses another challenge, both from the modeling and from the implementation (vectorization) points of view. Inspiring attempts involving these aspects have been carried out (e.g., in [13,48,57]). 15 × 1 × 8 cuboids which were mutually identical cuboids, and the fine mesh consisted of 24 × 1 × 12 cuboids. We plotted the values of the dominant component of stress (also visualized in Figure 6) at the end of loading in two important sets of nodes. Set A consists of all nodes with fixed coordinates y = 0.025, z = 0.1 in the reference configuration, and set B consists of all nodes with fixed coordinates x = 0.0, y = 0.025 in the reference configuration (cf. Figure 4). Hence, the two sets correspond to the top side of the beam and the plane of bending symmetry (most bent cross-section) (i.e., the left and top edges in Figure 6a). Let us note that the number of nodes in the sets increased with the refinement of the mesh, and the value in each node was obtained as a weighted average of the values in all node-adjacent elements. Figure A1a shows the spatial variation of the normal component of stress parallel with the x axis (in the reference configuration) on the top side of the bent beam for the three meshes. One can observe a very good mutual match with the highest variation close to the plane of bending symmetry, where the phase transformation occurred dominantly. For x > 0.05, the difference between the coarse and fine meshes did not exceed a few MPa. Figure A1b compares the same stress component values on the plane of bending symmetry, where a huge variation in the volume fraction of martensite developed from the bending (cf. the last plot in Figure 5). Whereas the variation of the maximum and minimum (i.e., values on the outer surfaces of the beam) was about 3% only, a higher variation can be observed in the central part (i.e., close to the neutral fiber (more precisely, the neutral plane)). Due to the expected high gradient of phase composition, an even more refined mesh in the z direction would be needed to investigate the strong phase gradients in this region, which was, however, beyond the scope of this paper.