Mathematical Problems in Engineering

Volume 2014 (2014), Article ID 497872, 13 pages

http://dx.doi.org/10.1155/2014/497872

Research Article

The Identification of Nonlinear Systems Using - Plane Coordinates

Ming-Fei Chen,¹ Dyi-Cheng Chen,² Chung-Heng Yang,¹ and Shih-Feng Tseng³

¹Department of Mechatronics Engineering, National Changhua University of Education, No. 2, Shi-Da Road, Changhua 500, Taiwan

²Department of Industrial Education and Technology, National Changhua University of Education, No. 2, Shi-Da Road, Changhua 500, Taiwan

³Instrument Technology Research Center, National Applied Research Laboratories, Hsinchu 30076, Taiwan

Received 3 October 2013; Revised 14 December 2013; Accepted 18 December 2013; Published 26 February 2014

Academic Editor: Hao-Chun Lu

Copyright © 2014 Ming-Fei Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In order to instantaneously distinguish the (coefficient of viscous damping) and (coefficient of stiffness), which are both functions of time in an M.C.K. nonlinear system, a new identification method is proposed in this paper. The graphs of the - are analyzed and the dynamic behavior of M.C.K. systems in a - coordinate plane is discussed. This method calculates two adjacent sampling data, the displacement, velocity, and acceleration (which are obtained from the responses of a pulse response experiment) and then distinguishes and of an instantaneous system. Finally, this method is used to identify the aerostatic bearing dynamic parameters, C and K.

1. Introduction

In order to control an M.C.K. system, the stiffness and damping must be known. However, in an actual M.C.K. system, the stiffness and damping are essentially nonlinear. Even though several studies consider the frequency domain [1–6], these analyses do not allow real-time identification. References [7–10] identify the time domain. Huang et al. [11] modified the logarithmic-decrement method to avoid disturbance from measured noise, but the second peak is still required to analyze the damping. As a result, a system that has at least two peaks is an underdamped system. In order to control the system, the displacement response must be almost critically damped and the peaks of its displacement response must be less than one. Therefore, the use of the logarithmic-decrement method is limited. For real-time control systems, the current method requires improvement.

A single degree-of-freedom linear oscillator with viscous damping has the equation of motion:

However, in reality, many such oscillators show some nonlinearity in stiffness and/or damping. This paper presents a simple method to identify this nonlinearity. If the mass is assumed to be known and constant, the motion can be measured, in response to a known external force, . The displacement, the velocity, and the acceleration can be measured and calculated. At any particular time , (1) has only two unknowns, and . The equation describes a straight line in the plane. If the system is exactly described by (1), then the set of straight lines at successive times, , must all pass through a single point in the plane, which then determines the values of and .

If the system is nonlinear or departs from the assumptions of (1) in some other way, such as having nonviscous damping, the successive lines do not all pass through a single point. However, if the variation is reasonably smooth, they still define an envelope curve. This curve in the plane gives the time-dependent values of the equivalent linearized and . This paper presents an algorithm to map this envelope and illustrates the method with numerical and experimental examples.

2. Theoretical Background and Methodology

2.1. The Coordinate Plane

If  kg,  (N·sec/m), and  (N/m), then are the rectangular Cartesian coordinates of a point that has as the abscissa and as the ordinate, as shown in Figure 1.

Figure 1: The coordinate plane.

Since , for M.C.K. system values are both functions of time, , form a trajectory on the plane of coordinate. This study terms the characteristic trajectory of an M.C.K. system, as shown in Figure 2. Sections 5.2 and 5.3 show a simulation and experimental examination of this situation. Only and are constant in the time domain. The characteristic trajectory is shown as Figure 1.

Figure 2: and are both functions of time on the coordinate plane.

2.2. The Limits of the Critical Damping

For a critically damped system, the damping ratio is

and the stiffness is derived as

If the mass in Figure 2 is 1 kg, then the critical damping limit is shown by the dotted line in Figure 3. If  (N/m) and  (N·sec/m), then (3) gives  (sec) for underdamping,  (sec) for critical damping, and  (sec) for overdamping.

Figure 3: The critical damping limits for Figure 2, with  kg.

3. Methodology

The general movement formula for an M.C.K. system is

where is a mass, is a damping coefficient, is a stiffness coefficient, and is an external force. Equation (1) represents a two-dimensional equation () at time, , on the coordinate plane:

In (5), , are taken to be functions of time and are generally considered constants in the definition of the time domain. It is known that , , , and for M.C.K. systems, the displacement response, , the velocity, , and the acceleration, , at any time, , are necessary for (5); namely,

It is known that at time, , the displacement response of the system is . After differentiation, and are known as are and , but and are unknown. Equation (6) can be rewritten as the slope-intercept from an equation of the line:

At another time, , represents the linear equation of the M.C.K. motion formula on the plane coordinates:

Figure 4 shows that from times, to , and for an M.C.K. system have a linear binary equation, and , at the locations of the plane coordinates. In terms of the dynamics, the intercept of the line is driven by the displacement response, by the acceleration, and by external forces from to , and the slope is driven by the displacement response and the velocity from to . , for a linear binary equation for the motion equation of an M.C.K. system change the location of to , which results in and , respectively. If , for the system are constant in the time domain, that is, , then mark the intersection point with and , as shown in Figure 4. , for a linear binary equation for the motion equation of an M.C.K. system rotate on the plane about the rotational center .

Figure 4: The dynamic behavior of the slope-intercept from the equation of the line on the plane coordinates, from times, to .

If and for the system are both functions of time, , in the continuous time domain, it is assumed that

As , it is inevitable that and satisfies and , such that is an intersection point with and , as shown in Figure 5. is the instantaneous center of rotation on the plane coordinates for the motion formula for an M.C.K. system. The center also forms a trajectory that defines the coordinates of the characteristic trajectory. Equations (7) and (8) are rewritten as follows:

Figure 5: The dynamic behavior of the slope-intercept from the equation of the line on the plane coordinates. and for the system are both functions of time, , in the continuous time domain.

is solved as the simultaneous solution of (10), so the time-domain displacement response (the velocity and the acceleration calculated by differentiating the displacement) obtained by experiment can be used to solve and for an M.C.K. system. The entire method is described as follows:

As is seen, the values of and from the last formula cannot be determined. If necessary, the number of samples can be easily adjusted to satisfy the demands of the analysis.

In theory, this method is very sensitive to noise. Since the velocity and acceleration are obtained by differentiating the displacement, the displacement, the velocity, and the acceleration must be differentiable and continuous. In order to ensure the feasibility of this method, the measurement data must be smoothed.

4. Inference

For example, for the impact response of the experiment to momentum , the initial conditions are and .

If represents no damping, then , , and are given as follows:

Consequently, Figure 5 gives

Equation (14) defines the relationship between the stiffness coefficient, , and the intercept, , on the axis.

If , so the system is underdamped; then , , and are given as follows:

Consequently, Figure 5 gives

If , (16) and (17) are, respectively, equal to (13) and (14).

If , which means that the system is overdamped, then , , and are given as follows:

Consequently, Figure 8 gives

Although this is a linear system, the nonlinear time-variant systems observed within a small time interval are almost all linear, time-invariant systems. Using the reasoning for , , the relationship between the damping ratio, , and the behavior of the system on the plane is known. This determines the system damping ratio on the plane. In the special case of the absence of damping, , directly indicates the value on the plane. This means that an M.C.K. system can be directly observed on the plane ( is a fixed value).

5. Application

The simulation response of an M.C.K. system is the measurement data obtained from this experiment. This study inversely calculates the two adjacent data points of the simulation response to determine the and values for the system. In Section 5.3 an experiment involving an aerostatic bearing system is performed to validate this method.

5.1. Theoretical Analysis of Pulse Response Simulation

Assuming that the parameters of the MCK are  (kg),  (N·sec/m),  (N/m),  (kg·m/sec), and  (N), the initial conditions are  (m),  (m/sec), time = 7 s, the displacement response of 71 points is taken, and the time space is 0.1 s. Equation (15) gives solutions for the displacement, the velocity, and the acceleration, as shown in (20), (21), and (22), respectively. Consider

Figures 6(a), 6(b), and 6(c) show the displacement, the velocity, and the acceleration, respectively. Equation (11) gives the dynamic behavior of the slope-intercept from an equation of the line on the plane coordinates, as shown in Figure 6(d). Each line has an instantaneous center of rotation at the values and L1–L8 are used to label this portion of the data only. Figure 6(a) shows that the and values for this system occur before the first peak and that the value is solved in Figure 6(e). That is, the value is constant in the time domain. The value is solved in Figure 6(f). This means that the value in the time domain is also a constant. Table 1 shows the portion of the data, , , for time functions and the slope of the M.C.K. system in Figure 6, in which , , , and so forth.

Table 1: The portion of the data, , , from Figure 6(d).

Lt
(sec)
(m)
(m/sec)
(m/sec²)

(N/m)
 (N·sec/m)

L1
0.2
−0.0016
−0.0061
0.018612
11.59676
−3.79838
4
2

L2
0.4
−0.00247
−0.00269
0.015259
6.173171
−1.08659
4
2

L3
0.6
−0.00273

0.011026
4.036969
−0.01848
4
2

L4
0.8
−0.00255
0.001723
0.006754
2.64883
0.675585
4
2

L5
1
−0.0021
0.002687
0.003011
1.436505
1.281747
4
2

L6
1.2
−0.00152
0.002984
0.000111
0.072932
1.963534
4
2

L7
1.4
−0.00094
0.002795
−0.00185
−1.97586
2.987929
4
2

L8
1.6
−0.00042
0.002304
−0.00292
−6.92301
5.461503
4
2

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Using the time-domain displacement response from the Pulse Response Simulation to determine the , values. (a) The time-domain displacement response is . (b) The velocity response is , which is the differential of . (c) The acceleration response is , which is the differential of . (d) The characteristic trajectory of the motion formula on the plane. (e) The coefficient of stiffness is solved as graph . (f) The coefficient of viscous damping is solvedas graph .

Figure 7(a) shows that and are solved by inverse calculation. That is, and for this system are constant, similar to the parameter settings for the response simulation. The , coordinate is , which is on the upper left of the critical damping limits, so this is an underdamped system. The calculated margin of error is less than , so it can be ignored. Figures 7(b), 7(c), and 7(d) compare and with and . In the plot, the last item is neglected because it cannot be determined at the end of the , value.

(a)

(b)

(c)

(d)

Figure 7: The theoretical analysis of the Pulse Response Simulation at the coordinate plane. (a) The characteristic trajectory of the motion formula and the critical damping limit for the M.C.K. system on the plane. (b) The error in the values distributed on the coordinate plane. (c) The time-domain error in . (d) The time-domain error in .

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8: Using the time-domain displacement response for the simulation experiment to solve the , values, which are both functions of time. (a) The time-domain displacement response is . (b) The velocity response is , which is the numerical differential of . (c) The acceleration response is , which is the numerical differential of . (d) The characteristic trajectory of the motion formula on the plane. (e) The coefficient of stiffness is solved as graph, , , and is not constant. (f) The coefficient of viscous damping is solvedas graph, , , and is not constant.

5.2. Numerical Analysis of a Simulated Pulse Response

In order to confirm the feasibility of this method in identifying nonlinear systems, it is assumed that a curve imitates an impact test process, as shown in Figure 8(a). This virtual curve must represent the trajectories of a nonlinear system, in order to verify the accuracy of the proposed method. For these reasons, a displacement theoretical equation for the virtual nonlinear system is presented in (23). The displacement trajectories for the virtual nonlinear systems are obtained from (23) and are shown in Figure 8(a). Using the finite difference method to obtain the speed in Figure 8(b) from the displacement curve in Figure 8(a), and the acceleration in Figure 8(c) from the speed in Figure 8(b). The methods already mentioned are then used to obtain the parameters, and . Using (23), (24), and (25), another set of theoretical solutions for the displacement, the velocity, and the acceleration are then calculated. The theoretical finite difference solutions are compared with the theoretical solutions to verify the accuracy of the method. Because this is an envelope approximation, the smaller the time intervals, the more accurate the solutions. The accuracy of the calculation process depends on the number of iterations.

Assuming that (kg) and , are unknown, the M.C.K. system receives a single pulse of force, with no other external forces. The simulation of the displacement response, defined by (23), is shown in Figure 8(a):

The simulation takes 4 s and the time space is 0.001 s. Although (23) can be deduced as the velocity of (24) and as the acceleration of (25), the displacement response is used by the numerical differential to obtain the velocity and the acceleration during the actual experiment, as shown in Figures 8(b) and 8(c). The initial conditions are obtained from the displacement response and the velocity, so the size of the pulse can be ignored. The status of other forces, which is (N), is the only other consideration. Figure 8(d) shows that (11) solves the , values as a characteristic trajectory on the plane. Table 2 shows the portion of the , values that represent time functions and the slope of the M.C.K. system in Figure 8, in which , , . Figures 8(e) and 8(f) show that the and values can change in the time domain. In the plot, the last value for , cannot be directly solved. However, by using the Central Difference Method, the velocity can be determined on the first application. The acceleration can be calculated from this velocity data, using the same method, so the first two and last three items are not considered. When the two adjacent sampling points have been calculated, the system can determine instantaneous , values. Data 2 in Figure 9(a) shows that (23) derives the velocity for (24) and the acceleration for (25) by using 4001 points to determine the , values of the system, and data 1 is obtained by the numerical differentiation of the velocity and the acceleration. Figures 9(b), 9(c), and 9(d) show that the error is reduced after  s of the initial response, which is seen by comparing data 1 with data 2. Data 3 in Figure 9(a) shows that  (kg) at the critical damping limit for the M.C.K. system. The characteristic trajectory of , increases with time and approaches the critical damping limits. That is, the damping ratio increases with time and approaches the critical damping ratio and the coordinate of the characteristic trajectory for , in the upper left boundary, so the M.C.K. system is underdamped at the time of the simulation. This is similar to the design of a shock absorber system that includes variable stiffness and damping. Firstly, the response of the system is determined and then the differentiation method previously discussed is used to determine the necessary stiffness and damping.

Table 2: The portion of the data, , for the time functions from Figure 8(d).

Lt
(sec)
(m)
(m/sec)
(m/sec²)

(N/m)
(N·sec/m)

L1
0.01
−0.05955
−5.89536
13.58601
228.1615
−99.0059
70.05914
1.596898

L2
0.03
−0.17433
−5.5664
18.55715
106.4511
−31.9311
35.84346
2.211251

L3
0.09
−0.47109
−4.28554
22.77945
48.35459
−9.09703
20.07025
3.109183

L4
0.19
−0.78703
−2.07592
20.35073
25.85776
−2.63768
15.25057
4.021417

L5
0.33
−0.90441
0.198556
11.87931
13.1349
0.219543
14.23512
5.011415

L6
0.51
−0.73229
1.434002
2.598733
3.548783
1.958248
15.49573
6.100836

L7
0.73
−0.40647
1.351987
−2.25806
−5.55523
3.326129
18.51146
7.235646

L8
0.99
−0.14247
0.67754
−2.37718
−16.6854
4.755633
22.61141
8.263207

(a)

(b)

(c)

(d)

Figure 9: The simulated pulse response on the plane. (a) The characteristic trajectory of the motion formula and the critical damping limit for the M.C.K. system on the plane. (b) The error in the , values on the coordinate plane. (c) The time-domain error in . (d)The time-domain error in .

5.3. An Aerostatic Bearing System

During the transient response interval, an aerostatic bearing has nonlinear stiffness and damping characteristics. Therefore an experiment to determine the dynamic parameters for an aerostatic bearing is proposed. The aerostatic bearing is a closed loop feedback system. When it is located at a balance point, the gas film gap remains constant. When the load changes, the aerostatic bearing air gap has a new equilibrium position.

In general, an aerostatic bearing’s dynamic performance depends on the dynamic parameters of the bearing, such as the stiffness and damping values. Since finding an exact solution for these values involves complex fluid mechanics problems, it is difficult to determine them. In order to verify the nonlinear dynamic parameter identification method presented in this paper, an experiment to measure stiffness and damping values of the aerostatic bearing was performed. Figures 10 and 11 show a schematic diagram and a photograph of the experimental system. A free-falling steel ball is the pulse force, which acts on the aerostatic bearing and excites transient responses. Using the transient responses and the presented method, the bearing dynamic characteristics, such as the stiffness and damping, can be determined. The experimental results are shown in Figure 12.

Figure 10: The impact response experiment for the aerostatic bearing system.

Figure 11: A photograph of the impact response experiment.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 12: The , values for the aerostatic bearing system on the plane.

For the experimental results of the displacement, the smoothing process is performed using the MATLAB software function, “polyfit,” as shown by the red line in Figure 12(a). On the time axis, the smooth curve deviates from the actual measured curve after 0.025 seconds, because there is insufficient data. As it is close to the steady state, it does not affect the determination of the transient response.

Figure 12(a) section [A, E] shows a complete waveform, which marks the major part of the transient response of this system. In section [A, E] the fitted curve and the actual measured values are in good agreement and this is the interval of interest. Intervals [A, B] and [D, E] are the range of the gas film squeeze process and interval [B, D] is the range of the gas film gap expansion process. Points A, C, and E are the aerostatic bearing equilibrium positions where the displacement is zero.

Figure 12(b) shows that interval [A, E] is a continuous curve whose track is formed by an envelope. Figures 12(c) and 12(d) show the change in the values of and in the time domain. Intervals [A, B] and [B, D] are smooth curves, which show that the change in is inversely proportional to the air gap shown in Figure 12(e). In these sections, the stiffness of the aerostatic bearing decreases. At point B, the air gap has the smallest value, but the stiffness, , has a maximum value. However, at point C, the air gap has the largest value, but the stiffness is a minimum.

The curve at interval [A, E], shown in Figure 12(f), shows the relationship between and the air gap, which is continuous and nonlinear. Interval [B, D] shows the rising state of the aerostatic bearing. During interval [B, C] it is seen that is inversely proportional to the absolute value of the air gap. However, during the interval [C, D] it is seen that is proportional to the absolute value of the air gap. Intervals [A, B] and [D, E] show the relationship between the air gas and during compression of the bearing.

6. Conclusions

This paper presents a novel and effective method for the identification of nonlinear systems, using plane coordinates. A new time-domain analysis is proposed. In the plane coordinates, the dynamic behavior of an M.C.K. system provides useful information. In plane coordinates, any motion formula for an M.C.K. system can be expressed as a line in a geometric plane (a two-dimensional space) and each line has an instantaneous center of rotation on the coordinate plane in this system. The two adjacent sampling data are calculated and then the values of and for the instantaneous system are determined. and are both functions of time, so in a traditional M.C.K. system, and are fixed constants, but there may be exceptions. In other words, the design that allows variable stiffness and damping can be identified in the time domain.

Nomenclature

| Coefficient of viscous damping
| Coefficient of stiffness
| Mass
| Rectilinear equations
| A single pulse force (impact force)
| Force functions of time
| Time
| Displacement of the M.C.K. system
| Velocity
| Acceleration
| The slope of the line ()
| The intercept of the line ().

Greek Symbols

| Damping ratio
| Natural frequency
| Damped natural frequency.

Superscripts

Right convergence.

Subscripts

| Node time
| Functions of time
| Node time.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

S. Adhikari and J. Woodhouse, “Identification of damping. Part 1: viscous damping,” Journal of Sound and Vibration, vol. 243, no. 1, pp. 43–61, 2001.

A. Srikantha Phani and J. Woodhouse, “Viscous damping identification in linear vibration,” Journal of Sound and Vibration, vol. 303, no. 3–5, pp. 475–500, 2007.

S. R. Ibrahim and E. C. Mikulcik, “A method for the direct identification of vibration parameters from the free response,” Shocked Vibration Bulletin, vol. 47, pp. 183–198, 1977.

S. R. Ibrahim, “An approach for reducing computational requirements in modal identification,” AIAA Journal, vol. 24, no. 10, pp. 1725–1727, 1986.

S. R. Ibrahim, “Random decrement technique for modal identification of structures,” Journal of Spacecraft and Rockets, vol. 14, no. 11, pp. 696–700, 1977.

M. E. Gaylard, “Identification of proportional and other sorts of damping matrices using a weighted response-integral method,” Mechanical Systems and Signal Processing, vol. 15, no. 2, pp. 245–256, 2001.

C. R. Zhou and C. S. Zhao, Mechanical Vibration Parameter Identification and Its Application, Beijing Science Press, Beijing, China, 1989 (Chinese).

C.-H. Lamarque, S. Pernot, and A. Cuer, “Damping identification in multi-degree-of-freedom systems via a wavelet-logarithmic decrement. Part 1: theory,” Journal of Sound and Vibration, vol. 235, no. 3, pp. 361–374, 2000.

W. J. Staszewski, “Identification of damping in MDOF systems using time-scale decomposition,” Journal of Sound and Vibration, vol. 203, no. 2, pp. 283–305, 1997.

J. Chen, Y. L. Xu, and R. C. Zhang, “Modal parameter identification of Tsing Ma suspension bridge under Typhoon Victor: EMD-HT method,” Journal of Wind Engineering and Industrial Aerodynamics, vol. 92, no. 10, pp. 805–827, 2004.

F.-L. Huang, X.-M. Wang, Z.-Q. Chen, X.-H. He, and Y.-Q. Ni, “A new approach to identification of structural damping ratios,” Journal of Sound and Vibration, vol. 303, no. 1-2, pp. 144–153, 2007.