Study of Self-Excited Thermomechanical Oscillator with Shape Memory Alloys

1. Introduction

Self-excited oscillations accompany many technological processes and, in most cases, have detrimental effects on them. The main reasons for their occurrence are mainly the dependence of friction forces on relative sliding, changes of normal loads, viscoelastic interactions or the influence of various geometrical effects [1,2,3].

It has been shown that purely temperature self-excited oscillations cannot exist because the heat conduction equation, unlike other equations of mechanics, is not invariant under time reversal. Periodic thermally forced self-excited oscillations that propagate as waves are possible [4]. The combination of thermal and mechanical vibrations is known by the term thermo-mechanical oscillations (TMO). The foundations of the theory of TMOs were laid in 1829 by A. Travelian [5], and was later further developed by M. Faraday, A. Tyndal, and J. Rayleigh. Some of the first theoretical and experimental studies of the conditions for the occurrence of self-excited thermomechanical oscillations were reported by Nesis E. I. [4]. Awrejcewicz and Pyryev [6] created a massless elastic element model of frictional self-excited vibrations and investigated the possibilities of their existence.

Thermomechanical self-excited oscillations have a negative influence on many technological processes. To show that the loss of stability in constant speed cutting is due to thermomechanical self-excited oscillations authors Astashev and Korendyasev [7] investigated the model of heat generation in the cutting region. Fel’dshtein [8] investigates self-excited thermomechanical oscillations of current-carrying wires by determining the existence conditions and describing them by numerical simulations. Toda et al. [9], using the Barenblatt model, investigate thermomechanical self-excited oscillations that occur during cold drawing of films. A technique to stabilize and suppress the oscillations in an opto-thermomechanical oscillator is demonstrated in [10].

In microelectromechanical systems, (MEMS) thermal processes occur at a much higher rate than in conventional macrosystems due to the significantly lower values of temperature time constants. And in these systems, one can note detrimental effects of self-excited thermomechanical oscillations, such as those induced by laser heating of double-clamped MEMS beams [11], or the induced frequency self-locking that occurs when laser illuminating thin MEMS resonators [12,13].

The micrometric dimensions of thermally-exposed components in MEMS create conditions for building fundamentally new devices in which self-excited thermomechanical oscillations are the basis of their operating principle. One such large class of MEMS are self-excited thermomechanical microresonators [14,15] which are in competition with their capacitive [16] and piezoelectric analogues [17].

Shape memory alloys (SMAs) are smart materials possessing the property of remembering a shape, which if changed within certain limits, can be fully recovered by heating. When the change of the shape is a stretch due to tensile stress, the SMA specimen tends to shorten to its original undeformed dimensions when it is restored. This behavior is used to create relatively low-cost actuators that possess several advantages such as generating large forces and displacements at high energy density and noiselessness. Thermal activation of SMAs creates conditions for cyclic control [18,19,20] and in many cases for self-excitation of thermomechanical oscillations [21,22,23,24].

Oscillating systems based on shape memory alloys are widely used in energy harvesting technologies. Several researchers use SMA substrates on which they deposit a piezoelectric layer to generate electricity by vibration [25,26,27,28]. Another interesting application of SMAs in energy harvesting systems is the two-step conversion of heat to mechanical and then mechanical to electrical energy [29,30,31,32]. A peculiarity of thermomechanical models with SMAs are the nonlinearities arising from the hysteresis in the stress-strain, stress-temperature, and temperature-strain relationships. The existence conditions of oscillatory modes of nonlinear thermomechanical self-excited systems with hysteresis are a challenge for modern research in the field of differential equations and control theory [33,34,35].

The aim of this paper is to create a mathematical model of a new self-excited thermomechanical oscillator with SMA filament, to validate the obtained results experimentally and to draw conclusions about the type and parameters of oscillations. In the development of the mathematical model, the aim is to describe the thermodynamic processes in detail, considering the hysteresis phenomena in the SMA filament.

2. Design Concept of the Self-Excited Thermomechanical Oscillator

The schematic of the self-excited thermomechanical oscillator is shown in Figure 1. The structure is symmetrical consisting of rigidly fixed at both ends SMA filament OABC on which two metal spheres in points А and B are arranged. Nitinol (NiTi) SMA was used for the filament [36]. A heater is fixedly located at the base of the thermomechanical oscillator.

The principle of the oscillator is as follows. Under the action of the weight of the two spheres, the SMA filament sags and approaches the heater at position OA₁B₁C. Тhe proximity of the filament to the heater raises its temperature and causes contraction while at the same time harden it. The contracted filament rises back vertically under the temperature, moving away from the heater at position OA₀B₀C. This is followed by a lowering of the temperature, softening and relaxation of the filament. As a result, under the weight of the spheres, the filament falls again and approaches the heater at position OA₁B₁C. This is followed by a cyclic repetition of the above processes, resulting in periodic oscillations of the system.

3. Building a Dynamic Model and Determining the Input Data for Its Solution

The SMA filament in the static position is pre-stretched between the OC supports. The distance between the supports is assumed to be OC = 2l and the length of the filament before tensioning is 2l_s. Thus, for half of the prestrain of the SMA filament Δls₀ is valid:

$$\Delta l_{s0}=l-l_s,$$

(1)

The x-axis of the coordinate system Oxy is chosen to lie on the segment OC and its origin is located in the middle of the same segment. Symmetrically to the y-axis on the OC segment, two points A₀ and B₀ are selected at which two spheres of equal mass m are fixed. The distances along the horizontal x-axis are denoted by с $l_{10}=O{A}_0$, $l_{20}=A_0{B}_0$ and $l_{30}=B_{0}C$. Since for the relative lengths, the notations

$$\frac{l_{01}}{l}=\frac{l_{03}}{l}={\lambda }_0,$$

(2)

are chosen then for the deformations from the prestressing $\Delta l_{s10}$ and $\Delta l_{s30}$ in the end sections OA0 and B0C will follow:

$$\Delta l_{s01}=\Delta l_{s03}=\Delta l_{s0}{\lambda }_0,$$

(3)

and in the middle section A₀B₀ this deformation will be

$$\Delta l_{s012}=2\Delta l_{s0}\left(1-{\lambda }_0\right),$$

(4)

It can be seen in Figure 1 that due to the symmetrical arrangement of the heater and the CMA filament, only the middle section of the AB will be heated. It is assumed that the temperature of this section is distributed uniformly along its length.

Since the dynamic model will be constructed using the Lagrange equations of the second kind, the kinetic energy of the oscillator is represented in the form

$$E_k=\frac{m}{2}\left({\dot{x}}_A^2+{\dot{x}}_B^2+{\dot{y}}_A^2+{\dot{y}}_B^2\right),$$

(5)

where ${\dot{x}}_A,\hspace{0.17em}\hspace{0.17em}{\dot{x}}_B$ are the velocities of points A and B on the x-axis and ${\dot{y}}_A,\hspace{0.17em}\hspace{0.17em}{\dot{y}}_B$ are the velocities of points A and B on the y-axis.

The end sections OA and BC are assumed not to heat up, unlike the middle AB section, which changes its temperature and Young’s modulus, so the potential energy takes the expression

$$E_p=\frac{k_l}{2}\left(\Delta l_1^2+\Delta l_3^2\right)+\frac{k_T}{2}\Delta l_2^2,$$

(6)

where $k_l$ is the elastic constant of the extreme cold sections, $k_T$is the elastic constant of the middle warm section,

$$\Delta l_1=\sqrt{{\left(x_A-x_O\right)}^2+{\left(y_A-y_O\right)}^2}-l_{10}+\Delta l_{s0}{\lambda }_0,$$

(7)

$$\Delta l_3=\sqrt{{\left(x_B-x_C\right)}^2+{\left(y_B-y_C\right)}^2}-l_{30}+\Delta l_{s0}{\lambda }_0,$$

(8)

are the deformations of the end sections, $x_O$, $x_A$, $x_B$, $x_c$ $y_O$, $y_A$, $y_B$, $y_c$ are the coordinates of points O, A, B and C.

For the middle section, the deformation depends on the mutual position of points A, B and the pre-deflection of the filament:

$$\Delta l_2=\sqrt{{\left(x_A-x_B\right)}^2+{\left(y_A-y_B\right)}^2}-l_{20}+2\Delta l_{s0}\left(1-{\lambda }_0\right),$$

(9)

The stiffness of the end sections is

$$k_l=\frac{E_{1}A}{l_{01}},$$

(10)

where A is the cross-sectional area, and the Young’s modulus $E_1$ is

$$E_1=E_m,$$

(11)

where $E_m$is the NiTi module in the martensitic state. This is the case if it is assumed that the temperature in the final sections is lower than the final martensitic temperature, i.e., T<Mf, the Young’s modulus according to [37] is equal to the modulus of NiTI in the fully martensite state.

For the middle section, the stiffness is

$$k_{Т}=\frac{E_{02}A}{l_{20}},$$

(12)

where the Young’s modulus $E_{02}$ is defined by the piecewise function

$$E_2=\{\hspace{0.17em}\begin{array}{c}\begin{array}{l}{E}_a-\left(E_a-E_m\right)R_m\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le {\epsilon }_{s2}\le {\epsilon }_m^y\\ E_a-\left(E_a-E_T\right)R_m\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_m^y\le {\epsilon }_{s2}\le {\epsilon }_m^d\\ E_a-\left(E_a-E_d\right)R_m\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_m^d\le {\epsilon }_{s2}\end{array}\end{array},$$

(13)

where $E_m$, $E_T$, $E_d$, and $E_a$ are Young’s modulus of the fully twined, partially twined, detwinned martensite and austenite respectively; ${\epsilon }_m^y$ is the yield strain of the twined martensite; ${\epsilon }_m^d$is the minimum strain of the detwinned martensite and

$${\epsilon }_{s2}=\frac{\Delta l_2}{l_{02}}=\frac{\sqrt{{\left(x_A-x_B\right)}^2+{\left(y_A-y_B\right)}^2}-l_{20}+2\Delta l_{s0}\left(1-{\lambda }_0\right)}{l_{20}},$$

(14)

is the current strain of the middle section of the filament.

The relative martensitic fraction $R_m$ of NiTi is calculated according to [32,37], and to investigate the major and minor hysteresis similarly to [18] it is modified in the form

$$R_m=\frac{S_R}{2}\left(1+\tanh \frac{T_m-T}{{\lambda }_T}\right)+R_C{S}_R,$$

(15)

where the temperature $T_m$ depends on whether the thermal process is heating or cooling and is obtained by the relation

$$T_m=\frac{1}{2}\{\begin{array}{c}{A}_S+A_f\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{dT}{dt}>0\\ M_s+M_f\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{dT}{dt}<0\end{array},$$

(16)

the coefficient ${\lambda }_T$has unit °C and here it has the form

$${\lambda }_T=\frac{1}{2}\{\begin{array}{c}\frac{A_S+A_f-b_a}{c_a}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{dT}{dt}>0\\ \frac{M_s+M_f-b_m}{c_m}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{dT}{dt}<0\end{array},$$

(17)

and As Af Ms Mf are the start and final temperatures for austenite and martensite of SMA, the limiting values of the coefficients, $b_a<0.5({А}_f-A_s)$, $b_m<0.5(M_s-M_f)$, $c_a<0.5$, and $c_m<0.5$ have been established experimentally [37].

The major and minor hystereses of the relative martensitic fraction are determined using the method described in [18], adapted for the model of [37]. In summary, this method is based on the points at which the temperature process changes from warming to cooling and vice versa. According to the adopted method, if the temperature gradient change point $T_f$ is in the transition periods when both martensite and austenite fractions are present in the СМА i.e.,

$${Т}_f\in \left[M_s,M_f\right]\cap \left[A_f,A_s\right],$$

(18)

the scale factor of the martensitic fraction acquires the values of the case of transition from heating to cooling equivalent to

$$\hspace{0.17em}{\dot{T}}_{f-}>0\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}{\dot{T}}_{f-}<0\hspace{0.17em}\hspace{0.17em},$$

(19)

where ${\dot{T}}_{f-}$ and ${\dot{T}}_{f-}$ are the temporal temperature gradients before and after the onset of the fluctuation, respectively the scale factor of the martensitic fraction has the form

$$S_R=\{\begin{array}{l}1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_S>T_f>M_f\\ R_m^f=\hspace{0.17em}\frac{1}{2}\left(1+\tanh \frac{T_m-T_f}{{\lambda }_T}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{A}_f>T_f>A_S\end{array},$$

(20)

аnd the fraction

$$R_C=\{\begin{array}{l}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_S>T_f>M_f\\ 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{A}_f>T_f>A_S\end{array},$$

(21)

If at the change point the temperature changes from cooling to warming

$$\hspace{0.17em}{\dot{T}}_{f-}<0\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}{\dot{T}}_{f-}>0\hspace{0.17em}\hspace{0.17em},$$

(22)

The scale coefficients are respectively:

$$S_R=\{\begin{array}{l}\hspace{0.17em}1\hspace{0.17em}-\hspace{0.17em}{R}_m^f\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_S>T_f>M_f\\ R_m^f=\hspace{0.17em}\frac{1}{2}\left(1+\tanh \frac{T_m-T_f}{{\lambda }_T}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_s\ge T_f\ge M_f\end{array},$$

(23)

and

$$R_C=\{\begin{array}{l}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_f>T_f>A_f\\ 1-R_m^f\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_s>T_f>M_f\end{array},$$

(24)

A sub minor hysteresis occurs if a minor hysteresis already exists after the point $T_f$. The temperature at which sub minor hysteresis occurs is denoted by $T_f^s$. The other important condition for sub minor hysteresis existing is when cooling ${М}_f<T_f^s<M_S$, and when warming ${А}_S<T_f^s<A_f$

For a sub minor hysteresis, the following scale factors are valid for the warming-to-cooling transition:

$$S_R=\{\begin{array}{l}1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_f>T_f^S>A_f\\ R_m^{fs}=\hspace{0.17em}\frac{1}{2}\left(1+\tanh \frac{T_m-T_f^S}{{\lambda }_T}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{A}_f>T_f^S>A_S\end{array},$$

(25)

and

$$R_C=\{\begin{array}{l}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_S>T_f^S>M_f\\ 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{A}_f>T_f^S>A_S\end{array},$$

(26)

At the transition from cooling to warming the scale coefficients are

$$S_R=\{\begin{array}{l}{R}_m^f\hspace{0.17em}-R_m^{fs}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_f>T_f^S>A_f\\ R_m^{fs}=\hspace{0.17em}\frac{1}{2}\left(1+\tanh \frac{T_m-T_f^S}{{\lambda }_T}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_s>T_f^S>M_f\end{array},$$

(27)

and

$$R_C=\{\begin{array}{l}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_f>T_f^S>M_s\\ 1-R_m^{fc}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{M}_s>T_f^S>M_f\end{array},$$

(28)

To determine the heater temperature distribution in the area where the SMA filament moves, the experimental setup shown in Figure 2 was designed. The temperature in the vicinity of the heater 3 was measured using a movable thermocouple 1, the position of which was set using an inductive position sensor 5 type APO-075-002-000 of TT Electronics Ltd. UK. The constant temperature of the heater is measured by a thermocouple 2 which is fixed in its body. A mechanical dial gauge 4 was also added to the experimental system to calibrate the data of the inductive sensor 5. The data from the two thermocouples is fed to a Data Acquisition System (DAQ) type cDAQ -9191, and the inductive sensor signal is fed into a DAQ type USB -9191 of National Instruments Ltd. The data is processed and recorded by a LabVIEW program. Different values of the constant temperature of the heater are set using a thermostat and by changing the position of the movable thermocouple, the type of heat field in the vicinity of the heater is determined.

Based on the experiments, the temperature of the filament in the middle part is assumed to be a stationary function of the y coordinate defined by the piecewise function

$$T=\{\begin{array}{l}\left(T_h-T_0\right)e^{-{\left(\frac{y-y_e}{\theta }\right)}^2}+T_0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}y<y_e\\ T_h\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}{y}_e\le y\le y_f\\ \left(T_h-T_0\right)e^{-{\left(\frac{y-y_f}{\theta }\right)}^2}+T_0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}y>y_f\end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},$$

(29)

where $y_e,\hspace{0.17em}{y}_f$ are coordinates dependent on the heater position determined experimentally (Figure 3), $T_h$is the maximum temperature received by the filament from the heater, $\theta$ is a conditional length defining the slope of the curve T=T(y) and $T_0$ is the room temperature.

In Figure 3 the experimental values of the hater’s thermal field for three maximal temperatures 60, 67. are 78 °C are presented with solid lines. The corresponding approximated thermal fields are presented with dashed lines.

The approximated heat field function along the vertical coordinate gives rise in this model to consider the relative martensitic fraction, as a complex function of the $y_{А}$ coordinate and time t, i.e.,

$$R_m=R_m\left(y_A\left(t\right)\right),$$

(30)

From which it follows that the Young’s modulus and stiffness of the SMA filament in the middle section are considered respectively as

$${Е}_2={Е}_2\left(x_A\left(t\right),y_A\left(t\right)\right),$$

(31)

and

$$k_T=k_T\left(x_A\left(t\right),y_A\left(t\right)\right),$$

(32)

The assumption about the temperature dependence of the middle section on the y position leads to the representation of the thermal stiffness as a function of time. The potential character for the temperature distribution near the heater provides a rationale for deriving the filament forces from the SMA to be derived via the second-order Lagrange equations.

It is assumed that the motion of the filament is such that the oscillator is only in its symmetric configuration, where $x_B=-x_{А}$, $y_B=y_{А}$, $x_O=-l$, $x_C=l$,$y_C=y_O$=0, $l_{30}=l_{10}$, $\Delta l_3=\Delta l_1$, and $l=OC/2$.

In this case the Lagrangian of the oscillator takes the form

$$L={Е}_K-E_P=m\left({\dot{x}}_A^2+{\dot{y}}_A^2\right)-k_l\Delta l_1^2-2k_T\left(x_A,y_A\right)\Delta l_{22}^2\left(y_A\right),$$

(33)

where the symmetry-transformed deformation of the middle section is rewritten as

$$\Delta l_2=2\Delta l_{22}\left(x_A\right),$$

(34)

Here with $\Delta l_{22}\left(x_A\right)$ expresses

$$\Delta l_{22}\left(x_A\right)=x_A-\frac{l_{20}}{2}+\Delta l_{s0}\left(1-{\lambda }_0\right),$$

(35)

Which can be rewritten as

$$\Delta l_{22}\left(x_A\right)=x_A+l_{22},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{l}_{22}=-\frac{l_{20}}{2}+\Delta l_{s0}\left(1-{\lambda }_0\right),$$

(36)

The strain of the middle segment is expressed

$${\epsilon }_{s2}=\frac{2(x_A+l_{22})}{l_{20}},$$

(37)

and the deformation of the end section is converted into:

$$\Delta l_1=\sqrt{{\left(x_A+l\right)}^2+y_A^2}-l_{10}+\Delta l_{s0}{\lambda }_0,$$

(38)

In the above formula the expression $-l_{10}+\Delta l_{s0}{\lambda }_0$ is denoted as $l_{11}$and the deformation of second section obtains the form

$$\Delta l_1=\sqrt{{\left(x_A+l\right)}^2+y_A^2}+l_{11},$$

(39)

The differentiation of the Lagrangian is performed with respect to the generalized coordinates $x_A$, $y_A$ and the generalized velocities ${\dot{x}}_A$and ${\dot{y}}_A$, according to the Lagrange equations of the second kind

$$\frac{d}{dt}\frac{\partial L}{\partial {\dot{q}}_i}-\frac{\partial L}{\partial q_i}=Q_i\hspace{0.17em}\hspace{0.17em}i=1,2,$$

(40)

where $\{q_1,q_2\}=\{x_A,y_A\}$, $\{{\dot{q}}_1,{\dot{q}}_2\}=\{{\dot{x}}_A,{\dot{y}}_A\}$и $\{Q_1,Q_1\}=\{-{\beta }_x\dot{x},-2mg-{\beta }_y\dot{y}\}$ are the generalized coordinate, generalized velocity, and generalized force, respectively, and are the damping coefficients along the x and y axes.

After differentiation according to (40) the system of differential equations is derived

$$|\begin{array}{c}m\ddot{x}+k_l\Delta l_1\left(x\right)\frac{d\Delta l_1\left(x,y\right)}{dx}+\frac{\partial k_T\left(x,y\right)}{\partial x}\Delta l_{22}{\left(x\right)}^2+2k_T\left(x,y\right)\Delta l_{22}\left(x\right)\frac{d\Delta l_{22}\left(x\right)}{dx}=-{\beta }_x\dot{x}\\ m\ddot{y}+\frac{\partial k_T\left(x,y\right)}{\partial y}\Delta l_{22}{\left(x\right)}^2-k_l\Delta l_1\left(x\right)\frac{d\Delta l_1\left(x,y\right)}{dy}+\frac{\partial k_T\left(x,y\right)}{\partial y}\Delta l_{22}{\left(x\right)}^2=-mg-{\beta }_x\dot{y}\end{array},$$

(41)

where the derivative of the filament end segment with respect to x is

$$\frac{d\Delta l_1\left(x\right)}{dx}=\frac{x+l}{\sqrt{{\left(x+l\right)}^2+y^2}},$$

(42)

the derivative of the middle filament segment with respect to x is

$$\frac{d\Delta l_{22}\left(x\right)}{dx}=1,$$

(43)

The stiffness $k_T$ do not depend on x and its derivative with respect to x is equal to zero, or

$$\frac{\partial k_T\left(x,y\right)}{\partial x}=0,$$

(44)

The derivative of $k_T$with respect to y is equal to

$$\frac{\partial k_T\left(x,y\right)}{\partial y}=\frac{\partial }{\partial y}\frac{E_{02}A}{l_{02}}=\frac{A}{l_{02}}\frac{\partial E_{02}}{\partial y},$$

(45)

In (41), the partial derivative of Young’s modulus with respect to y for the middle section has the form:

$$\frac{\partial E_2}{\partial y}=\{\hspace{0.17em}\begin{array}{c}\begin{array}{l}\left(E_m-E_a\right)\frac{d{R}_m}{dy}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le {\epsilon }_{s2}\le {\epsilon }_m^y\\ \left(E_T-E_a\right)\frac{d{R}_m}{dy}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_m^y\le {\epsilon }_{s2}\le {\epsilon }_m^d\\ \left(E_d-E_a\right)\hspace{0.17em}\frac{d{R}_m}{dy}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_m^d\le {\epsilon }_{s2}\end{array}\end{array},$$

(46)

The relative martensite fraction taking into account the heater heat field distribution function (29) is modified into

$$R_m=\{\begin{array}{c}\frac{S_R}{2}\left(1-\tanh \frac{T_m-T_h{e}^{-\frac{1}{{\theta }^2}\left(y-y_e\right)}+T_0}{{\lambda }_T}\right)+R_C{S}_R\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}y\le y_e\\ \frac{S_R}{2}\left(1-\tanh \frac{T_m-T_h+T_0}{{\lambda }_T}\right)+R_C{S}_R\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}{y}_e<y<y_f\\ \frac{S_R}{2}\left(1-\tanh \frac{T_m-T_h{e}^{-\frac{1}{{\theta }^2}\left(y-y_f\right)}+T_0}{{\lambda }_T}\right)+R_C{S}_R\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}y\ge y_f\end{array}\hspace{0.17em}\hspace{0.17em},$$

(47)

The derivative of the relative martensite fraction in this case is

$$\frac{d{R}_m}{dy}=\{\begin{array}{l}\frac{S_R{T}_h}{2{\lambda }_T{\theta }^2}{e}^{-\frac{1}{{\theta }^2}\left(y-y_e\right)}\left(1-\tanh \frac{T_m-T_h{e}^{-\frac{1}{{\theta }^2}\left(y-y_e\right)}+T_0}{{\lambda }_T}\right)\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}y\le y_e\\ 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}\hspace{0.17em}{y}_e<y<y_f\\ \frac{S_R{T}_h}{2{\lambda }_T{\theta }^2}{e}^{-\frac{1}{{\theta }^2}\left(y-y_f\right)}\left(1-\tanh \frac{T_m-T_h{e}^{-\frac{1}{{\theta }^2}\left(y-y_f\right)}+T_0}{{\lambda }_T}\right)\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}y>y_f\end{array},$$

(48)

After the simplifications described above, the dynamic model (41) is recast into the system of differential equations

$$|\begin{array}{c}m\ddot{x}+\frac{E_{m}A\left(x+l\right)\left(\sqrt{{\left(x+l\right)}^2+y}-l_{11}\right)}{l_{10}\sqrt{{\left(x+l\right)}^2+y}}+2k_T\left(x,y\right)\left(x+l_{22}\right)=-{\beta }_x\dot{x}\\ m\ddot{y}+\frac{\partial k_T\left(x,y\right)}{\partial y}{\left(x+l_{22}\right)}^2+\frac{E_{m}A\left(\sqrt{{\left(x+l\right)}^2+y}-l_{11}\right)y}{l_{10}\sqrt{{\left(x+l\right)}^2+y}}=-mg-{\beta }_y\dot{y}\end{array},$$

(49)

In order to find a numerical solution of the system (49) its order is lowered by the following substitutions

$$|\begin{array}{c}{\dot{v}}_x=-\frac{E_{m}A\left(x+l\right)\left(\sqrt{{\left(x+l\right)}^2+y}-l_{11}\right)}{m{l}_{01}\sqrt{{\left(x+l\right)}^2+y}}-\frac{2}{m}{k}_T\left(x,y\right)\left(x+l_{22}\right)-\frac{{\beta }_x}{2m}\dot{x}\\ \dot{x}=v_x\\ {\dot{v}}_y=\frac{E_{m}A\left(\sqrt{{\left(x-l\right)}^2+y}-l_{11}\right)y}{m{l}_{10}\sqrt{{\left(x-l\right)}^2+y}}-\frac{\partial k_T\left(x,y\right)}{m\partial y}{\left(x+l_{22}\right)}^2-g-\frac{{\beta }_y}{2m}\dot{y}\\ \dot{y}=v_y\end{array},$$

(50)

The numerical solution of the system was carried out using the parameters given in Table 1.

To solve the system (50), a numerical solution program was created using an implicit 5th order Euler method. The solutions of the system of differential equations (50) were obtained under zero initial conditions for the velocities $\dot{x}$ and $\dot{y}$, and the initial conditions for the coordinates of the masses were taken $x\left(0\right)=-0.1$ and $y\left(0\right)=\mathrm{0.01.}$

4. Conducting the Experiments and Comparing the Results

The displacements of the spheres along the x and y coordinates were measured synchronously in real time using the experimental system shown in Figure 4. A photograph of the system is shown in Figure 4a and its simplified block diagram is depicted in Figure 4b. Filament 1 is immovably clamped on two supports 6. Two spheres 4 are arranged symmetrically with respect to the heater 3 and the fixed supports. On the filament adjacent to the spheres are mounted movable sensing elements of non-contact inductive sensors 2 and 5 which are type APO-075-002-000 of TT Electronics Ltd. UK. In order to detect the simultaneous displacement of the spheres in the horizontal x-direction and the vertical y-direction, the movable sensing elements and the sensor bodies are oriented at 90 degrees to each other. In Figure 4, sensor 2 measures the x-axis sphere displacement and sensor 5 is oriented to measure the vertical y-displacement of the sphere. The signals from the two sensors are input to the Data Acquisition System DAQ USB - 6003 from National Instruments, then processed by a laptop using a LabVIEW program. The data from both sensors and the heater temperature are recorded in an Excel file.

Multiple experiments were conducted by varying the weights masses, the initial tension of the SMA filament and the heater temperature which was kept constant using a thermoregulator.

In the experiments it was found that for oscillations to occur, certain conditions concerning the pre-tension of the SMA filament, the difference in room and heater temperature needed to be met. The most significant influence for the occurrence of oscillations is the initial position of the filament relative to the heater. Once the conditions for the oscillations to occur are satisfied, stable oscillations of the filament with both masses are generally observed with unchanging energy.

Analysis of the experimental results shows that the self-excited oscillations have a stochastic character. The reason for the appearance of the random processes is due to the nonlinear hysteresis variations of the Young’s modulus of the SMA, which are a function of the temperature, the sign of its gradient, the type of deformation and the position of the SMA filament.

Figure 5a shows the plot of the experimental measurements of the x-axis position of the sphere. It can be seen that the oscillatory motion has a period which varies randomly within certain limits while the amplitude also varies according to a random distribution law.

An interesting impression of the x-axis oscillations is an extra modulated oscillation with approximately 20 times the frequency and approximately 5 times the amplitude. This modulation is also confirmed in Figure 5b, which represents the result of the theoretical model for the x-axis oscillations. The logical explanation for this phenomenon, according to the authors, is the existing a purely mechanical normal vibration of the SMA filament that occurs under the action of the variable axial force. This variable axial force is generated by the thermal contractions of the SMA filament. The difference in the frequencies of the two oscillations is due to the fact that the mechanical oscillations have a frequency which depends on the stiffness of the filament and mass, while the thermal oscillations depend on the vertical position of the filament, and primarily on the temperature time constant of the filament determining the heating cooling time. For the same reasons, the amplitudes of the longitudinal mechanical oscillations have a small value, while the thermal amplitudes can reach up to 8% of the filament length. Thus, the x-axis displacement can be assumed to be the sum of mechanical and thermal random oscillatory motions.

Figure 6 plots the displacements of the sphere along the y-axis. The experimental data is shown in Figure 6a and the data obtained from the dynamic model is plotted in Figure 6b. It can be observed that stochasticity is also present in the vertical oscillations both in terms of period and frequency. Modulated oscillations along the y-axis direction are not seen in either the experimental data or the data obtained from the theoretical model. The reason for this is that here the transverse oscillations are purely mechanical, induced by thermal contractions along the x-axis. The appearance of the stochasticity here is also explained by the change of Young’s modulus at different temperatures and strains.

5. Conclusions

A new system for self-excited thermomechanical oscillations based on shape memory alloys was created. The properties of this system are investigated theoretically by a dynamic model in which the hysteresis properties of the SPF are accounted for. The results obtained from the model are confirmed by direct measurements of displacements along two perpendicular directions. Using this experimental method, displacements were obtained with relatively high accuracy, since integration errors typical of inertial measurement systems such as MEMS accelerometers or gyroscopes were eliminated. The obtained results show that the oscillations of the system have a stochastic character in terms of frequency and vibration. It is found that due to the birefringent nature of the oscillations along the x-axis, there are two types of oscillations. The main ones called thermal have an amplitude dependent on the thermal contractions of the filament and a frequency determined by the temperature time constant of the filament from the SPF. The second oscillations called purely mechanical have a frequency dependent on the ratio of the mass of the sphere and the stiffness of the system and an amplitude that is a function of the longitudinal thermal contractions, the natural frequency and the damping of the SPF filament.