Abstract
This paper considers and reviews a number of known phenomenological models, used to describe hysteretic effects of various natures. Such models consider hysteresis system as a “black box” with experimentally known input and output, related via formal mathematical dependence to parameters obtained from the best fit to experimental data. In particular, we focus on the broadly used Bouc-Wen and similar phenomenological models. The current paper shows the conditions which the Bouc-Wen model must meet. An alternative mathematical model is suggested where the force and kinematic parameters are related by a first-order differential equation. In contrast to the Bouc-Wen model, the right hand side is a polynomial with two variables representing hysteresis trajectories in the process diagram. This approach ensures correct asymptotic approximation of the solution to the enclosing hysteresis cycle curves. The coefficients in the right side are also determined experimentally from the hysteresis cycle data during stable oscillations. The proposed approach allows us to describe hysteretic trajectory with an arbitrary starting point within the enclosed cycle using only one differential equation. The model is applied to the description of forced vibrations of a low-frequency pendulum damper.
1. Introduction
The hysteresis features of a given process (shape of the loop trajectories, asymptotical resemblance, symmetry of forward and reverse processes, etc.) are determined by physical nature [1–11] of the process. Often the dynamics of complex mechanical systems is not straightforward and not easy to model and may depend on a large number of interactions.
Alternative approach is to view the system as a “black box,” with known input and output parameters. Then the interactions between them can be established phenomenologically and their values can be identified experimentally [12–14]. Phenomenological approach is effective in creating a general mathematical model, which can be “migrated” from one field to another. Not only does this show the “depth” of mathematics per se, but more specifically, the overall similarity of various hysteretic processes in nature.
Among widespread phenomenological models that are used to describe hysteresis in various fields of science and technology, the models based on the spectral factorization by relay nonlinearities are particularly important. Originally such approach was proposed in 1935 by German physicist Preisach [15–17], who treated the magnetization process as statistical result of alternating magnetization of individual elementary regions, domains. It is believed that a domain can be in a state of saturation with the direction of magnetization along or against the external field. Accordingly, the magnetization of each domain is described by the switch functions determining the rectangular hysteresis loop. An important component in the model is the particle distribution function, which determines values of magnetization in an arbitrary field.
Later Krasnosel’skii and Pokrovskii and their followers proposed a rigorous mathematical algorithm [18] based on Preisach’s ideas. Similar phenomenological concepts were drawn up and are being developed in different fields of physics and mechanics [19–22]. A common drawback of these models is their parameter identification which requires complex experiments and data analysis.
In 1967, Bouc suggested a solution to the problem of induced vibrations in a mechanical system whose restoring force hysteresis was displacement-dependent [23]. Hysteretic trajectory is described by standard nonlinear first-order differential equation, whose coefficients are identified experimentally. In 1971, he proposed a model for an abstract physical system which he treated as a “black box” with known input and output parameters [24]. In 1976, the model was generalized by Wen [25, 26] and it is known as the Bouc-Wen model since that.
The differential Bouc-Wen “black box” model [12] became very popular due to the fact that it allows analytical description of various hysteretic loops of a damped system [27]. The model has been successfully applied to piezoelectric elements [28], magnetorheological shock absorbers [29], wooden junctions [30], isolation of structure foundations [31], and many more cases.
Let us consider a hysteretic system that transforms a time-dependent input signal into an output signal . In accordance with the Bouc-Wen model, a first-order nonlinear differential equation with switch functions is used to depict . The general form that links input and output signals is [32]where is the selected piecewise-smooth function, which is experimentally identified by a reference signal.
Equation (1) is a part of a general system of dynamic equations, where and are the unknown functions, and it can be used for description of a mechanical system with hysteretic energy dissipation. Let us consider a one-dimensional oscillator as an illustration. In this case, the dynamic equations can be written as [33, 34]The sum is damping force , where is viscous damping coefficient, is stiffness coefficient, is variable describing the hysteretic trajectory, and is constant components of . The values of , and are all identified experimentally. Various modifications of (2) are known, which allow the description of many physical hysteretic processes with relative accuracy [12, 27–37].
Using (1) and (2), it becomes possible to build a dependency that determines the piecewise-smooth continuous hysteretic trajectory. The value of , when , forms the continuity of points , where is the sequence number that grows with time t from the starting point. As it crosses these points, the derivative successively changes sign and the hysteresis branches shift, as it is depicted in Figure 1. Bold arrows in the 2nd and 4th quadrants represent the forward and reverse hysteresis directions, corresponding to the rise and fall of , when, respectively, or .
The identification method is useful for determining model parameters. The error (divergence) between the experimentally acquired output values and those calculated by model algorithms remain relatively small. Calculations are carried out for the reference input signal, after which hysteresis is modelled for other input signals.
There are known cases when the calculated Bouc-Wen model parameters do not agree with the experimental results for other input signals. This demonstrates the ambiguity of identification, which can lead to model instability, in relation to the input signal.
There are established conditions, to which the Bouc-Wen model must adhere [12, 32–34]. The most important are the adequacy of the chosen mathematical model to the physical process and its stability. A model is considered stable on input if the bounding of the input signal corresponds to the bounding of the output signal (BIBO stability).
In the present paper, we establish the system parameters describing a hysteresis process given by (1) with the rhs written as a two-variable polynomial. The polynomial coefficients are determined from experimental data for the envelope cycle curves, which encloses all possible hysteretic trajectories. Figure 1 shows the envelope cycle by two dashed lines and , corresponding to the forward and reverse processes.
We use computational experiments to show that the proposed model features asymptotical stability and allows description of complex hysteretic trajectories that are adequate to a real process.
2. Hysteresis Process Properties
It is well known that dependencies between force and kinematic parameters in a hysteretic system are of cyclic nature. The deformation diagrams show that the path of each cycle is a loop shape formed by two curves (branches), which correspond to the increase and decrease of an external action (load) parameter with time . The initial cycle starting point is determined by the prehistory of the transient process and can be anywhere in the space of variations of the examined parameters. If the oscillations are destabilized, their hysteresis loops can be both opened or closed, differing from each other in shape and relative position.
Experiments show that, in a broad range of hysteresis systems, single type branches of local cycles tend asymptotically to the relative curves of the envelope cycle during monotone process parameter change.
The physical model used to describe such hysteresis by the differential approach (1) adheres to the following conditions.
Condition 1. It is assumed that any physically possible hysteresis trajectory is bound within the space of the envelope cycle, which is experimentally constructed for the maximum range of process parameter changes .
Condition 2. The whole hysteresis process is considered to be independent of frequency. The envelope cycle curves can be thus obtained from quasi-static experiments.
Condition 3. It is also assumed that all local forward process curves are asymptotically related to the curve, and, respectively, the reverse process curves are asymptotically related to . Thus, during a monotone change of parameter , each local curve will tend to the respective envelope cycle curve. For example, Figure 1 shows a curve with starting points at that tends to as decreases. Similarly, a curve with starting points tends to as rises. Generally, all of the local curves, when extrapolated, or , will meet at and in accordance with process direction.
The suggested approach is using a standard first-order differential equation (1) that establishes the dependency for each branch of the hysteresis trajectory. The rhs is selected from a class of functions that provide asymptotical approximations of solution to the respective or of the enclosed cycle. In this case, it is possible to describe an infinite amount of direct and reverse process branches, with differing starting points, but all tending towards or as dictated by process direction.
3. Hysteresis Differential Equation
A standard first-order differential equation is used to model the two hysteresis trajectory branches, with the rhs being a polynomial of two variables and :
coefficients are determined by methods of approach using experimental data for the two branches of the envelope cycle. The matrix has thus two sets of values, and , respectively, with their selection being
Equation (3) must be connected to the system motion equations for time integration under the selected starting conditions. For this, the derivative is transformed into and is added to the lhs of (3). Multiplying both sides by gives the following differential relation:which complements the system’s equations of motion. The sign determines hysteresis direction as well as the choice of constants in conjunction with (4).
The direct and reverse processes are described by (4). At a point , where , the direction is reversed with the continuity of preserved. This change alternates during the time integration.
Because of (4), the resulting two equations of (5) can be united by a switch function , where , which represents speed. Therefore, a single equation, a particular form of (1), replaces (5) and (4), which models the process in both directions.
In most cases, the forward and reverse process curves are identical, differing only in the direction of . The envelope cycle curves are thus symmetrical with regard to the starting point at the coordinate system zero, as illustrated in Figure 2. The diagram shows that for a point at on the forward process there is an opposite point at on the reverse process curve and vice versa.
The reverse process can thus be described by the forward process equation
With variables , . The reverse transformation to and , instead of (7), gives
In this case, the equations of direct and reverse processes can be united and written as
4. Evaluating Parameter of the Model
As stated above, coefficients from (3) are evaluated using the method of the best fit. We minimize the discrepancy representation of into a polynomial with experimentally found and variables that make up envelope cycle curves.
Let us assume that experimental data gave us a sequence of points , which represent one of the envelope curves. To make the sequence , it is possible to use a finite difference expression for the derivative. However, such an approach becomes unfeasible if there are too few data points or their spread is too large.
Instead, we use a reverse method, analytically approximating the hysteretic branch’s derivative, by use of a polynomial
Minimizing the square functional gives us a system of equations to determine where ,. Moreover, if , then .
Thus, from (10), it follows that , where is determined by solving (11). Assuming that , where , we arrive at the desired amount of points .
Let be the value of . We can now acquire a set of these values using the method outlined above.
In accordance with the ordinary least squares method, we construct a discrepancy function as quadratic functionalwhere is number of data points.
Minimization of (12) gives a system of algebraic equations in relation to values of where
Equation (14) can be simplified by introducing a vector of unknowns with elementsand a rhs vector with elementsThus, instead of (14), we have where the quadratic matrix elements are
From (18), the vector whose elements are coefficients in accordance with (16) can be determined. Calculation of is done for the forward and reverse envelope curves with known and . the upper limits of polynomial summation in (3).
and parameters are obtained from experiments to allow the best approximation of the envelope cycle curves and the rate at which local curves approach them. Calculations are undertaken for each combination of and , limiting them with several maximum allowed values and . are thus found experimentally for each combination of and . Equation (3) is then integrated by , which is taken as an independent variable. A starting point for an arbitrary branch is set within or at the border of the envelope cycle. In case of the latter, integration of (3) gives a better approximation for the envelope curves if compared to the quadratic polynomial approximation (10). Examples of such computational experiments are described in [11, 36].
5. Pendulum Torsional Damper and Detuner
One practical example of natural hysteresis is the galloping of the power lines, which is a major concern for both engineers and utilities [14, 35–37]. It is caused by strong winds and icing, which invoke self-oscillations, reminiscent of aircraft wing flutter.
To counteract the potentially dangerous phenomenon, a range of devices are used. Let us examine a pendulum Torsional Damper and Detuner (TDD) which dampens and mismatches vertical and torsional oscillations with energy dissipation hysteresis.
Figure 3 depicts two TDD versions for two- and three-phase lines. In this paper, we will analyse the latter (Figure 3(b)). A general schematic of this TDD is shown in Figure 4(a), and a kinematic diagram is presented in Figure 4(b).
(a)
(b)
(a)
(b)
The main elements, which make up the TDD, are the damping unit 1; curvilinear pendulum rods 4 with weights 2 and counterweights 3; and upper rods 5 to suspend the dampener from the wires, by anchor type fittings 6.
The damping unit (shown disassembled in Figure 5) consists of one or more driving discs 7, suspended from the wires and a number of driven discs 8 that are joined to the pendulum. A central axle 9 allows the discs to rotate about each other. The discs inner surface features special cavities 10 which house elastomeric balls 11 that link the driving and driven discs. Wire movement is passed via rod 5 to the driving discs. Due to the pendulum’s static inertia, the driving disks rotate about the driven ones. This causes displacement of the balls in the cavities between the discs. Because displacement is limited by the cavity shape, motion causes ball deformation and friction which leads to hysteresis.
In spite of the device’s simple design, the energy dissipation mechanism is rather complex as it consists of several simultaneously operating independent processes. These include nonlinear elastomeric deformations, friction, and sliding displacement which generate heat. As said above, it is possible to avoid this laborious task of describing TDD’s hysteresis based on the fundamental physical processes involved by applying the phenomenological method, which we will do in the following chapter.
6. Experiment Results
The experimental research of the TDD prototype was carried out at JSC “Elektrosetstroyproyekt” [14].
Experimentally, we have obtained a set of trajectories where the torsional moment is hysterically dependent on , the angle of disk rotation. The experiment sets were grouped by installation parameters, which characterize the “stiffness” of the damping units. These parameters are the amount of disk “sandwiches,” the amount of rubber balls and their distribution in the cavities, and the torque with which they are compressed between the disks. The damping unit was loaded both manually and automatically.
Experiments showed the following:
Hysteretic trajectories make up the envelope cycle region formed by the forward and reverse process curves.
Hysteresis loops are independent of TDD’s pendulum oscillation frequency. Damping is thus done mainly by hysteresis whilst viscosity contribution to energy dissipation is small.
During monotone load parameter change, every local curve tends to the corresponding envelope cycle curve. All of the local curves in the same direction meet at angle , its limiting value.
Forward and reverse processes are symmetrical.
The rate at which local curves asymptotically approach the envelope depends on installation parameters and physical properties of the damping elements.
A selection of hysteretic dependencies, obtained from quasi-static TDD tests, is shown in Figure 6.
From the experiments, we obtained data for the envelope cycle curves and , which turned out to be identical, that is, symmetrical about the origin. It is thus possible to model hysteresis using (9).
coefficients in (4) were found by minimizing the rhs discrepancy in (3) as a polynomial to the set of values we obtained from the experiments. All of the computations were carried out for different combinations of and , as described in (16)–(19). We eventually selected and , as these allow best approximation of the experimental data. The quadratic minimization of potential (12) gave us the following coefficients: ; = 1422.2; = 10.027; = 1.088; = 0.103; = 0.015; = −23.9; = −4.1; = −2.0; = −0.35; = 143.2; = 419.8. If the angle is measured in radians, then coefficients will be in [N·m], whilst are unitless.
It should be noted that the upper borders of summing and determine the different possibilities of hysteresis curve approximations. The least squares method of their different values give the different values of coefficients. Test calculations were carried out with different combinations of (from 2 to 8) and (from 2 to 4). The results were controlled by the magnitude of the discrepancy function and, consequently, the analytical curve of experimental data. Good approximation is achieved with and combinations. and provide the simplest approximation, which was chosen for comparison with others. The calculations showed that discrepancy functions with differ from other discrepancy functions no less than 10–12%, which is acceptable for engineering calculations.
7. TDD Mathematical Model
A principle TDD schematic is shown in Figure 7 (numeration is continued from Figure 5). Disk 7 is suspended from the wires that require damping and is thus the driving element. Physically it shares a rigid axle with disk 12 that allow them to rotate about one another. Compressed between them are rubber balls 11, which resist this rotation. Attached to disk 8 is the pendulum. Disks 8 and 12 are connected to each other via an elastic element 13, which on this example is a coil spring.
, , and , respectively, represent the turning angles of the driving and driven discs, and the pendulum, which are measured from the vertical in the positive direction counterclockwise.
The generated torque is a result of interaction between the driving disc 7 and the driven disc 12 via the system of damping elements 11. Let us depict it as , where . Torque is the result of the curling element 13 and is considered proportional to angle . By Hooke’s law , where is the torsional stiffness of 13 and is its length. During rotation, there are inertial moments of disks 12 and 8, , respectively.
From the above considerations, the vibration equations of the coupled discs can be written aswhere is combined mass of pendulum and disk 8, is center of mass (CoM) eccentricity, and is acceleration due to gravity.
To model hysteresis and evaluate TDD efficiency in energy dissipation, it is necessary to define the motion of the driving disk. It is accepted that , where and are the amplitude and angular frequency of forced harmonic oscillations, with , where is frequency in Hz.
We can now rewrite (20) along with (9):where . This whole system (21) can be numerically integrated at starting conditions
Damper efficiency can be evaluated in terms of energy dissipation power.
Differentiating (23) with respect to time gives , from which we obtain an equation that is easy to integrate over time along with (21):
The solution of (24) has a horizontal asymptote that corresponds to the dissipation power of stabilized forced vibrations.
8. Modelling Results
Below are the TDD energy dissipation analysis results for two design concepts shown in Figure 8. They are labelled in accordance with their functions I, inertial, and G, gravitational dampers.
For the I-damper, the differences between axis of rotation and weighs’ center of mass are m, m, with weights masses being kg, kg. For the G-damper the lengths are m with mass of kg. Driving disk moment of inertia is taken as kg·.
Figures 9(a)–14(a) show the integration results of (21) with the starting conditions (22). Integration is done over the s range, which was selected to match the time interval of stable vibrations.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
Figures 9(b)–12(b) show the hysteretic trajectory modelling under harmonic oscillations of driving disk with amplitude at rad and at frequencies = 0.2, 0.4, 0.6 and 0.8 Hz. The smooth approximations of the envelope cycle and curves are shown by dashed lines. Figures 9(b)–12(b) show plots of pendulum deviation angle from the vertical.
Figure 13 depicts the integration of (23). It is clear that dependency curves approach asymptotes with time, the ordinate of which is the energy dissipation power of steady oscillations in I- and G-dampers.
Power calculations of at different frequencies led to dependencies. These are shown in Figure 14(a). Up to 0.5 Hz the G-damper appears to be more effective than the I-damper. However, the power of both dampers never exceeds 1 W, which suggests low effectiveness of both versions. Between 0.65 Hz and 0.85 Hz, the energy dissipation of the I-damper is sufficiently larger than the G-damper. Up to 0.75 Hz G-damper is practically nonfunctioning, whereas I-damper dissipation is already 40 W. Above 1 Hz, the dissipation power of both dampers is identical. It can be thus concluded that the I-damper is effective in the Hz frequency range, where wire galloping begins.
The results are illustrated in Figure 14(b), where the transfer function depends on the frequency of forced harmonic oscillations in the driving disk. By transfer function we mean the relation of pendulum oscillation amplitude in regard to the driving disk (or the support frame) to the amplitude of driving disk’s harmonic oscillations. At 0.65 Hz the amplitudes of I-damper’s driving and driven disks equalize, whereas the G-dampers disks are almost synchronized.
Figure 15 depicts a diagram of a twin pendulum damper, whose rods are at right angle to each other. The driving disc’s moment of inertia is the same as before kg·. For each weight, m, kg.
This damper was also built and tested at JSC “Elektrosetstroyproyekt.” Its theoretical analysis was based on the method outlined above. Also the free oscillations of the driven disc and pendulums, with the driving disc immobilized, were studied.
Figure 16 shows the results with initial deviation of the driving disc at 0.4 rad without a staring speed. The dependence of the disc rotation angle on time is shown in Figure 16(a), where the straight and dashed lines represent the smooth approximation from experimental data and numerical integration of (21). Figure 16(b) shows the hysteresis .
(a)
(b)
9. Conclusions
The paper develops a kinematic approach to the description of hysteresis in a range of nonstationary processes. The approach is related to the famous Bouc-Wen method, where a mathematical model is based on a standard first-order differential equation. However, unlike the Bouc-Wen, the rhs is depicted as a two-variable polynomial, the time-dependent hysteresis parameter and a function of this parameter. These variables can be translation (angular) and corresponding force (torque) to real-life scenario. The polynomial coefficients are found by analytical approximation of the envelope cycle curves during stabilized oscillations.
The suggested model of physical dependencies is analytical, which makes it suitable for description of nonlinear mechanical systems.
The approach can be used to solve problems related to nonstationary oscillations in mechanisms with hysteretic type of energy dissipation. This is demonstrated on the example of a pendulum type low-vibration damper. An algorithm to analyze its effectiveness in counteracting wire galloping was also outlined.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was carried out with financial support from the Russian Ministry of Education and Science (Federal Programme no. 14.604.21.0188).