Abstract
Porous Shape Memory Alloys (SMAs) are of particular interest for many industrial applications, as they combine intrinsic SMA (shape memory effect and superelasticity) and foam characteristics. The computational cost of direct porous material modeling is however extremely high, and so designing porous SMA structure poses a considerable challenge. In this study, an attempt is made to simulate the superelastic behavior of porous materials via the modeling of fully dense structures with material properties modified using a porous/bulk density ratio scaling relation. Using this approach, direct modeling of the porous microstructure is avoided, and only the macroscale response of the model is considered which contributes to a drastic reduction of the computational cost. Foam structures with a gradient of porosity are also studied, and the prediction made using the fully dense material model is shown to be in agreement with the mesoscale porous material model.
1. Introduction
Shape Memory Alloys (SMAs) exhibit unusual mechanical properties such as shape memory and superelasticity, which make them very attractive for a wide range of industrial applications, from aerospace to medicine [1, 2]. For more than a decade, not only fully dense but also porous or foamed forms of SMA have been studied because of their additional benefits: low density, high permeability, and energy dissipation properties [3]. The level of foam porosity is selected as a function of the final use of a material while structural applications of porous SMA require low-to-medium porosity foams (pore volume fraction (PVF) ≤ 40%); biomedical applications, such as bone implants [4], need highly porous material with PVF of up to 70%.
The properties of porous SMAs are strongly dependent on their porous microstructure, whose length scale is much smaller than that associated with the macroscale response of the whole material. From a numerical point of view, modeling this micro/macro behavior has a tremendous numerical cost. To alleviate the complexity of the micro/macroapproach, that is, the explicit representation of the porous microstructure, researchers have chosen different routes, such as micromechanical averaging techniques [6, 7] or Unit Cell approach [8, 9]. However, those strategies are based on assumptions (low porosity, regular pore distribution, spherically shaped pores, etc.) that are not fully compatible with the biomedical application foams [5] that we are studying in the present paper.
To overcome the micro/macro numerical cost while reconciling the needs of biomedical foams, we have decided to follow the strategy adopted by Gibson and Ashby [10] by defining scaling relations. In this approach, a property of the foam material may be obtained as a function of the relative density using the following power law function: where the subscript denotes the solid, that is, fully dense, material, the density, and and the coefficients to be determined.
It is clear that this approach focuses only on the macroscale response and, when averaging the material behavior, that it explicitly excludes the microscale effects (stress concentration, cracks, local plastic deformation, etc.) appearing in a nonuniform, porous pattern. The scaling relations thus obtained could then be useful for designing mechanical parts such as bone implants—giving a fair estimation of their overall stress-strain behavior while drastically reducing the associated numerical effort.
The paper is organized as follows. In Section 2 we present the models of porous SMA being studied, followed by the numerical simulation results in Section 3. Our conclusions are given in Section 4.
2. Computational Setup
2.1. SMA Constitutive Model
All the computations are performed using the commercial finite element software Ansys. The SMA model used is the one implemented in Ansys [11], based on the work of Auricchio et al. [12]. The model is briefly summarized below, and the interested reader is invited to refer to Auricchio et al. [12] and references therein for more information.
The material model considers the nonlinear behavior that occurs during the reversible phase transformation from Austenite to Martensite ( in Figure 1) or from Martensite to Austenite ( in Figure 1). For the sake of clarity, we adopt hereafter the following convention concerning exponents and subscripts: , , , and refer to, respectively, Austenite, Martensite, Austenite to Martensite, and Martensite to Austenite phase transformations. It is assumed that the material behavior does not show any permanent strain and that the material is perfectly isotropic. Figure 1 depicts the five constants used to model the SMA behavior:(i) is the starting stress for the phase transformation from Austenite to Martensite; (ii) is the final stress for the phase transformation from Austenite to Martensite; (iii) is the starting stress for the phase transformation from Martensite to Austenite; (iv) is the final stress for the phase transformation from Martensite to Austenite; (v) is the maximum superelastic strain.
Note that, along with the SMA constitutive model, the isotropic elastic Austenite model parameters (Poisson’s coefficient and Young’s modulus, the latter depicted by in Figure 1) must be inserted. To complete the model presentation, , , , and moduli can also be used for the model description as shown in Section 3.2.
2.1.1. Note on the Model
Auricchio’s model as implemented in Ansys software may be considered somewhat weak, as more sophisticated models have been published since by different research groups [13, 14] and even by the Auricchio’s research group [15]. However this model presents two fundamental assets which make it very suitable for the purpose of the present work.(i)Firstly, the model focuses on the superelastic behaviour of the SMA using a low number of constitutive parameters, and it is successfully validated in a large number of works.(ii)Secondly, its direct implementation in the commercial software Ansys excludes the necessity to apply USER Material Subroutine (USERMAT) and therefore allows the reader to easily verify the results obtained in this work.
With this in mind, the use of a more sophisticated model is not justified in light of these assets.
2.2. Computational Model and Finite Element Discretization
We here define the geometrical models used for modeling SMA foams. Most of the manufacturing techniques (sintering and self-propagating) produce a nonregular porosity, as illustrated in Figure 2: the pores are randomly created in terms of their shape, orientation, and distribution.
To model such a porous material, it is clear that the numerical model geometry must also be generated in a random way to fit as much as possible the characteristics of the experimental sample. However, we will also investigate regular porous patterns in order to evaluate how close this approach can reproduce the behavior of a randomly foamed material. For this reason, both randomly and periodically distributed porosity models will be addressed, keeping in mind that the definition of the scaling relations is the main purpose of the study.
2.2.1. Randomly Distributed Porosity: RVE Approach (Mesoscale Model)
Numerical modeling of foam material is a nonroutine problem, particularly in the case of randomly shaped and randomly distributed porosity. Indeed, as we try to capture the macroscale behavior, the existence of random porosity in terms of pore shape, size, and orientation leads to micro/macroanalysis which is very expensive from a numerical point of view. It is thus appealing to define a much smaller model that is representative of the macroscale behavior: the so-called Representative Volume Element (RVE) or mesoscale model.
The determination of the RVE size has been widely studied [16, 17] and, more specifically, for titanium foam in Shen and Brinson [18]. The main idea of this approach lies in applying specific boundary conditions on a set of randomly generated unique-sized RVEs where the final property is given by averaging the results obtained for each RVE. The number of RVEs needed depends on the nature of the property to be modeled and on the porosity and the size of each RVE. However, since we limit our study to a two-dimensional model, the computational effort is reduced, and so the averaging approach is not needed. The appropriate RVE size is determined when the convergence of elastic responses is reached.
To ensure a good connectivity between the constitutive elements of the mesoscale model, we consider a hexagonal lattice where the porosity is created by connecting grains of matter (hexagonal elements) together: each element of the lattice may be considered either as an SMA grain or as a pore when filled with matter or empty, respectively. Initially all the elements on the edges of the domain are considered as matter elements, as illustrated in Figure 3 (in red). Inside this frame, “paths” of matter are created in the following way: one of the existing matter elements is selected at random, and one of its adjacent elements is assigned as an SMA element. The assignment of new matter elements is repeated until the desired pore volume fraction (PVF) is reached. To ensure the uniformity of the final pore distribution, the entire lattice domain is divided into a set of subdomains, and the random choice of an element to start a path is then restricted in a subdomain with greater porosity. Typical pore patterns obtained with the algorithm are depicted in Figure 4.
(a) Pore volume fraction of 0.3
(b) Pore volume fraction of 0.7
This model definition strategy presents several assets. (i)A good homogeneity of pore size, shape, and distribution is ensured, even in the high-porosity cases (PVF = 0.7). (ii)Each matter element of the lattice is connected to the main cluster by at least one edge. (iii)The lattice is easily meshed by dividing each hexagonal element into three equally quadrilateral-shaped elements (see Figure 5).
2.2.2. Periodically Distributed Porosity: UCFEM Approach
Using regular-shaped porosity is an attractive way to decrease the numerical cost of modeling randomly distributed porous material because it makes it much simpler for the problem to be well posed (contrary to the RVE model): for each porosity, one and only one geometry exists. Furthermore, the use of a periodic structure decreases the numerical cost of the model since only a part of this structure (as in [9]) or a small number of initial structural elements are needed.
This approach, called the Unit Cell Finite Element Method (UCFEM), allows a small part of the model to be accounted for, due to the symmetry properties of the geometry. The following unit cell geometries are retained in Figures 6(a)–6(d): “square with squared hole,” “square with rounded hole,” “hexagon with hexagonal hole,” and “circle with circled hole.” Note that the outer diameter of the “circle with circled hole” cell was squared on its left, right, bottom, and top-oriented edges to avoid geometrical singularities. Note also that this last pattern is the only geometry to exhibit a bidisperse porosity distribution in terms of shape and size: both rounded- (inside the cell) and starred- (between four adjacent cells) shaped pores appear when assembled.
(a) Square
(b) Hole
(c) Hexagonal
(d) Circle
3. Numerical Results
3.1. Quasistatic Elastic Properties
We first look at the quasistatic elastic behavior exhibited by both the porous RVE and UCFEM models during a simple tensile test: a displacement corresponding to a nominal strain of 0.05% is imposed on the top edge of the geometry modeled while a null displacement condition is imposed at the bottom. On the other edges of the geometry, symmetry conditions are imposed to model the UCFEM cluster. The model is meshed by bidimensional quadratic finite elements (8 nodes) with full integration method, and both plane stress and small deflection conditions are assumed.
As stated in the previous section, we focus on the macroscale response of the porous models. Therefore, for the sake of clarity, we always refer to the macroscale or apparent properties of the foam (if not stated otherwise). Figure 7 depicts the evolution of the elastic modulus for both RVE and UCFEM models as a function of PVF, starting from a fully dense Young’s modulus of 42 GPa.
(a) Porous mesoscale model
(b) UCFEM models
To evaluate the porous mesoscale model (RVE), ten geometries were randomly generated for each porosity type, for three different lattice sizes (100 × 100, 400 × 400, and 500 × 500) (Figure 7(a)). The solid line corresponds to the mean value obtained for a given PVF, and the error bars correspond to the standard deviation. We observe that the variance of the elastic modulus obtained depends both on the initial lattice size and on the PVF considered: the greater the porosity, the greater the deviation from the mean value. However we can note a significant reduction of the deviation, at equal porosity, between 100 × 100 and 400 × 400 lattices, as well as the slight difference observed between 400 × 400 and 500 × 500 ones. These observations allow us to assume that the results have both geometrically and numerically converged, starting from the 400 × 400 lattice size (see insert of Figure 7).
Considering the finite element model retained, we have carefully verified that the model converged both from the numerical (macroscale stress-strain response) and geometrical (randomness of the porosity) points of view. On a more numerical point of view, the final finite element model is defined by meshing each hexagonal lattice element with 3 quadrilateral finite elements. Consequently the fully dense mesoscale finite element model is defined by roughly 600,000 quadrilateral finite elements which is assumed to be fine enough to properly catch both the macroscale and microscale (local stress concentration, etc.) responses. Furthermore, in order to ensure a good accuracy of the results obtained, elements with quadratic interpolation are used.
As expected, the UCFEM results show a broad range of values depending on the cell geometry: the elastic modulus exhibited by the square cell geometry (Figure 6(a)) is several times higher than that of the circle cell geometry (Figure 6(d)), for all the PVFs considered (Figure 7(b)). Although the circle UCFEM pattern has the closest response to the RVE model, among all the cell geometries, the two approaches give rather diverging results for high porosities (Figure 8).
Since we are following the Gibson-Ashby approach [10], we should define a scaling relation for the elastic modulus as a power law function of PVF given by (1), with PVF = . For the elastic response, the best power law fits for RVE and UFCEM (indicated by the plain lines in Figure 8) are obtained with and coefficients collected in Table 1. The value corresponding to the RVE model, contrary to the UCFEM model, differs significantly from that given by Gibson and Ashby [10]. It has already been highlighted, however (by Roberts and Garboczi [19] for 3D FEM and by Wang et al. [20] with experiments), that for three-dimensional structures, effectively ranges from 1 to 3, indicating a much more complex dependence than what would result from the periodic cell theory, on which the Gibson-Ashby approach is based.
3.2. Quasistatic Superelastic Properties
As stated earlier, the results of interest are the macroscale superelastic properties of a sufficiently large RVE, allowing the definition of a scaling relation for these properties. The SMA calculations are performed in Ansys commercial software with the built-in SMA constitutive model briefly presented in Section 2.1. Since this model does not support plane stress option, the geometry computed for SMA material is identical to the elastic one with one-element thickness and a large deflection assumption. The model is then meshed using brick finite elements exhibiting quadratic interpolation (20 nodes), and full integration method is set. The model parameters are summarized in Table 2. A constant value of the Poisson ratio is used even if this parameter is different whether the material is fully Austenitic, fully Martensitic, or being transformed. It has been numerically verified that the results of a simulation carried out with a Poisson ratio of 0.3 are practically identical (deviation less than 1%) to the results of the same simulation carried out with a Poisson ratio of 0.45. Note that to reduce the number of iterations required to complete phase transformation, is set to 0.01, which is approximately 1/10 of the Ti-Ni superelastic range. The stress-strain plot of a simple tensile test for a fully dense sample is depicted in Figure 9.
As seen in the previous section, the results obtained for a 400 × 400 lattice have a deviation to the mean which is small enough to assume the geometrical independence of the results. Only one calculation is then performed for each PVF from 0.3 to 0.7, and the results are shown in Figure 10(a). For comparison, the same tests are realized using the UCFEM circle cell pattern (Figure 10(b)).
(a) Random
(b) Periodic
We observe that for SMA modeling, a UCFEM model, as compared to an RVE model, demonstrates the same trend as before: whereas at low porosity (PVF = 0.3) the UCFEM response is relatively close to the RVE response, at high porosity (PVF = 0.7), they diverge; the UCFEM model is much stiffer than the RVE one. Thus, as shown by DeGiorgi and Qidwai [21], when comparing 2-dimensional RVE and 3-dimensional UCFEM models, the circle cell periodic pattern cannot reproduce the response of a highly irregular porous structure. Therefore, in the subsequent paragraphs we focus exclusively on the RVE modeling approach to reproduce porous-like behavior, using the fully dense material characteristics determined from the Gibson-Ashby power law scaling relation.
Determining the SMA constitutive model parameters by using an RVE approach is not a straightforward procedure. Indeed, only and can be easily determined with a stress-strain curve, as they are the intersections of, respectively, an Austenite elastic line with an phase transformation line and an line with an Austenite line. The determination of and and is quite problematic, because in a porous sample, contrarily to a bulk sample, Martensitic transformation is never entirely completed at the end of a test. To support this point, we present in Figure 11 (Figures 11(a)–11(c)) the average volume fraction of each phase (, , or ) at the end of test (). For PVF = 0.3, only 30% of the material is fully transformed to Martensite, whereas for PVF = 0.7, this quantity is much lower (1.7%).
(a) PVF = 0.3
(b) PVF = 0.5
(c) PVF = 0.7
These results are in agreement with Panico and Brinson [22], who showed that a full Martensitic material cannot be reached for a PVF of 0.4 and greater, illustrating the difficulty of accurately defining , , and . To circumvent this obstacle while complying with Ansys constitutive model, we will evaluate the evolution of the four slopes of the bilinear SMA model (denoted by , , , and , resp., for Austenite, , Martensite, and phase states; see Figure 1) as a function of the relative density. Moreover, because of the irregularity of the porous sample, the initial sharply defined behavior exhibited by the fully dense material (Figure 1) is smoothed using transitions between the four different phase states (, , , and ).
For the purpose of illustration, Figure 12 depicts the equivalent nodal stress state reached for a global strain of for PVF = 0.3, 0.5, and 0.7. The maximum of the scale is set at MPa; that is, all the elements having an equivalent stress greater than 140 MPa are plotted in pitch red and are assumed to be fully transformed to Martensite.
(a) PVF = 0.3
(b) PVF = 0.5
(c) PVF = 0.7
Figure 13 illustrates the evolution as a function of PVF of, respectively, the moduli (, , , and ) and the critical stresses (, ). Just as we did in the previous section for the elastic modulus, the results are fitted by a power law defined in (1) (see Table 3 for values and ).
(a)
(b)
The remaining parameters of the SMA constitutive model implemented in Ansys are then defined using the moduli and the previously defined critical stresses’ scaling relations. It is then possible to model the behavior of the porous material using the fully dense model with the Gibson-Ashby power law (Table 3). Figures 14(a)–14(c) depict the stress-strain diagrams obtained with such an approach (dotted lines), which show an excellent agreement with the mesoscale RVE stress-strain plots (solid lines) for all the PVF values considered.
(a) PVF = 0.3
(b) PVF = 0.5
(c) PVF = 0.7
The authors would like to emphasize the drastic decrease in computation time offered by the scaling relation: in the previous figure, the fully dense model is defined with only 7600 elements versus 144000, 240000, and 336000 elements for mesoscale models with a PVF of, respectively, 0.7, 0.5, and 0.3. For mesh-based numerical methods such as FEM, the complexity of the algorithm grows exponentially with the number of the dimensions and then with the number of elements. In the present case, such a reduction in the number of elements allows a reduction of the computational time from typically 8–20 hours (depending on the convergence and on the PVF value) to less than 5 minutes on a standard desktop computer.
3.3. Model with Gradient of Porosity
In many application areas, to fulfill specific functions it is necessary to create a heterogeneous porosity over all of a structure. Bioinspired materials that try to mimic the mechanical behavior of bones are an excellent example: as the bone structure changes from a dense external structure to a porous internal one, the manufactured part must exhibit such a porosity gradient [23]. To evaluate the influence of the porosity gradation on the mechanical response of foams, two mesoscale (RVE) models with homogeneous (Figure 15(a)) and heterogeneous (Figure 15(b)) porosities presenting the same mean PVF are set up. For the sake of simplicity, the porosity gradation is achieved here using a pile-up of homogeneous layers of different porosity connected together to form a gradated porosity model. Note that for the homogeneous mesoscale model, porosity layers are directly set during the geometry generation, resulting in a smooth connection between them (no visible delimitation). The heterogeneous mesoscale model presents five layers, with a PVF ranging from 0.3 to 0.7 as indicated by the tinting in Figure 15(b). Each layer exhibits 0.1 PVF difference with its neighbours.
(a) Homogeneous mesoscale model
(b) Heterogeneous mesoscale model
Both models are made of SMA (Table 2) and subjected to an identical tensile test up to a total strain of . Figure 16 depicts the stress-strain diagrams obtained. As we can observe in Figure 16, the gradation of porosity has a pronounced effect on the overall response: the homogeneous model presents a stiffer behavior and it can reach a maximum stress that is twice that of the heterogeneous model.
To evaluate the ability of the Gibson-Ashby approach to reproduce the gradation of porosity, three heterogeneous porous models with a mean PVF of 0.5, ranging from 0.55 (at the bottom) to 0.45 (at the top), are set. The difference between these models lies in the number of homogeneous layers used to model the porosity gradation: (i)10 homogeneous layers, referred to as case 1 (PVF difference between adjacent layers: 0.011); (ii)5 homogeneous layers, referred to as case 2 (PVF difference between adjacent layers: 0.025); (iii)2 homogeneous layers, referred to as case 3 (PVF difference between adjacent layers: 0.1).
Figure 17 depicts the biporous case 3, that is, porosity modelled using two homogeneous layers with different PVF values.
(a) Heterogeneous mesoscale model
(b) Fully dense equivalent model
We observe in Figure 18 that the agreement between the fully dense model and the mesoscale one depends on the PVF difference between two adjacent layers. The stress-strain fully dense response is almost superimposed with the mesoscale one until a macroscale strain of 0.006, 0.0055, and 0.0045, respectively, for cases 1, 2, and 3. For higher stains, the fully dense responses seem to diverge quickly from the mesoscale response, but the smaller the difference between adjacent layers, the smaller the divergence. As for the hysteretic behavior of the response, the three fully dense models match the mesoscale responses in the unloading stage in the same way that they did for the loading one: almost superimposed curves for small strains, depending on the PVF difference between the models’ homogeneous layers.
(a) Case #1
(b) Case #2
(c) Case #3
It can therefore be concluded that, as previously stated, it is possible to accurately model the porosity gradation with a fully dense model. The accuracy of the results obtained depends both on the numbers of homogeneous layers and on the maximum macroscale strain achieved; the PVF difference between adjacent layers must be reduced as we go to higher macroscale strain levels to achieve good agreement with mesoscale results.
4. Conclusion
The use of porous forms of SMA is growing in various areas due to the attractive features exhibited. However, designing such a material is still a challenging process as the properties of the macrostructure depend nonlinearly on the microstructure (shape, distribution, etc.). In this paper we set up scaling relations for the superelastic behavior of two-dimensional SMAs as a function of the global porosity of a structure. Even though the homogeneous behavior law thus defined does not account for microscale effects (cracks, local plasticity, etc.), it does offer the possibility to predict SMA foam behavior with a drastically reduced numerical cost and a good level of accuracy. Furthermore, the scaling relations allow SMA foams to be modelled accurately, presenting a gradation of porosity. Finally, as in Shen and Brinson [24], 2-dimensional results cannot be extended to 3-dimensional SMA structures; this topic will be addressed in a future work.