Induction heating of a shape memory alloy beam

The shape memory alloy heating by eddy currents is a quick solution for the shape change. Then, the analysis of the temperature field as a function of the shape is important to build a mechanical model in large deformation. Even if the temperature can be obtained by experiment, a computational model is useful. The computation of the induced currents in a nickel–titanium shape memory alloy beam is here considered with a T V model adapted to thin shells with the help of a change of coordinates. It allows us to take into account the shape change, without the need of remeshing, as a function of the temperature. Experiments are carried out to validate the model.


Introduction
One of the main types of shape memory alloys (SMA) is nickel-titanium (NiTi). It is used in various applications such as actuators [1]. However, the actuation frequency is low. To overcome this limitation, body heating with eddy current is investigated.
The general operating principle of the SMA is as follows. At low temperatures, the alloy is in its martensitic state. A mechanical constraint is applied to distort the piece permanently. If this metal is then heated beyond its transformation temperature, it changes from martensitic to austenitic state. Provided the deformation does not include folds with too low curvature, the alloy comes back to its initial state (Fig. 1).
The main factor of the process is the heating time or cooling time to change phase. Joule heating is a good candidate to decrease the heating time [2]. The induction heating has been evaluated for wires, as well as tubes in torsion and recently considered in the case of beams [3] but only in the case of small deformations (i.e. with linear strain-displacement equations). A complete mechanical model requires semi-empirical data during the phase change: the control of the rise of temperature with respect to the time as well as its spatial homogeneity is important to determine numerically the phase of the alloy.
Here, the aim of the paper is to use a flexible model of computation of induced currents to study the sensitivity of the different parameters (beam geometry, frequency, decentring) for a given shape. The change of shape is included in the case of large deformations (the linear straindisplacement equations are then no longer valid).
To leave the full mechanical model aside, the deformation of the beam is given as a function of time with the help of the well-known elastica model [4]. It allows us to compute the induced currents and the temperature for each given shape as a function of time.

Geometric parameterization of the change of shape
The mechanical assumption is that the parametric equation of the beam (Fig. 1) is the one of the elastica [5]: where ' ∈ [ À c, c]. c is determined by the constraint that the beam length is constant: This parametric equation comes from the elasticity of a beam in large deformation [4]. Briefly, if a cantilever beam (Fig. 2) is considered (i.e. the beam is fixed at one side and a force F n is normally applied at the other side), the deflection equation is as follows: where u is the tangent angle with respect to the horizontal axis, E is the Young's modulus, I the bending moment of inertia, s the curvilinear abscissa, and u f is the value of the angle u at the end side. Then, equation (3) after integration becomes where u max is such that _ uðsÞ ¼ 0 (Fig. 2). Equation (1) can be found with The same deflection is obtained if the fixed end is replaced by a symmetric axis.
It has been experimentally noticed on top of Figure 1 that the distorted beam has globally the shape of an elastica, which corresponds to a relaxation of the elasticity stress tensor after a plastic deformation (see Fig. 2 when F n decreases).
The change of shape during the heating is taken into account by using the elastica parameterization (Eq. (1)): the decrement of parameter k allows continuously following the change of shape, from the initial shape obtained after relaxation to the final plane original form (corresponding to k ! À 1).

Computational model 3.1 Eddy current computation
The beam region is defined as a domain D of the Euclidean space E 3 . To find a position X in this domain, an application F is built such that any X corresponds to x in a reference space D with F(x) = X.
The coordinates (relative to the frame vectors (k x , k y , k z )) in the space E 3 being given in vector form, we define the map L as The arguments (x, y) of L lie in the subset S ¼ ½Àc; c Â ½À1=2; 1=2 of a dimensionless space where c satisfies the length constraint (Eq. (2)) and the transverse length is l y . (X, Z) are the functions (Eq. (1)) where m is a fixed value (0.02) and k is the parameter that allows the change of shape. The image of S by L is labelled S. The normal N to this surface is used to build the application between x = (x, y, z) and X = (X, Y, Z). H is the height of the beam. And the beam region, noted D, is the image Different formulations can be used to compute eddy currents. The chosen one is a T À V formulation, since it allows to treat only the beam (and not the surrounding domain). The relative permeability of Ni-Ti being 1.002, it is approximated to 1. The T À V equations for conducting domain D are as follows: T vanishes outside D. d is the skin depth (d ≃ 1.4 mm at 100 kHz for the Ni-Ti). H s the magnetic induction source field (of the inductor). It is obtained by a Biot-Savart computation. T is solution of the weak form: The boundary integral obtained in the weak form by the integration by parts vanishes due to the condition T Â N = 0.∇⋅T reduces to the jump of T ⋅ N on ∂D, since T is free divergence in D and vanishes outside D. V is defined in the whole space. Since the permeability is m 0 , to avoid considering the outer mesh of the beam, V is computed with the Biot-Savart formula: Equations (11) and (12) are coupled equations. As the geometry is evolutionary, the remeshing is necessary. For simplicity reasons, an alternative Lagrangian approach is chosen, where the computation of T by  equation (11) is done in the space of reference D, at the cost of a more complex expression of the differential operators. The change of coordinates is done with the parameterization (Eq. (7)). In summary, the following table can be displayed: The electric and magnetic potentials ( In the same way, where n is in the reference space (0, 0, 1/H). The Gram matrix g is introduced as The Gram determinant is For a beam of small thickness with respect to the skin depth (here H ∼ d/10), the main component of T is directed along the normal N, and its variation with respect to N is low. Then, the electric potential vector t can reasonably be approximated by tðxÞ ≃ tðx; yÞn: With the transforms, equation (11) becomes where the source term h s ⋅ n is n(∇ F(x)) À1 ⋅ H s i.e. (H s ⋅ N). Equation (12) is computed by a Biot-Savartlike integral as the unique contribution of t to the computation of v is the jump of t on the upper and lower boundaries of the domain D. w is computed as follows: Equation (19) contains singularities, and truncation errors appear due to the proximity of the upper and lower boundaries. Then, a specific analytic treatment is set up to compute with enough accuracy w and its derivative with respect to z.
To avoid the resolution of a linear system with a nonsparse matrix (the CPU time of a LU decomposition rises as the number of nodes cubed), the computation of t and ∂ z v is alternated: equation (18) is computed for a given ∂ z v, and ∂ z v is computed again with equation (19). It allows to get a convergence for a limited number of iterations. The parallel solver MUMPS (MUltifrontal Massively Parallel Sparse direct Solver) is used to reduce the computation time.

Joule power and thermal model
The source field h s is obtained by Biot-Savart's formula. Joule power losses can be computed after solving the eddy current problem as with its associated resistance R and power density _ q ( t * is the conjugate of t). The electric conductivity s depends on the temperature (1.21 MS/m in the martensitic state and 1.31 MS/m in the austenitic state).
The Fourier number is about 250 for the normal direction to the beam, but 0.25 for the tangential directions, which means that the temperature is almost independent of this normal direction, but that the transient temperature distribution has to be considered in S. If, at time t the temperature is labelled u n and at time t + Dt as u n+1 , the variational form of the heat equation is obtained in S as where (l,C p ) are the thermal conductivity and the specific heat (both weakly temperature-dependent, i.e. about 2% of variation within the considered temperature range). r is the mass density, u ∞ the air temperature, h the total exchange coefficient. This exchange coefficient considers convection and linearized radiation. The latent heat is negligible.
In summary, the flow chart of Figure 3 displays the thermal-electromagnetic coupling. Due to the low parameter variation as a function of the temperature, the coefficients update is made when the temperature variation is beyond a certain value (chosen at 5°C). The shape change, driven by the given parameter k, occurs when the temperature reaches the transformation temperature (here, 50°C).

Results
The beam dimensions are l X = 80 mm, H = 0.12 mm, and l Y between 12 and 40 mm (which are lower and upper bounds for the inductor diameter). Figure 4 shows the two-wire inductor geometry. The inductor frequency is 200 kHz.

Plane beam
At first, the plane beam case is considered. Figure 5 displays the resistance corresponding to the induced currents as a function of the airgap e (i.e. the vertical distance between the plane beam and the top plane of the inductor). The resistance injected in the beam is measured with a rlc-meter (plain lines with crosses). The measure is carried out at the inductor terminals, and the impedance is zeroed without the beam, to measure effectively the resistance injected in the beam. It is compared to the one computed with equation (20) (dashed lines).
The model accuracy is quite good for a width between 12 and 40 mm.

Elastica beam
Then, the elastica beam is considered (Fig. 6): the parameter k is changed into Dz, the height of the beam (the plane beam k ! À 1 corresponding to Dz = 0). A system to get and maintain the right shape as in Figure 2, but at both ends has been elaborated. The measurement of the resistance corresponding to the power injected in the beam is carried out in the same way as for the plane beam. The noticed uncertainty is about ±1mV. Then, to get a valuable comparison between experiment and measurement, only the 40-mm-wide beam is considered, with various airgaps e (in the range of 1.5À5 mm). There is a good accordance between both, even when Dz is large.    The decentring of the elastica beam y c , with respect to the y-axis (i.e. perpendicular to the length of the beam) is investigated (Fig. 7) at 200 kHz for various airgaps. Due to the uncertainty measurement, the 40-mm-wide beam, as in Figure 6, is considered for a medium value of Dz (11 mm). The measurements have been carried out for three positions: cantered, half-coil radius decentring, and coilradius decentring. Computed and measured resistances are very similar. The decentring of the elastica beam x c , with respect to the x-axis is also studied (Fig. 8). The x c decentring is less sensitive, as the inductor remains under the beam.
The resistance of 40-mm-wide elastica beam (Dz = 11 mm) is computed as a function of the airgap e for various frequencies. The measured values are slightly above the computed ones for frequencies above 200 kHz, which can be explained by both the uncertainty measurements and the fact that H/d increases.

Full deflection model
The beam of width l Y = 12 mm is heated by the two-wire inductor (Fig. 4). The leads are perpendicular to the beam, so they do not heat the beam. The current and the frequency are kept constant (30 A, 200 kHz). The temperature distribution is computed with equation (21) as a function of the time (Fig. 10). The initial shape of the beam (t = 0) is the one on top left of Figure 10. This shape remains the same until the temperature reaches the transformation temperature (50°C here). The simulation shows that this temperature is reached at 0.9 s: it is consistent with the experimental beam deformation, which happens at 1 s (under the same conditions). The shape    No-p5 change is done with the physical support of the elastica from the original shape to the final shape (top-down left column, then top-down right one). The temperature is maximum at the centre, where the inductor is located. Since the diffusion time is high with respect to the length of the beam, the ends of the beam are not heated. After 2 s, the maximum temperature is about 75°C.

Conclusion
A T À V model adapted to thin shells with the help of a change of coordinates has been introduced and confronted with experiment. The elastica, which is used as physical support allows to find the eddy current and temperature distribution as a function of time. It can be used to build a more elaborate mechanical model to take into account the large deformation.