International Journal of Plasticity 16 (2000) 541±562 www.elsevier.com/locate/ijplas
Simulations of localized thermo-mechanical behavior in a NiTi shape memory alloy John A. Shaw * Department of Aerospace Engineering, The University of Michigan, MI, USA Received in final revised form 5 October 1999
Abstract Previous experiments have shown that stress-induced martensitic transformation in certain polycrystalline NiTi shape memory alloys can lead to strain localization and propagation phenomena when loaded in uniaxial tension. The number of nucleation events and kinetics of transformation fronts were found to be sensitive to the nature of the ambient media and imposed loading rate due to the release/absorption of latent heat and the material's inherent temperature sensitivity of the transformation stress. A special plasticity-based constitutive model used within a 3-D ®nite element framework has previously been shown to capture the isothermal, purely mechanical front features seen in experiments of thin uniaxial NiTi strips. This paper extends the approach to include the thermo-mechanical coupling of the material with its environment. The simulations successfully capture the nucleation and evolution of fronts and the corresponding temperature ®elds seen during the experiments. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: Localization; Phase transformation; Shape memory alloy; Constitutive behavior; Finite elements
1. Introduction There are a number of remarkable materials, such as shape memory alloys (SMAs), piezoelectric materials, magnetostrictive materials, etc., which exhibit strong coupling between their mechanical behavior and other ®elds, such as thermal, electric or magnetic ®elds, respectively. These materials are often referred to as ``interactive'', or (to stretch the point) ``smart'' materials, in that they sense a change in their environment and respond in a mechanical way. The unique properties of these materials are not new, having been discovered 30±60 years ago, yet their use for * Fax:+1-734-763-0578. E-mail address:
[email protected] (J.A. Shaw). 0749-6419/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S0749-6419(99)00075-3
542
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
application in active structures is relatively new and recently seems to be gaining momentum. SMAs, such as NiTi (or Nitinol), exhibit two remarkable properties, the shape memory eect and pseudoelasticity (see Fig. 1). The shape memory eect ( to in Fig. 1) is the material's ability to recover large mechanically-induced strains (up to 8%) by moderate increases in temperature (10±20 C). Pseudoelasticity ( to ) refers to the ability of the material in a somewhat higher temperature regime to accommodate strains of this magnitude during loading and then recover upon unloading (via a hysteresis loop). The underlying mechanism is a reversible martensitic transformation between solid-state phases, often occurring near room temperature. The transformation can be induced by changes in temperature or by changes in stress causing a strong thermo-mechanical coupling in the material behavior. The material also has very nonlinear mechanical behavior, high internal damping, and high yield stresses. All of these properties make NiTi a promising candidate for novel structural applications (see Perkins, 1975; Funakubo, 1987; Otsuka and Wayman, 1988; Duerig et al, 1990). NiTi's remarkable behavior arises from the interplay of two phases (see lattice sketches in Fig. 2), a high temperature phase (austenite), having a cubic lattice structure, and a low temperature phase (martensite), having a monoclinic structure (Otsuka et al., 1971). Due to its low degree of symmetry, the martensite phase exists either as a randomly twinned structure (low temperature, low stress state) or a stress-induced detwinned structure that can accommodate relatively large strains without permanent deformation.
Fig. 1. Thermo-mechanical response of NiTi wire in water: shape memory eect response ! (experiment taken from Shaw, 1997).
! ; pseudoelastic
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
543
Fig. 2. Dierential scanning calorimetry of NiTi strip.
Despite the fact that the material was discovered nearly 40 years ago (Buehler et al., 1963), constitutive models have developed rather slowly, hampered by the complexity of the material behavior and the somewhat limited experimental basis for many years. The materials science literature is rich on the subject, and the understanding of the micromechanical aspects has reached a mature level. However, the analytic bridge between microscopic and macroscopic behavior is quite complex and is a recent area of research (see, for example, Ball and James, 1987; Batthacharya and Kohn, 1996; Siredey et al., 1999). The last decade has seen some notable constitutive modeling eorts (see Tanaka et al., 1986; Brinson, 1993; Levitas, 1994; Boyd and Lagoudas 1994; Patoor et al., 1995) and some new experimental ®ndings (see Leo et al., 1993; Shaw and Kyriakides, 1995; Sittner et al., 1995; Gall et al., 1999), providing new impetus for design and application. Yet, it is fair to say that reliable constitutive models suitable for many engineering applications are not yet available, especially under cyclic loading conditions. In particular, few of the SMA constitutive models acknowledge the material instabilities which have been observed in pseudoelastic NiTi (Shaw and Kyriakides, 1997). [For notable exceptions see James (1983) and the 1D thermodynamic framework of Abeyaratne and Knowles (1993) and Knowles (1999)]. A continuum-level plasticity approach with a special trilinear eective stress±strain response was recently shown to capture many features of the localized deformation ®elds seen during unstable stress-induced transformation in uniaxially loaded thin strip NiTi (Shaw and Kyriakides, 1998). The approach was limited to quasi-static, isothermal, irreversible behavior, representative of the material response for very slow loading rates in a convective medium, such as water. In this paper the approach will be
544
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
extended to include the thermal interactions between a mechanically loaded NiTi specimen and its environment to capture the nonisothermal response at higher loading rates and in a less convective medium, such as air. Some experimental observations from Shaw and Kyriakides (1997) are ®rst reviewed and then the modeling approach is demonstrated through some ®nite element simulations. 2. Experimental observations This paper focuses on the material instabilities which occur in pseudoelastic NiTi loaded in uniaxial tension. Experimental studies of the phenomena of initiation and kinetics of phase transformation fronts and the associated thermal sensitivities are discussed in detail in Shaw and Kyriakides (1995, 1997) and Shaw (1997). This section brie¯y reviews some general experimental ®ndings and discusses in particular two experiments that will be simulated in the following section. First, it was shown experimentally that the ``rate eect'' in the material is not the usual viscoelastic phenomena, but rather, a strong interaction between the latent heat of transformation and the material's extreme (rate-independent) sensitivity to temperature. A dierential scanning calorimetry (DSC) response is given in Fig. 2 for the polycrystalline, near-equiatomic NiTi alloy used in Shaw and Kyriakides (1997). The area under each power peak/valley separating the three solid-state phases, austenite (A), martensite (M), and rhombohedral-phase1(R), represents the latent heat of transformation. A set of displacement controlled, isothermal mechanical responses at various temperatures is provided in Fig. 3a. Notice the nucleation peaks that must be surmounted at the beginning of the stress plateaus (during loading). Fig. 3b shows how the nucleation and transformation stresses (stress plateaus) increase signi®cantly with increasing ambient temperature. The exothermic A!M transition during loading2 tends to cause self-heating, which in turn raises the material's underlying transformation stress according to Fig. 3b. Consequently, the nature of the ambient media, gas versus liquid, plays a surprising role in the apparent material response due to the dierences in the prevailing heat transfer environment (primarily the convective properties). This issue is an unusual, and perhaps counter-intuitive eect among typical structural materials. For example, a distinct stress plateau (unstable behavior) is observed during transformations if the specimen is kept isothermal, i.e. slow loading rate and a convective medium like water. A seemingly stabilized stress±strain behavior (positive tangent modulus) is observed under adiabatic conditions, i.e. faster loading rate or a relatively nonconvective medium, like air. It should be noted that this behavior is opposite to the destabilizing eects of temperature rise during adiabatic shear banding in more conventional materials (Molinari and Clifton, 1987), owing to the dierent microstructural deformation mechanisms. 1
Note, the R-phase is not encountered in the remainder of this paper. The endothermic transition (M!A) during unloading tends to cause self-cooling, but this transformation is not the subject of this paper. 2
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
545
Since transitions in untrained NiTi can occur in a mechanically unstable manner, the transformation may occur in a localized way, i.e. through nucleation events and subsequent propagation of distinct phase fronts along the length of a uniaxially loaded specimen. Fig. 4 shows an experiment taken from Shaw and Kyriakides (1997) in which a NiTi thin strip specimen (gage section: t=0.4, w=2.5, L=39mm, and free length: Lf=50.8mm) was stretched at a slow elongation rate : (f =L 10ÿ4 sÿ1 ) in room temperature air (25 C). The specimen was initially austenite, and unstable transformation to martensite began once a critical stress level was reached. During this transformation, photographs were taken at 40 s time intervals (see Fig. 4a). The large changes in strain caused an observable color change of the specimen surface due to the disturbance of a naturally occurring brittle oxide layer. In this way the distinctly inhomogeneous evolution of the transformation was
: Fig. 3. (a) Pseudoelastic responses of NiTi strip at f =L 10ÿ4 sÿ1 ; (b) ®t of nucleation and propagation stresses for A!M transformation (loading).
546
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
: Fig. 4. Experiment 1. Evolution of A!M transformation (loading) at f =L 10ÿ4 sÿ1 in room air: (a) photographs, (b) infrared thermal images, [experiment taken from Shaw and Kyriakides (1997)].
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
547
tracked optically. Simultaneously, an infrared thermal imaging system provided synchronized temperature of the specimen (see Fig. 4b where the color temperature legend spans 5 C centered around the room temperature). The motion of transformation fronts in Fig. 4a started with a nucleation event3 near the top of the gage section where there was a slight stress concentration at the taper. A single transformation front, seen as a shear band, moved down the specimen axis as the transformation progressed from time to time . The exothermic transition is seen in the thermal sequence in Fig. 4b as a yellow spot, then small red spot, moving down the length of the specimen. Once the temperature of the front rose by about 3.5 C a second nucleation occurred near the lower end of the gage section at time . At this time there were two converging fronts, they shared the overall transformation rate (prescribed by the end displacement rate) and travelled half as fast as when there was a single front. Consequently, a smaller temperature peak existed (green spots) until time . As the fronts neared each another, however, they began to thermally interact, and a temperature rise was observed between time and time . Once the specimen gage section was completely transformed at time the specimen temperature returned to ambient. Fig. 5: shows the same type of experiment conducted at a loading rate 10 times faster (f =L 10ÿ3 sÿ1 ). In this case there were four nucleation events, two staring at the ends of the gage length at time and time and two occurring in the interior before time and time , leaving as many as six simultaneously traveling fronts in the specimen. The thermal sequence showed severe self-heating throughout the experiment (note the larger 20 C range in the thermal legend). This allowed multiple nucleation barriers to be overcome. In fact, front temperatures exceeded 37 C during latter parts of the transformation, due to the higher rate of transformation and the lack of time for convection to occur. During these temperature ¯uctuations, the stress changed according to the propagating stress±temperature trends of Fig. 3b, interrupted by small stress dips whenever new fronts were nucleated or when two fronts coalesced. Equilibrium required the average axial stress to be constant along the length. Consequently, the nucleation stress of a colder region could be surmounted if the local temperature rise of a propagating front caused its propagation stress level to rise by at least the size of the nucleation peak. Therefore, more fronts were observed at higher loading rates where local self-heating was more severe. To summarize, the number of propagating fronts was determined by the number of nucleation events that could occur. This, in turn, was determined by the amount of excess self-heating and nonuniform temperature ®elds necessary to overcome the mechanical nucleation barriers. In addition, at these low to moderate strain rates the speed: of all fronts :were nearly the same according to the kinematical relation c =
n, where is the elongation rate, n is the number of currently traveling 3 Note: nucleation in this context refers to nucleation of transformation fronts, not nucleation of martensite. It has been correctly identi®ed that some A!M transformation may procede and follow the stress plateau (Liu et al., 1998). For example, one can see the homogenous self-heating due to some early A!M transformation at time in Fig. 4b before the start of localized deformation.
548
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
: Fig. 5. Experiment 2. Evolution of A!M transformation (loading) at f =L 10ÿ3 sÿ1 in room air: (a) photographs, (b) infrared thermal images [experiment taken from Shaw and Kyriakides (1997)].
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
549
fronts, and is the strain jump along the plateau. See, for example, the change in front speed after in Fig. 4a when the second nucleation event occurred. Furthermore, the kinematics of such an inhomogeneous deformation has a signi®cant eect on the aforementioned temperature sensitivity of the material, since the actively transforming region is local to the traveling front rather than uniformly distributed along the length.: This induces a hypersensitive ``rate'' eect, occurring at global strain rates, such as f =L 10ÿ3 sÿ1 , that one might normally consider quasistatic. 3. Finite element simulations The ®nite element approach of Shaw and Kyriakides (1998) was a ®rst attempt to capture the 3-dimensional deformation ®elds during the nucleation and propagation in thin strips under uniaxial tension. It used a conventional J2 plasticity model that was calibrated to a special trilinear stress±strain curve. In an isothermal, purely mechanical, setting it successfully captured details of the transformation front features that have been observed in experiments both on NiTi (martensitic transformation) and ®ne grained mild steel (Luders bands). Since it was a plasticity-based model, it ignored the reversibility of the deformation in NiTi upon unloading. This limitation will also be accepted here by focusing only on the A!M transformation. The approach will be extended, however, to include the thermo-mechanical coupling with the ambient environment. This introduces a time scale into the calculations (that was absent before) which is needed to predict the number of fronts and the stress history for experiments conducted at typical loading rates in air. The numerical approach will ®rst be outlined and then comparisons will be made between the experiments already presented and two corresponding numerical simulations. 3.1. Numerical approach The existence of distinctly inhomogeneous deformation ®elds has important implications on modeling. The global force±displacement (engineering stress±strain) response must lose (momentarily) its positive slope, since during propagation, multiple strain states are possible for a single axial stress state. Several investigators have studied the behavior of a 1-D continuum solid with a nonconvex strain energy density function (see Erikson, 1975; Falk, 1980; Ericksen, 1991; Abeyaratne and Knowles, 1993). They often treat the ensuing phase front as a mathematical discontinuity in strain and develop appropriate jump conditions for the moving boundary. Although the problems simulated here are uniaxial in nature (actually generalized plane stress), the response is investigated in a 3-D continuum setting. This precludes abrupt strain discontinuities (the front is a propagating ``neck'' with a ®nite width and pro®le shape) and allows one to follow the evolving deformation ®elds without any assumptions on front morphology. The numerical simulations that will be presented were performed with the ABAQUS ®nite element software (HKS, 1997) using a fully 3-D, transient, coupled
550
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
displacement temperature analysis. From a mechanical point of view, it is a quasistatic analysis since inertial eects are neglected, yet it is a transient calculation from a heat transfer standpoint. The constitutive model for the material is based on traditional rate-independent, ®nitely deforming J2 ¯ow theory with isotropic hardening (available as a built-in material type in ABAQUS). In this way the material is assumed to be homogeneous, isotropic, plastically incompressible, and plasticity irreversible. The homogeneity assumption is reasonable, since NiTi is ®ne grained (grain size typically 5m). The isotropic assumption is reasonable if there is no signi®cant crystallographic texture. It is recognized this may not be the case for heavily drawn or rolled NiTi (see Gall et al., 1999). Since the stress states calculated here are largely uniaxial tension (with only a small bending and shear component), the assumption does not seem to aect the results adversely. The irreversibility of the deformation does not aect the validity of the results provided that unloading does not occur: below the critical stress for the reverse transformation. Only monotonic loading, > 0, will be analyzed. The key ingredient in the constitutive model is a trilinear uniaxial true stress±log strain curve having an up±down±up shape. This is used to calibrate the plasticity ¯ow rule. Using the ®t of the measured nucleation stresses and propagation stresses of Fig. 3b, a series of shifted/scaled trilinear curves is constructed to model the temperature dependence of the isothermal material behavior. The derived temperature dependent stress±strain model is shown in Fig. 6. Each trilinear curve is
Fig. 6. Constitutive model: nominal stress±strain and true stress±log strain.
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
551
constructed to model the nucleation stress, N , propagation stress, P (i.e. the Maxwell stress of the trilinear model), and transformation strain, (see Table 1). The intermediate branch of the trilinear curve has a negative slope in both the engineering (reference) stress and true (Cauchy) stress representations. This is necessary to simultaneously model the as-measured nucleation peak and transformation strain of the material. (Note, this is distinct from the approach of Hutchinson and Neale (1983), in which the engineering stress±strain response has an unstable branch but the true stress±strain response is everywhere stable.) However, it poses some theoretical and numerical diculties in the form of a loss of ellipticity of the underlying incremental boundary value problem and the potential for a meshdependent ®nite element calculation. In our case, however, this does not adversely aect the results as will be discussed at the end of this section. Table 1 Trilinear ®t of true stress±log strain T ( C)
E1 (GPa)
1n(1+E1)
true;1 (MPa)
1n(1+E2)
true;2 (MPa)
E3 (GPa)
15 25 35 45 55
52.0 55.5 59.1 62.6 66.1
0.006511 0.007621 0.008600 0.009470 0.01025
338.6 423.2 507.9 592.6 677.5
0.05350 0.05399 0.05450 0.05504 0.05559
322.5 389.9 457.6 525.6 593.8
15.0 15.0 15.0 15.0 15.0
The coupled deformation-temperature analysis available in ABAQUS solves simultaneously for equilibrium and the heat equation. It allows one to choose a portion
f of the speci®c inelastic work to be converted to an internal heat source as follows, 1 : :
1 q s f p ; : : where qs is the heat generation rate per unit mass, is the mass density, and p = is the speci®c plastic work rate. This is used here to simulate the release of latent heat during the A!M transformation. The phase transformation, however, generates latent heat through a speci®c enthalpy change
h, which is decomposed here into stress independent and stress dependent parts as 1 A!M A!M ÿ hO ÿ P ;
2 ql ÿh !M is the stress free speci®c where ql is the speci®c latent heat change and hA O enthalpy change. Integrating Eq. (1) along the Maxwell stress, neglecting the small change in elastic energy, and equating with Eq. (2) leads to the temperature dependent factor of Eq. (3).
f1ÿ
A!M hO p
3
552
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
!M This factor is greater than one (since hA is negative) and changes with temO perature according to the temperature dependent propagation stress. Since ABAQUS does not permit values greater than one, a value of unity was assigned and the thermal parameters (speci®c heat, thermal conductivity, and convective ®lm coecient) were scaled down by this factor as a function of temperature. This results in the correct temperature ®eld due to self heating without abandoning the built-in capabilities of ABAQUS. The environment is modeled by assigning a convective ®lm coecient to each free surface of the specimen (a typical value for stagnant air was chosen). Heat conduction and radiation to the environment were judged to be minor heat transfer mechanisms and were neglected. As a ®rst approximation all thermal parameters, speci®c heat, thermal conductivity, ®lm coecient, and thermal expansion coecient, were modeled to be constant and independent of strain or temperature. (It is known that thermal conductivity and expansion coecient, in particular, may depend on whether the material is austenite or martensite). Table 2 shows the chosen thermal and physical parameters. The latent heat was a measured DSC value for the NiTi alloy used. The thermal conductivity, k, speci®c heat, C, thermal expansion coecient, , and mass density, , are typical values taken from vendor literature. A few dierent values of the convective ®lm coecient, h, were tried between 2 and 20 W/m2K typical of free convection in stagnant air (Incropera and DeWitt, 1996), and a value of 4 seemed to give good results. Table 3 provides the ®tted nucleation and propagation stresses and the calculated scale factor f as a function of temperature. The temperature of the specimen ends is ®xed at the ambient temperature (25 C), modeling the grips as perfect heat sinks. This is a reasonable assumption for the rather large metallic grips used. One end of the free length of the specimen is ®xed
Table 2 Thermal and physical properties
Enthalpy change (zero stress) Thermal conductivity Convective ®lm coecient Speci®c heat Thermal expansion coecient Density
Parameter
Value
!M hA O
ÿ12.3 18 4 837 10 6.5
(J/g) k (W/m K) h (W/m2 K) C (J/kg K) (10-6/K) (g/cc)
Table 3 Mechanical parameters T ( C)
N (MPa)
P (MPa)
f
15 25 35 45 55
336.4 420.0 503.5 587.0 670.6
320.6 393.8 467.0 540.2 613.4
0.04990 0.05018 0.05061 0.05115 0.05178
6.00 5.05 4.38 3.89 3.52
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
553
while the other end is displaced by f . Other surfaces are traction free. In addition, only half of the specimen thickness is modeled by enforcing symmetry conditions, zero normal displacement and zero normal heat ¯ux, at the mid-thickness plane. These boundary conditions are speci®ed in Eq. (4). ÿ w
x; y; 0 0 u
0; y; z 0 uÿ Lf ; y; z f ; v
0; 0; z 0 v Lf ; 0; z 0
4 ÿ @T
x; y; 0 0: T
0; y; z Tamb; T Lf ; y; z Tamb; @z A few comments should be made regarding possible mesh sensitivity of the ®nite element calculation when using a true stress±strain model with a negative slope. It has been well established that admitting an unstable branch in a true stress±strain response leads to discontinuous deformation gradients and ®ne phase structures which depend on the numerical mesh (see Tvergaard et al., 1981; Silling, 1988). This diculty can be resolved by introducing a ``penalty'' in the strain energy density in the form of strain gradient terms (see for example, Triantafyllidis and Aifantis, 1986; Triantafyllidis and Bardenhagen, 1993; Hutchinson and Fleck, 1996). This eectively introduces a length scale into the calculation that must be speci®ed a priori. Since the strip thickness is the minimum continuum length scale of interest, it appears that by keeping one ®nite element through the half thickness the same is accomplished here. Any ®ne phase structures are thereby suppressed in the calculations, and this is acceptable, since that level of detail is not of interest. The homogenization of the ®ne scale does not seem to adversely aect the calculations, and the good results below seem to con®rm this view. The only practical issue was occasional convergence problems at the point where the stress±strain model slope suddenly turns negative. This was resolved by introducing a small geometric imperfection to initiate the ``transformation'' and by adjusting the time step increment at certain critical points in the calculation. Furthermore, as discussed in Shaw and Kyriakides (1998) only the local neighborhood of a transformation front where material is actively transforming is aected, since material elsewhere exists on a stable branch of the stress±strain curve. For the isothermal case only minor mesh sensitivity was seen in the details of the neck pro®le when re®ning the in-plane mesh (Shaw, 1997). Adding the temperature coupling here seems to improve the situation further, since self-heating has a stabilizing eect on the material behavior. A mesh re®nement study was conducted on a simpler version of the current problem, and no mesh dependency was observed when increasing the planar density of the mesh by nearly two orders of magnitude. It appears that a very ®ne mesh indeed would be needed to resolve ®ne phase structures within a front pro®le. 3.2. Simulation 1 The ®nite element mesh used in both simulations is shown in Fig. 7. The elements are eight-node brick continuum elements, with linear strain and temperature
554
J.A. Shaw / International Journal of Plasticity 16 (2000) 541±562
interpolation (ABAQUS type C3D8HT). The gage section of the model consists of one element through the half-thickness, 12 elements across the width, and 182 elements along the length. The mesh was chosen to give elements which are nearly cubic in shape to avoid introducing any directional bias into the calculation. A slight dent is located along the side near the top of the gage length to control the location of the ®rst nucleation (shown with the ``