A Computer Simulation of the Nitinol Thermal Expansion under Fast Varying Working Conditions

We discuss the setup of a simulation on ANSYS to predict the thermal expansion of parts made of Nitinol. A simulation is justified for working conditions in which the part heating is not homogeneous originating a temperature gradient across the part section such that an analytical estimate for the part expansion cannot be calculated. We apply the simulation to the topological optimization of a square section geometry and a bullet geometry for water assisted injection molding. For the topological optimization we consider as parameter the wall thickness and consider both the cases of fast varying temperature and fast varying temperature and pressure.


Introduction
Nitinol is a Nickel Titanium Alloy (NiTi) that was originally manufactured at the Naval Ordnance Laboratory [1]. It exhibits both shape memory effect and superelasticity as well as being very resistance to corrosion and fatigue. Due to the molecular organization of the metal alloy it has two distinct phases, Martensitic phase (M ) at lower temperatures and Austenitic phase (A) at higher temperatures. The phase transition M → A (rising temperature) takes place in between the start temperature A s and finish temperature A f > A s and the phase transition A → M (lowering temperature) takes place between the start temperature M s and the finish temperature M f < M s such that A s > M f and A f > M s . The stress-free Martensitic phase (temperature dependent) is characterized by a twinned multi-variant crystallography structure, in addition there is a stress induced Martensitic phase (S, stress dependent) characterized by a detwinned configuration with a single variant crystallographic structure.
The shape memory effect consists in deforming the material below the temperature M s (M phase), when it is heated above the temperature A s (A phase) it recovers its original (recorded) shape. As for superelasticity allows for about 15 times the elastic motion of steel and consists in that, when the material is above the temperature A s , under mechanical loads the stress induces a phase transition A → S, once the load is removed the phase transition S → A also restores its original (recorded) phase. Hence both for the thermal and stress induced phase transitions the material response shows a hysteresis cycle.
Although originally the study of NiTi was carried mostly for the 1:1 Nickel weight percentage of 55%wt corresponding to the atomic percentage of 49.9%a, soon it was realized that the specific values of the transition temperatures and mechanical properties highly depend on the Nickel percentage [2] (see references therein), specifically for NiTi alloys the Nickel atomic percentage ranges from 49.0%a to 50.7%a corresponding to M s between +100 o C to −50 o C (respectively) and the ductility is reduced as Nickel percentage is increased. As for the addition of other elements to the alloy, Platinum and Palladium allow to change M s by −10 o C to +250 o C, Niobium allows to increased the hysteresis width and Copper allows to reduce the hysteresis width. Also both the transition temperatures and mechanical properties of the alloys are highly dependent on the manufacturing and treatment of the parts [3,4,5,6], in particular it is claimed that for lower Nickel percentages the cold working is the predominant factor in setting the material properties while for higher Nickel percentages the heating treatment is the predominant factor. Although industry standards exist [7], it is still hard to fully compare existing results and discussions in the literature as the specific alloy percentages, part manufacturing process and treatment are often omitted. In this work our main objective is to analyze how NiTi part dimensions vary with temperature. In particular, for a given NiTi part, we study how the dimensions at ambient temperature are related with the dimensions at working temperature and pressure. We note that for relatively long operation times, homogeneous temperature and pressure can be assume. Hence the steady state approximation is valid and it is possible to analytically compute the part dimensions at working conditions. However for fast varying working temperatures and pressures we have to account for thermal conduction and non homogeneous expansions. For that purpose we are going to setup a simulation in ANSYS. Reciprocally the simulation also allows to relate dimensions at manufacturing temperature to dimensions at ambient temperature, not exclusively but also including, digital directly manufactured parts.

Simulation Setup
In addition to the NiTi properties dependence on the alloy composition and part manufacturing and treatment, we note that generally the material properties are distinct in both Martensitic and Austenitic phases. Of relevance to the present study we note that the Coefficient of Thermal Expansion (CTE) has been measured experimentally for distinct samples with distinct treatments [3,5] and generally varies in between two distinct stable values for both NiTi phases. Moreover the CTE during the phase transition is generally non-linear and as expected due to the shape memory effect gives rise to a thermal hysteresis cycle [3].
The simulation of the shape memory effect and superelasticity is today available in commercial softwares, for the particular case of simulation using the software ANSYS see for instance [8]. To setup a simulation with specific NiTi material properties including the shape memory effect let us consider the following phase transition temperatures There is already a pre-defined material for shape memory effect, for which we consider the parameters listed in table 1 with reference temperature matching M f . The specific model employed is defined an discussed in ANSYS documentation [9]. The plasticity and material dilatation for the shape memory effect are modeled based on the stress matrix. The Temperature Scaling Parameter β controls the phase transition through the Maxwell stress describing the dependence on the temperature of the shift of the stress hysteresis cycle. However it does not describe neither the standard plasticity nor the standard thermal expansion of the material which must be supplemented independently of the shape memory effect. For the simulation we consider the material properties listed in table 2 for the Martensitic and Austenitic phases.
To fully reproduce the thermal hysteresis cycle we consider as approximation a linear interpolation between the values of each of the material properties for both temperatures A s and A f for rising temperatures (M → A) and M s and M f for lowering temperatures (A → M ). In particular, for both rising and lowering temperature branches of the hysteresis cycle, the material expansion is modeled by the Isotropic Secant Coefficient of Thermal Expansion (see [9]) with the same Zero-Thermal-Strain Reference Temperature T ref for both directions of the cycle which ensures that below M f and above A f the modeled material expansion is exactly matched independently of the temperature gradient direction. For the simulation we consider the reference temperature T ref to match the environment To test the setup we consider a simple geometry corresponding to a bar of square section 4 mm × 4 mm (aligned with the yz-plane) and length of l 0 = 10 mm (aligned with the x-axis). The geometry dimensions are quoted for the environment temperature. We simulate only 1/4 of the bar such that the point at (0, 0, 0) mm is fixed. For the evaluation of displacements and strain we consider the point along the x bar symmetry axis at x = l = 10 mm, (0, 0, 10) mm. For the simulation setup see figure 1.  We recall that the coefficient of thermal expansion is defined with respect to the original length at environment temperature T ref (it corresponds to the secant for each temperature), such that for experimentally measured values at each of the NiTi phases for the coefficient

Direct Digital Manufacturing and Polymers
the material expansion can be analytically computed ∆l = l 0 α (T − T 0 ). For a given part under a known pressure and thermal distribution, both to test the simulation and evaluate the effective value of α at each time step i we obtain the discretized values of the coefficient of thermal expansion to be To test the simulation we consider that the bar is subject to a homogeneous pressure of 0.  As for the expansion along the x-axis and numerical CTE (4) as functions of the temperature are plotted in figure 3 representing the thermal hysteresis cycle.
The shape memory effect does not play a direct role in the present study, however let us show that it is also properly modeled for the setup just described. Considering the same geometry representing  figure 4. As for the stress hysteresis cycle is plotted in figure 5 for the stress and deviatory stress component xx.

166
Direct Digital Manufacturing and Polymers

Fast Varying Working Conditions
In the previous section we have considered instantaneous homogeneous heating of the full volume. In real working conditions often the Nitinol parts are subject to fast variation of temperature and pressure. As it takes some time for the bulk temperature to rise to the surface temperature, there is a temperature gradient across the section and the part expansion only matches the surface temperature after some time. This is due to the heat conduction across the part section (from the surface boundary to the bulk) not being instantaneous, it is ruled by the standard heat conduction equations [9]. As an example let us consider the same geometry and a surface temperature profile rising from T (t = 0.1s) = 22 o C to T (t ≥ 1s) = 150 o C for the 2 external faces of constant y and z such that the simulated geometry is interpreted as a section of a larger part (the remaining faces are, for simulation purposes considered perfectly insulated). The time evolution of the maximum and minimum temperature obtained are plotted in figure 6, as well as a temperature contour at time t = 1s. We note that only at time t ≈ 2.5s the part is fully homogeneously heated to the temperature T = 150 o C. Also, only at t ≈ 2.5s the part expansion will match the expected value for T = 150 o C while, for instance at t = 1s, will expand only up to 80% of that value. We present in figure 7 the expansion along the 3 axis at t = 1s and the plot of the external faces expansion (y, z-axis deformation) and of the length expansion (x-axis deformation) as a function of time.
We note that, as the temperature rise is not homogeneous, both the Austenitic phase and Martensitic phase will generally be present while the part heats such that the coefficient of thermal expansion is varying across the part section. In addition when superficial treatments for the part are being considered the NiTi material properties may also be varying across the section (even for homogeneous heating) such that it is harder to carry an analytical computation for the part expansion.
For larger parts than our simple test geometry the heating time is larger and the temperature gradient becomes more significant. Our simulation is useful for parts intended to work on relatively high temperature environments for short time periods (without preheating of the part) and for which the exact dimensions at working conditions are known. In the next section, as an example, we consider a topological optimization for a simple hollow part also of square section.

Topological Optimization -An Example
As an example of applicability of the simulation presented in this work let us consider a part of square section with dimensions 10 mm × 10 mm and undetermined length. For the simulation we are considering a length of 10 mm and 1/4 of the part with the symmetry axis aligned with the xaxis. The part is intended to work for about t = 1.2s under a temperature profile that rises from environment temperature T (t = 0.1s) = 22 o C to T (t = 1s) = 150 o C and a pressure profile that rises from P (t = 0.05s) = 0.1 MP a to P (t = 0.2s) = P max . We note that the homogeneous heating of the part to the work temperature would correspond to a expansion of ∆ x = ∆ y = 0.00601600 mm with respect to the part symmetry axis (x-axis), however for the specified working conditions this value will never be reached as the material in the bulk of the part will not have time to reach the external surface temperature.
We are further assuming that the part can be hollow with the inner cavity filled with air which we consider at t = 0s to be at a pressure of 0.1 MP a and environment temperature of 22 o C. We will aim at optimizing both the part sectional deformation (y, z-axis) and the material amount (volume) required to manufacture the part. For this purpose the only parameter of the topological optimization will be the wall thickness of the geometry l w ∈ [0.5, 5] mm where l w = 5 mm corresponds to a solid part and for analysis parameters we are considering the y, z-axis expansion of the part due to both temperature and pressure at the time t = 1.2s. We consider two distinct cases:

168
Direct Digital Manufacturing and Polymers (i) P max = 0.1 MP a demanding a maximum working transversal deformation at t = 1.2s of |∆ y | = |∆ z | < 0.005 mm which corresponds simply to an optimization on the thermal expansion.
For the temperature profile simulation we are considering a setup with the external surfaces heated to the working temperature and the remaining surfaces considered to be perfect insulated (representing the part symmetry and sectioning). As for the structural analysis we are considering a setup with the inner surfaces subject to a pressure of 0.1 MP a and the outer surfaces subject to the pressure profile P (t). Both these setups are represented in figure 8 for a solid and a hollow part. An example of the temperature profile for the part at t = 1.2s is plotted in figure 9 and for the part deformation along the z-axis is plotted in figure 10.
Running several thermal and structural simulations for distinct values of the wall thickness t w we obtain both the part thermal expansion and deformation under pressure at time t = 1.2s as a function of this parameter. We plot the results obtained in figure 11.
The topological optimization process should be iterative such that for the points for which the optimization parameters are near the desired parameter value (aim) further optimization points must be generated (using for example the bisection method). For the cases we are analyzing we obtain for case (i) the value of t w = 1.875 mm corresponding to a deformation of |∆ x (t = 12s)| = |∆ y (t = 12s)| = 0.00492400 ± 0.00000492 mm and for the case (ii) the value t w = 2.375 mm corresponding to a deformation of |∆ x (t = 12s)| = |∆ y (t = 12s)| = 0.0175020 ± 0.0001750 mm. We note that we have use as only parameter the part deformation strictly at the time t = 1.2s, as already mentioned the value of the deformation is not uniform over time and depending on the application the topological optimization may include as additional parameters the deformation at several distinct times. Also generally and depending on the applications additional analysis are often considered such as the part fatigue under the working cycle of the part.

Direct Digital Manufacturing and Polymers
Applicability to Direct Digital Manufacturing -Injection Bullet As a possible example of applicability of the simulation and topological optimization developed here let us consider one of the bullet geometries discussed in [10] and intended to water assisted injection molding. These bullets aim at hollowing the plastic parts at the molding phase of compactification and the process consists on injecting water in the back of the bullet such that it is propelled forward inside the molded plastic. Often these bullets are, for research purposes, directly manufactured by extrusion molding as described in [10]. For this specific simulation we are considering a bullet made of Nitinol with cylindrical symmetry around the x-axis, maximum radius of 9 mm and length of 23 mm, and for the simulation setup we consider only 1/4 of the bullet geometry as represented in figure 12. Let us further consider the same temperature profile of the previous section rising from the temperature T (t = 0.1s) = 22 o C to T (t = 1s) = 150 o C. To simulate the injected water we consider the pressure profile normal to the back (inside) surface of the bullet rising from P in (t = 0.05s) = 0.1MP a to P in (t = 0.2s) = 25MP a and to simulate both the plastic pressure and dynamic friction we consider (respectively) a pressure profile normal to the front (outside) surface of the bullet rising from P out (t = 0.05s) = 0.1MP a to P out (t = 0.2s) = 16MP a and a directional pressure profile parallel to the bullet symmetry axis applied to the front (exterior) surface of the bullet rising from P out(dyn) (t = 0.05s) = 0.0MP a to P out(dyn) (t = 0.2s) = 5MP a (with direction from the front to the back of the bullet). For topological optimization of the Nitinol bullet we are considering as only parameter the wall thickness t w ∈ [0.25, 2.5] mm, as in the previous section we are aiming at optimizing the volume maintaining the radial deformation at the time t = 1.2s below a critical threshold. Specifically we consider the two distinct cases: (i) P in = P out = 0.1 MP a and P out(dyn) = 0 MP a demanding a maximum radial deformation over the time period t ∈ [0.2, 1.2]s of |∆ y | = |∆ z | < 0.01 mm which corresponds simply to an optimization on the thermal expansion.
As an example of the simulations results we present in figure 13 the deformation for the bullet geometry with wall thickness t w = 1 mm.
Running the simulations for several wall thickness and employing the bisection method we obtain for the topological optimization (i) t w = 1.00 mm corresponding to a maximum radial deformation of |∆ y | = |∆ z | = 0.00984510 ± 0.00000985 mm and for the topological optimization (ii) t w = 0.563 mm corresponding to a maximum radial deformation of |∆ y | = |∆ z | = 0.0197790 ± 0.0001978 mm. We plot the results obtained for both optimization cases in figure 14.

Conclusions
In this study we have setup a simulation in ANSYS to model the Nitinol thermal expansion reproducing the thermal hysteresis cycle. We have also discussed possible applications for parts modeling under fast varying working conditions, in particular we have presented an explicit topological optimization example for both a square section geometry and a bullet geometry in the context of directly manufactured parts for water assisted injection molding. We recall that as discussed in the introduction both the alloy percentages, the manufacturing process and the treatment affect the Nitinol material properties for both phases, hence either the material to be simulated should already have known properties or the properties must be measured before setting up a simulation that matches the real part.

172
Direct Digital Manufacturing and Polymers