- Email: [email protected]
A method for simultaneous identification of the full nonlinear damping and the phase shift and amplitude of the external harmonic excitation in a forced nonlinear oscillator
Computers and Structures 120 (2013) 77–85
Contents lists available at SciVerse ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
A method for simultaneous identification of the full nonlinear damping and the phase shift and amplitude of the external harmonic excitation in a forced nonlinear oscillator T.S. Jang ⇑ Naval Architecture and Ocean Engineering, Pusan National University, Busan 609-735, Republic of Korea
a r t i c l e
i n f o
Article history: Received 4 October 2012 Accepted 14 February 2013 Available online 15 March 2013 Keywords: Simultaneous identification Non-parametric identification Full nonlinear damping Phase shift and amplitude Lack of stability Regularization
a b s t r a c t We propose a new general method that is able to clarify the simultaneous identification of the full nonlinear damping characteristics as well as the phase shift and amplitude of the external harmonic excitation in a forced nonlinear single degree system. The method is the non-parametric identification based on the concept of zero-crossings regarding the system responses. It involves an integral equation of the first kind, resulting in the lack of stability properties. The difficulty of instability is overcome by the stabilization technique. The workability of the proposed method is demonstrated through a numerical experiment consisting of a highly nonlinear system. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Associated with the mathematical modeling of dynamical systems, a considerable amount of attention has been given to the system identification of nonlinear physical systems in various branches of science and engineering. In fact, an appropriate and correct mathematical modeling from measured response data through identification plays an essential role in achieving accurate predictions of motion responses to various loading environments, which is also important for improving the system performance such as the effective control and the estimation of the structure lifetime. Especially, the damping identification is extremely crucial to the forced systems in the field of structural engineering, mainly because it can be directly applied to the early-stage design of nonlinear physical systems. A number of methods have been proposed in the literature sourced for the damping identification of nonlinear systems arising in mechanics, whose primary objective involves the formulation of a mathematical model on the basis of measured data. Iourtchenko et al. [1] successfully presented a novel technique for an identification scheme in the case of parametric excitation. Iourtchenko et al. [2] further proposed a mathematical procedure based on the method of stochastic
⇑ Tel.: +82 51 510 2789; fax: +82 51 581 3718. E-mail address: [email protected] 0045-7949/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruc.2013.02.008
averaging for in-service identification of the damping characteristic from measured random vibration. Recently, Kang et al. [3] estimated the stiffness and damping parameters of a structure by using the measured acceleration. Lee and Law [4] proposed a procedure for the viscous damping model adopting the inverse modal transformation. With regard to the modeling of nonlinear damping systems, there were also the challenging issues on other studies [5–8]. What this paper involves is the simultaneous identification of nonlinear physical systems, i.e., a study to identify more than a ‘‘single’’ physical element of a system. This is in contrast to the usual current identification. In fact, most related studies on the system identification, however, are primarily limited to the identification of a ‘‘single’’ physical element of a system, say, either a system’s damping force or restoring force. Even though only few attempts have so far been made at the simultaneous identification, a recent work may be found in the paper of Jang et al. [9], where they newly proposed a general procedure for identifying both a full nonlinear restoring force and a full nonlinear damping function in _ in which D(y) is a general damping function of the form, DðyÞ y, the displacement y and y_ is the velocity in nonlinear oscillation systems. In addition, Jang [10] modified the general procedure _ as well [9] to identify a functional form of nonlinear damping, DðyÞ, as a full nonlinear restoring force. In this paper, we propose a new mathematical procedure to simultaneously measure not only the full nonlinear damping in _ but the phase shift and amplitude of external the form, DðyÞ y; harmonic excitation in a forced nonlinear single degree system.
78
T.S. Jang / Computers and Structures 120 (2013) 77–85
The procedure proposed here is thus regarded as a simultaneous identification of two physical elements of a system. However, it is different from the methods of the simultaneous identification [9,10] in that the two physical elements to be identified in the present study are the damping and the external force of a system. Of course, the procedure is also distinct from the (usual) force identification of a single physical element of the external force of a system. A method for force identification is, for example, illustrated in Jang et al. [11], who proposed a new method to find nonharmonic periodic excitation forces in nonlinear damped systems. When we often encounter vibrations due to the harmonic excitation, it may be hard to perceive the phase shift and amplitude of external harmonic loading beforehand [12]. Hence, in practice, the present work may be the one that has a wider applicable range than the existing methods identifying only the single damping element of a system. Furthermore, the procedure proposed is a fully nonlinear method. Consequently, it does not depend on a small parameter and thus overcome the disadvantage and limitations of the parametric methods such as a perturbation expansion procedure. Finally, it is found that a numerical experiment demonstrates the workability of the proposed method of the measuring procedure in nonlinear oscillations by using a forced van der Pol equation with cubic nonlinear restoring force. 2. The equation of motion We consider a nonlinear damped system, subject to the sum of two external harmonic loads, asin(xt) and bsin(xt), where the positive constants a, b and x are amplitudes and frequency. System dynamics considered are governed by Newton’s equation of motion, expressed as in the following form
€ þ DðyÞy_ þ RðyÞ ¼ a sinðxtÞ þ b cosðxtÞ for t > 0: my
ð1Þ
Mathematically, Eq. (1) is a second-order nonlinear ordinary differential equation, where the constant m > 0 represents the mass of a particle of an oscillator, the term involving D(y) the general form of nonlinear damping and the term R(q) a general nonlinear restoring force. The harmonic excitation in Eq. (1) can be, alternatively, represented in terms of the phase shift a and the amplitude r as follows
a sinðxtÞ þ b cosðxtÞ ¼ r sinðxt þ aÞ
ð2Þ
where
r¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a2 þ b
3.1. Pseudo-spring constant We attempt to develop a normal differential form for Eq. (1), which plays an essential role in deriving an integral equation in this study. To this end, we first introduce a pseudo-spring constant, kp > 0, and add the linear spring force of the form kpy to both sides of Eq. (1) to have
€ þ kp y ¼ fDðyÞy_ RðyÞ þ f ðtÞg þ kp y; my
t>0
ð5Þ
where
f ðtÞ a sinðxtÞ þ b cosðxtÞ:
ð6Þ
Eq. (5) is then simply rewritten as
€ þ kp y ¼ FðtÞ my
ð7Þ
with
FðtÞ DðyÞy_ RðyÞ þ f ðtÞ þ kp y:
ð8Þ
Eq. (7) is viewed as a normal (differential) form associated with Eq. (1) in this paper, which corresponds to an equivalent integral equation, as will be discussed in Section 3.2. 3.2. Volterra integral transform For simplicity, we assume that the motion is initially at rest for time t < 0 such that
_ yð0Þ ¼ 0:
yð0Þ ¼ 0;
ð9Þ
The displacement y(t) in Eq. (7) is then expressed as the Volterra integral transform of F in Eq. (8) with the kernel G(t, s:kp) [10]
yðtÞ ¼
Z
t
Gðt; s : kp ÞFðsÞds;
ð10Þ
0
and
Gðt; s : kp Þ ¼
1 y1 ðsÞy2 ðtÞ y1 ðtÞy2 ðsÞ : m y1 ðsÞy_ 2 ðsÞ y_ 1 ðsÞy2 ðsÞ
ð11Þ
ð3Þ
In Eq. (11), y1(t) and y2(t) are the motion responses of the linear oscillator with the pseudo-spring constant kp:
ð4Þ
€1 ðtÞ þ kp y1 ðtÞ ¼ 0; my €2 ðtÞ þ kp y2 ðtÞ ¼ 0; my
and
tan a ¼ b=a:
study: Given a motion response of a forced nonlinear system, is it possible to identify the nonlinear damping of the system and the phase shift and amplitude of harmonic excitation at the same time? To answer the question, we begin with an integral equation.
As mentioned briefly in introduction, the correct and accurate measurement of nonlinear damping of a system plays an important role in applications in a variety of physical problems in different disciplines. However, as far as a forced system is concerned, it is practically difficult to know the phase shift a in Eq. (4) and amplitude r in Eq. (3) of external loading directly from experiment in advance. In fact, it is not easy to install transducers within a moving structure during its construction, so that there are only few cases in which we can directly determine external loading in forced oscillations [11]. To solve the difficulty, a new general method will be proposed for measuring both the full nonlinear damping and the phase shift a and amplitude r of external harmonic load during forced nonlinear oscillations. 3. Integral equation The present study may be considered as an inverse problem in the form of a question, which is the basic idea underlying this
y1 ð0Þ ¼ 1; y_ 1 ð0Þ ¼ 0 y2 ð0Þ ¼ 0; y_ 2 ð0Þ ¼ 1:
ð12Þ
In general, the identification of nonlinear systems relies on measuring system’s responses. In this work, we assume that the displace_ ment y(t) and velocity yðtÞ of the system are measured during oscillations. Mathematically, this implies that the transform of Eq. (10) may also be regarded as an integral equation for F(t) in Eq. (8), because the left hand side of Eq. (10) is known from the measurement. Regarding the integral equation, it is meaningful to check whether Eq. (10) has a unique solution F(t), which was discussed and proved in Jang et al. [9]. It is noted that if Eq. (8) is rewritten as
f ðtÞ DðyÞy_ ¼ FðtÞ þ RðyÞ kp y;
ð13Þ
the right hand side of Eq. (13) becomes then known, since the (unique) F can be determined by solving the integral equation of Eq. (10). Denoting the known by u(t),
uðtÞ FðtÞ þ RðyÞ kp y; we have an alternative expression for (13), namely
ð14Þ
T.S. Jang / Computers and Structures 120 (2013) 77–85
_ uðtÞ ¼ f ðtÞ DðyÞy:
ð15Þ
4. Simultaneously measuring procedure
79
Then, the value of the nonlinear damping function D, evaluated at y = y(ti), can be explicitly calculated as
Dðyðt i ÞÞ ¼
f ðt i Þ uðti Þ _ iÞ yðt
ð24Þ
In this section, we present a mathematical procedure for the simultaneous measuring of both the nonlinear damping and the phase shift and amplitude of external excitation, the main idea of which is provided with aid of Eq. (15).
according to Eq. (15). To identify D(y) in more points of y, the above process can be repeated as many times as necessary for each ti, resulting in the desired measurement of the unknown D.
4.1. Phase shift and amplitude of harmonic excitation
5. Lack of stability
Notice that the external force can be balanced by the (known) u(t), whenever the damping term in Eq. (15) disappears. Actually, _ this happens when the velocity yðtÞ of the system response vanishes during oscillations. For a systematic analysis, let us define a _ zero-crossing time t = t1, t2,. . . in terms of the velocity yðtÞ of the system response, that is, it satisfies
Even though we have constructed a measuring procedure for the phase shift and amplitude of harmonic excitation and the functional form of nonlinear damping in the previous section, it is still necessary to discover u(t) in Eq. (15), equivalently, F(t) in Eq. (10). However, Eq. (10) is classified as an integral equation of the ‘‘first’’ kind, which leaves the question of stability or instability open.
_ i Þ ¼ 0 for 0 < t1 < t2 < . . . : yðt
5.1. Hilbert–Schmidt integral operator K
ð16Þ
Then, the damping force becomes zero
_ iÞ ¼ 0 Dðyðt i ÞÞ yðt
ð17Þ
when t = t1, t2,. . .. This immediately yields
uðti Þ ¼ f ðt i Þ;
t ¼ t 1 ; t2;...
y ¼ KðFÞ; ð18Þ
from Eq. (15).If we choose two first zero-crossing times in terms _ of the velocity yðtÞ, say, t = t1, t2, then Eq. (18) turns out to be a linear matrix equation for the unknown a and b appearing in Eq. (6):
uðt 1 Þ uðt 2 Þ
¼M
a
ð19Þ
b
where M denotes the two by two matrix
M¼
ð25Þ
where K is defined as the integral operator
KðFÞ
Z
t
Gðt; s : kp ÞFðsÞds:
ð20Þ
Assuming that the matrix M in Eq. (20) is nonsingular, the two amplitudes of a and b can be calculated by using the inverse M1:
ð26Þ
0
In what follows, we shall investigate the topological structure of the operator K in Eq. (26). We start by noticing that the responses of y1(t) and y2(t) are harmonic functions and the denominator of the kernel G in Eq. (11) is a Wronskian which is a non-zero constant [9]:
W ¼ y1 ðtÞy_ 2 ðtÞ y_ 1 ðtÞy2 ðtÞ ¼ const:
sinðxt1 Þ; cosðxt1 Þ : sinðxt2 Þ; cosðxt2 Þ
For convenience, Eq. (10) can be rewritten as the operator equation in a compact way
ð27Þ
This immediately leads to the inequality with regard to the kernel G
Z 0
T
Z
T
G2 ðt; s : kp Þdsdt < 1
ð28Þ
0
The amplitudes, a and b, so calculated, in turn, enable the phase shift a and the amplitude r in Eqs. (3) and (4) to be completely determined.
for an arbitrary (large) constant T > 0. We recall that a kernel satisfying the inequality in Eq. (28) is called a Hilbert–Schmidt kernel and the integral operator with a Hilbert–Schmidt kernel is known as a Hilbert–Schmidt integral operator [13]. Thus, it is found that the kernel G in Eq. (11) and the operator K in Eq. (26) are a Hilbert–Schmidt kernel and a Hilbert–Schmidt integral operator, respectively.
4.2. Nonlinear damping
5.2. Unbounded operator K1: discontinuity
Since the forcing f(t) in Eq. (6) is shown to be determined through Eq. (21) in Section 4.1, it is possible to estimate the following ratio
According to functional analysis, a Hilbert–Schmidt integral operator is shown to be a compact operator, that is, it transforms bounded sets into compact sets. Furthermore, a compact operator, defined on infinite dimensional spaces such as function spaces, is continuous, while its inverse is not continuous. Therefore, the operator K in Eq. (26) is compact and thus K1 is a discontinuous operator (or unbounded operator). Accordingly, F(t) in Eq. (10), represented in terms of the displacement y(t) as the below equation, does not depend continuously on the y(t) in Eq. (7), lacking stability properties of solutions [14,15]:
a b
¼ M 1
uðt 1 Þ
uðt 2 Þ
f ðtÞ uðtÞ ; _ yðtÞ
ð21Þ
ð22Þ
_ provided that yðtÞ – 0. In addition, we note that, from Eq. (15), the ratio is identical to the damping function D(y) that we are looking for. Thus, the estimation of the ratio implies the measurement of D(y), which is accomplished by using a proper time with regard to system responses. In contrast to a zero-crossing time t = t1, t2,. . . as in Section 4.1, we shall select a non-zero-crossing time in terms of the velocity y_ of the system response during oscillations, denoted as t = t1, t2, . . ., such that
_ i Þ – 0 for 0 < t 1 < t 2 < . . . : yðt
ð23Þ
F ¼ K 1 y
ð29Þ
Due to the discontinuity of K1, a small perturbation in y(t) or a small noisy error can be extremely amplified so that unreliable solutions may result from a small amount of noisy data. In fact, although Eq. (10) may be directly discretized into a numerical system of simultaneous linear equations (for example, by using the
80
T.S. Jang / Computers and Structures 120 (2013) 77–85
Fig. 1. Two system responses of the dynamic oscillatory system with the initial conditions of Eq. (35): c = 4.
Fig. 3. A typical graphical illustration of the L-curve: c = 4 and d = 0.01. Fig. 2. Phase diagram corresponding to the dynamic responses in Fig. 1: c = 4.
Table 1 Recovered a and b (or, r and a) of excitations for different noise level d: c = 4. Recovered
Exact d = 0.01 d = 0.05 d = 0.10 d = 0.15
a
b
4 3.7667 4.4980 3.5749 4.5311
8 8.0082 7.8518 7.5785 8.8896
r
a
Average error (%)
8.9443 8.8442 9.0489 8.3791 9.9778
1.1071 1.1308 1.0506 1.1300 1.0994
– 2.9 6.8 7.9 12.2
trapezoidal numerical integration rule), the directly discretized numerical system turns out to be ill-conditioned [14].
we thus would like an approximation for the solution to depend continuously on y(t). Tikhonov’s regularization is the most commonly used method of regularization [15]. In this study, however, we suggest the use of the Landweber’s regularization, which is known as a simple but effective iterative method. Thereby, the solution F in Eq. (10) is obtained as a limit of the sequence {Fj(t)} in the following iterative form of operator equation
F j ðtÞ ¼ F j1 ðtÞ kL fLðF j1 ðtÞÞg þ kL ðyðtÞÞ;
In order to prevent the lack of stability from affecting the performance of the present study of simultaneous measurement, we will suppress the instability through a stabilization (or regularization) technique. The main idea of the method of regularization is to replace an (unstable) integral equation of the ‘‘first’’ kind by a nearby (stable) integral equation of the ‘‘second’’ kind. Mathematically,
ð30Þ
for a real positive constant, k > 0. In Eq. (30), the operator L is the Volterra integral transform in the form in Eq. (10)
LðzÞ 5.3. Landweber’s regularization: a stabilization technique
j ¼ 1; 2 . . .
Z
t
Gðt; s : kp ÞzðsÞds
ð31Þ
0
for an arbitrary function z(t). The symbol L⁄ stands for the adjoint operator of L, and I the identity operator [15]. Based on the fact that Landweber’s regularization converges to the solution for an arbitrary initial guess F0(t), we start in this study with the initial guess of zero function, i.e., F0(t) 0, without loss of generality [15]. According to the theory of regularizations, there exist an optimal number of iterations in Eq. (30) for finding an accurate
T.S. Jang / Computers and Structures 120 (2013) 77–85
81
Fig. 4. Convergence behaviors of Fj(t): c = 4 and d = 0.01 (a) j = 104, (b) j = 105, (c) j = 106 and (d) j = 108.
and stable solution in the presence of noisy data. We shall apply Lcurve criterion [16] to determine the optimal number of iterations in Eq. (30), which will be discussed in more detail in Section 6. 6. Numerical experiment In this section, the workability of the proposed method will be demonstrated for measuring simultaneously the full nonlinear damping and the phase shift and amplitude of the external harmonic excitation during forced nonlinear oscillations. To this end, a model equation is chosen for the numerical experiment. 6.1. Forced Van der Pol equation
Fig. 5. Recovered f(t) asin(xt) + bcos(xt) compared with the exact one: c = 4 and d = 0.01.
We introduce a forced Van der Pol equation having a cubic nonlinear restoring to the numerical experiment, which evolves in time according to the second order differential equation.
€ þ cðy2 1Þy_ þ k1 y þ k2 y3 ¼ a sinðxtÞ þ b cosðxtÞ: my
ð32Þ
In the case of k2 = a = b = 0, Eq. (32) reduce to the Van der Pol equation which has a long history of being used in physical sciences, as non-conservative systems with nonlinear damping. The immediate comparison of equations (1) and (32) reveals the functional forms of D(y) and R(y) in Eq. (1), represented as
DðyÞ ¼ cðy2 1Þ and RðyÞ ¼ k1 y þ k2 y3 :
ð33Þ
For convenience and without loss of generality, the mass m in Eq. (32) is normalized as unit. The values for the spring constants k1 and k2 are taken as 50 and 5, respectively, and the amplitudes of excitation are set to be
a ¼ 4 and b ¼ 8:
Fig. 6. Recovered D(y): c = 4 and d = 0.01.
ð34Þ
We further select the interval of time domain as 0 < t < 1.2. Figs. 1 and 2 depict the time evolution solutions for Eq. (32) and their phase diagrams, respectively: Runge–Kutta type integration
82
T.S. Jang / Computers and Structures 120 (2013) 77–85
Fig. 7. Recovered f(t) asin(xt) + bcos(xt) compared with the exact ones for four different various frequencies: c = 4 and d = 0.01.
Fig. 8. Recovered D(y) compared with the exact one for four different various frequencies: c = 4 and d = 0.01.
schemes are used as numerical integration methods, in which the value for the coefficient c is taken as 4 and the initial condition is imposed by
yð0Þ ¼ 0;
_ yð0Þ ¼ 0:
ð35Þ
With data depicted in Figs. 1 and 2, we are now ready to follow the (inverse) measuring procedure proposed in the previous sections. By using the data of Figs. 1 and 2, we will attempt to identify the damping and forcing of a forced Van der Pol equation in Eq. (32) and the results are to be compared with the corresponding exact expressions, D(y) = c(y2 1), in Eq. (33) and f(t) in Eq. (6) with the coefficients of Eq. (34), respectively.
Fig. 9. Recovered f(t) asin(xt) + bcos(xt) and D(y) for 3 cases of pseudo-spring constant kp = 1, 3 and 5.
T.S. Jang / Computers and Structures 120 (2013) 77–85
83
Fig. 10. Recovered f(t) asin(xt) + bcos(xt) and D(y) for four different noise levels when jopt = 106: c = 4.
6.2. Noisy data and L-curve criterion As mentioned earlier, for the present simultaneous measurement, we are required to measure the displacement and velocity data of motion responses beforehand. In practice, however, the measured data are always deteriorated to some extent by noisy error. This implies that the left-hand side in the integral Eq. (10) is never exactly known but only up to an error of, say, noise level
d > 0. We thus assume that we know, from measuring devices, d > 0 and noisy data yd(t) such that
kyðtÞ yd ðtÞk2 6 d;
ð36Þ
in which the symbol kk2 denotes L2 norm [9–11,17]. In this study, four cases of noise levels of noisy data are randomly generated, which are tabulated in Table 1. Generally, in
84
T.S. Jang / Computers and Structures 120 (2013) 77–85
Table 2 Recovered a and b (or, r and a) of excitations for different damping coefficients c when d = 0.01. c
Case 1 Case 2 Case 3
2 3 4
Recovered a
b
4.426 4.037 3.767
7.623 8.292 8.002
r
a
Average error (%)
8.814 9.222 8.442
1.045 1.118 1.130
7.7 2.3 2.9
the presence of noisy data, Landweber’s regularization approaches the true solution in the initial stage of iteration, but it potentially deviate far away at a large number of iterations. As a result, the iterated solution of Landweber’s regularization goes to a useless solution because of the amplification of hidden noise. Depending on a given noise level d > 0, the accuracy of the Landweber’s regularization can be influenced by the number of iterations. One of the most effective ways to find the optimal number of iterations may be of an L-curve criterion [16], which mainly concerns how to select an appropriate iterated solution. A brief sketch is as follows. If we consider the two quantities
rj ¼ LðF j ðtÞÞ yd ðtÞ2
and sj ¼ F j ðtÞ2
j ¼ 1; 2 . . .
ð37Þ
which correspond to the norms of the residual and the iterated solution, respectively, then, the residual rj tends to decrease, while the value of sj increases as the number of iterations increase. This can be best illustrated by introducing a plot of the pair of (rj, sj), which may yield a curve exhibiting a typical ‘‘L’’ shape. According to L-curve criterion [16], an optimal number of iterations are considered to be the one corresponding to the corner of an L-shaped curve. Fig. 3 shows a typical shape of the L-curve in the present study, where the convergence constant k > 0 is chosen as k ¼ 1: here, c = 4 and d = 0.01. In relating the calculation of Eq. (31), essential for the Landweber’s iteration of Eq. (30), the pseudo-spring constant is set to be unit in this section, i.e., kp = 1. The optimal iteration number jopt appears to be clear at the corner of the L-curves as shown in Fig. 3, which is jopt = 106. The convergence behaviors of iterative solutions Fj(t) in the iteration Eq. (30) are illustrated in Fig. 4. Here, the Fj(t) are compared with the exact one, F(t) in Eq. (8), corresponding to Eq. (33). We see that not-enough iteration
Fig. 11. Recovered f(t) asin(xt) + bcos(xt) and D(y) for three different coefficients of c when jopt = 106.
T.S. Jang / Computers and Structures 120 (2013) 77–85
produces a poor but stable approximation, while too many iterations result in a worse (unstable) solution. 6.3. Measuring D(y), a and r Following the simultaneous measuring procedure as proposed in Section 4, we try to calculate the first two zero-crossing times t1, t2 in terms of the velocity y_ from data in Fig. 1. The results are as follows
t 1 ¼ 0:46 s; t 2 ¼ 0:91 s:
ð38Þ
Using the two times t1, t2, we can determine the amplitudes of the two unknowns a and b in Eq. (6) by solving the matrix equation of Eq. (21) and thus find out the phase shift a and the amplitude r from Eqs. (3) and (4). The recovered ones are a = 3.7667, b = 8.0082, equivalently, a = 1.1308rad and r = 8.8442, which have the average error about 3%. Fig. 5 depicts the recovered external harmonic excitation f(t) asin(xt) + bcos(xt) in Eq. (6) compared to the exact one, which is f(t) in Eq. (6) by using Eq. (34). Next, we select 25 non-zero-crossing times of Eq. (23) in terms _ of the velocity yðtÞ of the system response during oscillation. With the selected non-zero-crossing times, we apply Eq. (24) to recover the nonlinear damping function D(y). The results are shown in Fig. 6, where it is relatively accurate when compared with the exact results D(y) = c(y2 1) of Eq. (33). Figs. 7 and 8 also depict the identified f(t) in Eq. (6) and D(y) for the optimal iteration number jopt = 106, respectively. However, four different various frequencies are used for the results, which are in a good agreement with the exact ones of f(t) in Eq. (6) with Eq. (34) and D(y) = c(y2 1) of Eq. (33), respectively. We notice that it is clear from Eq. (5) that D(y) and f(t) are invariant with respect to the pseudo-spring constant kp, that is, they do not depend on a choice of kp. This is also numerically confirmed as shown in Fig. 9, where D(y) and f(t) are calculated based on the choices of the three different kp = 1, 3 and 5. It may be hard to see any effect of a pseudo-spring constant kp on D(y) and f(t), as expected. In Fig. 10, we examine f(t) and D(y) according to the four different noise levels d as listed in Table 1. Finally, we pick up three cases of different coefficients for c in Table 2, to check their dependency on f(t) and D(y). As exhibited Fig. 11, the results are compared with exact values of f(t) in Eq. (6) with Eq. (34), which show reasonable solutions. 7. Conclusion In the study reported here, an inverse problem is formulated for the simultaneous identification of the functional form of nonlinear damping as well as the phase shift and amplitude of the external harmonic forcing from measured data of a forced oscillatory motion. However, mathematically, it is involved in solving the firstkind integral equation. Thus, the identifying process is ill-posed in the sense of stability, which leads to a numerical instability that will influence the performance of the simultaneous identification. This implies that some small amount of noisy response data can be arbitrarily amplified, resulting in unreliable solutions. The difficulty is remedied by the use of Landweber’s regularization method, where the appropriate number of iterations for the regularization is determined by applying the L-curve criterion in the presence of noisy data. As a model equation, a forced van der Pol equation is used to determine the workability of the identification. Even though, the example provided herein is a highly nonlinear system, the method proposed in this paper for the simultaneous identification yields a good approximation of the exact expression in mea-
85
suring the full nonlinear damping and the phase shift and amplitude of the harmonic excitation in a stable but accurate manner. That is, we have shown that through the numerical experiment the proposed identification procedure works well, however, in this study, a real experiment was not carried out with the purpose of proving the workability, because it may be at times extremely difficult to do. For future work, it is necessary to test the workability of the proposed identification with a real experiment or an alternative from a different source. It is also required to extend the identification, limited to a forced nonlinear single degree of freedom system, to more complex systems, e.g., multiple degree of freedom systems. Acknowledgement The present paper is supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Grant no.: 2011-0010090). In addition, this research was supported by the Lloyd’s Register Educational Trust Research Centre of Excellence at Pusan National University and also by the National Research Foundation (NRF) grant funded by the Korea government (MEST) (Grant no.: K20901000005-09E0100-00510). Lloyd’s Register Educational Trust is an independent charity working to achieve advances in transportation, science, engineering and technology education, and training and research worldwide for the benefit of all. Finally, the author would like to extend special thanks to the student, Mr. Jinsoo Park in the Department of Naval Architecture and Ocean Engineering, Pusan National University for his help for the manuscript.
References [1] Iourtchenko DV, Duval L, Dimentberg MF, The damping identification for certain SDOF systems, developments in theoretical and applied mechanics, In: Proceedings of SECT AM XX conference, Auburn university, 2000. [2] Iourtchenko DV, Dimentberg MF. In-service identification of non-liner damping from measured random vibration. J Sound Vib 2002;255:549–54. [3] Kang JS, Park SK, Shin SB, Lee HS. Structural system identification in time domain using measured acceleration. J Sound Vib 2005;288:215–34. [4] Lee XY, Law SS. Identification of structural damping in time domain. J Sound Vib 2009;328:71–84. [5] Mohammad KS, Worden K, Tomlinson GR. Direct parameter estimation of linear and non-linear structures. J Sound Vib 1992;152:471–99. [6] Tomme JV, Evaluation of damping measurements in materials and structures, In: 13th International seminar on modal analysis, 1988. [7] Lazan BJ. Damping of materials and members in structural mechanics. Oxford Pergamon Press; 1968. [8] Naprstek J, Identification of linear and non-linear dynamic systems with parametric noises, in: Proceedings of the 4th international conference on stochastic structural dynamics, 1999. [9] Jang TS, Choi HangS, Han SL. A new method for detecting nonlinear damping and restoring forces in nonlinear oscillation systems from transient data. Int J Non-Linear Mech 2009;44:801–8. [10] Jang TS. Non-parametric simultaneous identification of both the nonlinear damping and restoring characteristics of nonlinear systems whose dampings depend on velocity alone. Mech Syst Signal Process 2011;25:1159–73. [11] Jang TS, Baek HS, Choi HS, Lee SG. A new method for measuring nonharmonic periodic excitation forces in nonlinear damped systems. Mech Syst Signal Process 2011;25–6:2219–28. [12] Ugol’kov VN. Methods of measuring the phase shift and amplitude of harmonic signals using integral samples. Meas Tech 2003;46:495–501. [13] Stakgold I. Boundary value problems of mathematical physics, vol. 1. SIAM; 1967. [14] Jang TS, Han SL. Application of Tikhonov’s regularization to an unstable two dimensional water waves: spectrum with compact support. Ships Offshore Struct 2008;3:41–7. [15] Groetsch CW. Inverse problems in the mathematical sciences. Vieweg; 1993. [16] Hansen PC. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev 1992;34:561–80. [17] Kirsch A. An introduction to the mathematical theory of inverse problems. Berlin: Springer; 1996.
Contents lists available at SciVerse ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
A method for simultaneous identification of the full nonlinear damping and the phase shift and amplitude of the external harmonic excitation in a forced nonlinear oscillator T.S. Jang ⇑ Naval Architecture and Ocean Engineering, Pusan National University, Busan 609-735, Republic of Korea
a r t i c l e
i n f o
Article history: Received 4 October 2012 Accepted 14 February 2013 Available online 15 March 2013 Keywords: Simultaneous identification Non-parametric identification Full nonlinear damping Phase shift and amplitude Lack of stability Regularization
a b s t r a c t We propose a new general method that is able to clarify the simultaneous identification of the full nonlinear damping characteristics as well as the phase shift and amplitude of the external harmonic excitation in a forced nonlinear single degree system. The method is the non-parametric identification based on the concept of zero-crossings regarding the system responses. It involves an integral equation of the first kind, resulting in the lack of stability properties. The difficulty of instability is overcome by the stabilization technique. The workability of the proposed method is demonstrated through a numerical experiment consisting of a highly nonlinear system. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Associated with the mathematical modeling of dynamical systems, a considerable amount of attention has been given to the system identification of nonlinear physical systems in various branches of science and engineering. In fact, an appropriate and correct mathematical modeling from measured response data through identification plays an essential role in achieving accurate predictions of motion responses to various loading environments, which is also important for improving the system performance such as the effective control and the estimation of the structure lifetime. Especially, the damping identification is extremely crucial to the forced systems in the field of structural engineering, mainly because it can be directly applied to the early-stage design of nonlinear physical systems. A number of methods have been proposed in the literature sourced for the damping identification of nonlinear systems arising in mechanics, whose primary objective involves the formulation of a mathematical model on the basis of measured data. Iourtchenko et al. [1] successfully presented a novel technique for an identification scheme in the case of parametric excitation. Iourtchenko et al. [2] further proposed a mathematical procedure based on the method of stochastic
⇑ Tel.: +82 51 510 2789; fax: +82 51 581 3718. E-mail address: [email protected] 0045-7949/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruc.2013.02.008
averaging for in-service identification of the damping characteristic from measured random vibration. Recently, Kang et al. [3] estimated the stiffness and damping parameters of a structure by using the measured acceleration. Lee and Law [4] proposed a procedure for the viscous damping model adopting the inverse modal transformation. With regard to the modeling of nonlinear damping systems, there were also the challenging issues on other studies [5–8]. What this paper involves is the simultaneous identification of nonlinear physical systems, i.e., a study to identify more than a ‘‘single’’ physical element of a system. This is in contrast to the usual current identification. In fact, most related studies on the system identification, however, are primarily limited to the identification of a ‘‘single’’ physical element of a system, say, either a system’s damping force or restoring force. Even though only few attempts have so far been made at the simultaneous identification, a recent work may be found in the paper of Jang et al. [9], where they newly proposed a general procedure for identifying both a full nonlinear restoring force and a full nonlinear damping function in _ in which D(y) is a general damping function of the form, DðyÞ y, the displacement y and y_ is the velocity in nonlinear oscillation systems. In addition, Jang [10] modified the general procedure _ as well [9] to identify a functional form of nonlinear damping, DðyÞ, as a full nonlinear restoring force. In this paper, we propose a new mathematical procedure to simultaneously measure not only the full nonlinear damping in _ but the phase shift and amplitude of external the form, DðyÞ y; harmonic excitation in a forced nonlinear single degree system.
78
T.S. Jang / Computers and Structures 120 (2013) 77–85
The procedure proposed here is thus regarded as a simultaneous identification of two physical elements of a system. However, it is different from the methods of the simultaneous identification [9,10] in that the two physical elements to be identified in the present study are the damping and the external force of a system. Of course, the procedure is also distinct from the (usual) force identification of a single physical element of the external force of a system. A method for force identification is, for example, illustrated in Jang et al. [11], who proposed a new method to find nonharmonic periodic excitation forces in nonlinear damped systems. When we often encounter vibrations due to the harmonic excitation, it may be hard to perceive the phase shift and amplitude of external harmonic loading beforehand [12]. Hence, in practice, the present work may be the one that has a wider applicable range than the existing methods identifying only the single damping element of a system. Furthermore, the procedure proposed is a fully nonlinear method. Consequently, it does not depend on a small parameter and thus overcome the disadvantage and limitations of the parametric methods such as a perturbation expansion procedure. Finally, it is found that a numerical experiment demonstrates the workability of the proposed method of the measuring procedure in nonlinear oscillations by using a forced van der Pol equation with cubic nonlinear restoring force. 2. The equation of motion We consider a nonlinear damped system, subject to the sum of two external harmonic loads, asin(xt) and bsin(xt), where the positive constants a, b and x are amplitudes and frequency. System dynamics considered are governed by Newton’s equation of motion, expressed as in the following form
€ þ DðyÞy_ þ RðyÞ ¼ a sinðxtÞ þ b cosðxtÞ for t > 0: my
ð1Þ
Mathematically, Eq. (1) is a second-order nonlinear ordinary differential equation, where the constant m > 0 represents the mass of a particle of an oscillator, the term involving D(y) the general form of nonlinear damping and the term R(q) a general nonlinear restoring force. The harmonic excitation in Eq. (1) can be, alternatively, represented in terms of the phase shift a and the amplitude r as follows
a sinðxtÞ þ b cosðxtÞ ¼ r sinðxt þ aÞ
ð2Þ
where
r¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a2 þ b
3.1. Pseudo-spring constant We attempt to develop a normal differential form for Eq. (1), which plays an essential role in deriving an integral equation in this study. To this end, we first introduce a pseudo-spring constant, kp > 0, and add the linear spring force of the form kpy to both sides of Eq. (1) to have
€ þ kp y ¼ fDðyÞy_ RðyÞ þ f ðtÞg þ kp y; my
t>0
ð5Þ
where
f ðtÞ a sinðxtÞ þ b cosðxtÞ:
ð6Þ
Eq. (5) is then simply rewritten as
€ þ kp y ¼ FðtÞ my
ð7Þ
with
FðtÞ DðyÞy_ RðyÞ þ f ðtÞ þ kp y:
ð8Þ
Eq. (7) is viewed as a normal (differential) form associated with Eq. (1) in this paper, which corresponds to an equivalent integral equation, as will be discussed in Section 3.2. 3.2. Volterra integral transform For simplicity, we assume that the motion is initially at rest for time t < 0 such that
_ yð0Þ ¼ 0:
yð0Þ ¼ 0;
ð9Þ
The displacement y(t) in Eq. (7) is then expressed as the Volterra integral transform of F in Eq. (8) with the kernel G(t, s:kp) [10]
yðtÞ ¼
Z
t
Gðt; s : kp ÞFðsÞds;
ð10Þ
0
and
Gðt; s : kp Þ ¼
1 y1 ðsÞy2 ðtÞ y1 ðtÞy2 ðsÞ : m y1 ðsÞy_ 2 ðsÞ y_ 1 ðsÞy2 ðsÞ
ð11Þ
ð3Þ
In Eq. (11), y1(t) and y2(t) are the motion responses of the linear oscillator with the pseudo-spring constant kp:
ð4Þ
€1 ðtÞ þ kp y1 ðtÞ ¼ 0; my €2 ðtÞ þ kp y2 ðtÞ ¼ 0; my
and
tan a ¼ b=a:
study: Given a motion response of a forced nonlinear system, is it possible to identify the nonlinear damping of the system and the phase shift and amplitude of harmonic excitation at the same time? To answer the question, we begin with an integral equation.
As mentioned briefly in introduction, the correct and accurate measurement of nonlinear damping of a system plays an important role in applications in a variety of physical problems in different disciplines. However, as far as a forced system is concerned, it is practically difficult to know the phase shift a in Eq. (4) and amplitude r in Eq. (3) of external loading directly from experiment in advance. In fact, it is not easy to install transducers within a moving structure during its construction, so that there are only few cases in which we can directly determine external loading in forced oscillations [11]. To solve the difficulty, a new general method will be proposed for measuring both the full nonlinear damping and the phase shift a and amplitude r of external harmonic load during forced nonlinear oscillations. 3. Integral equation The present study may be considered as an inverse problem in the form of a question, which is the basic idea underlying this
y1 ð0Þ ¼ 1; y_ 1 ð0Þ ¼ 0 y2 ð0Þ ¼ 0; y_ 2 ð0Þ ¼ 1:
ð12Þ
In general, the identification of nonlinear systems relies on measuring system’s responses. In this work, we assume that the displace_ ment y(t) and velocity yðtÞ of the system are measured during oscillations. Mathematically, this implies that the transform of Eq. (10) may also be regarded as an integral equation for F(t) in Eq. (8), because the left hand side of Eq. (10) is known from the measurement. Regarding the integral equation, it is meaningful to check whether Eq. (10) has a unique solution F(t), which was discussed and proved in Jang et al. [9]. It is noted that if Eq. (8) is rewritten as
f ðtÞ DðyÞy_ ¼ FðtÞ þ RðyÞ kp y;
ð13Þ
the right hand side of Eq. (13) becomes then known, since the (unique) F can be determined by solving the integral equation of Eq. (10). Denoting the known by u(t),
uðtÞ FðtÞ þ RðyÞ kp y; we have an alternative expression for (13), namely
ð14Þ
T.S. Jang / Computers and Structures 120 (2013) 77–85
_ uðtÞ ¼ f ðtÞ DðyÞy:
ð15Þ
4. Simultaneously measuring procedure
79
Then, the value of the nonlinear damping function D, evaluated at y = y(ti), can be explicitly calculated as
Dðyðt i ÞÞ ¼
f ðt i Þ uðti Þ _ iÞ yðt
ð24Þ
In this section, we present a mathematical procedure for the simultaneous measuring of both the nonlinear damping and the phase shift and amplitude of external excitation, the main idea of which is provided with aid of Eq. (15).
according to Eq. (15). To identify D(y) in more points of y, the above process can be repeated as many times as necessary for each ti, resulting in the desired measurement of the unknown D.
4.1. Phase shift and amplitude of harmonic excitation
5. Lack of stability
Notice that the external force can be balanced by the (known) u(t), whenever the damping term in Eq. (15) disappears. Actually, _ this happens when the velocity yðtÞ of the system response vanishes during oscillations. For a systematic analysis, let us define a _ zero-crossing time t = t1, t2,. . . in terms of the velocity yðtÞ of the system response, that is, it satisfies
Even though we have constructed a measuring procedure for the phase shift and amplitude of harmonic excitation and the functional form of nonlinear damping in the previous section, it is still necessary to discover u(t) in Eq. (15), equivalently, F(t) in Eq. (10). However, Eq. (10) is classified as an integral equation of the ‘‘first’’ kind, which leaves the question of stability or instability open.
_ i Þ ¼ 0 for 0 < t1 < t2 < . . . : yðt
5.1. Hilbert–Schmidt integral operator K
ð16Þ
Then, the damping force becomes zero
_ iÞ ¼ 0 Dðyðt i ÞÞ yðt
ð17Þ
when t = t1, t2,. . .. This immediately yields
uðti Þ ¼ f ðt i Þ;
t ¼ t 1 ; t2;...
y ¼ KðFÞ; ð18Þ
from Eq. (15).If we choose two first zero-crossing times in terms _ of the velocity yðtÞ, say, t = t1, t2, then Eq. (18) turns out to be a linear matrix equation for the unknown a and b appearing in Eq. (6):
uðt 1 Þ uðt 2 Þ
¼M
a
ð19Þ
b
where M denotes the two by two matrix
M¼
ð25Þ
where K is defined as the integral operator
KðFÞ
Z
t
Gðt; s : kp ÞFðsÞds:
ð20Þ
Assuming that the matrix M in Eq. (20) is nonsingular, the two amplitudes of a and b can be calculated by using the inverse M1:
ð26Þ
0
In what follows, we shall investigate the topological structure of the operator K in Eq. (26). We start by noticing that the responses of y1(t) and y2(t) are harmonic functions and the denominator of the kernel G in Eq. (11) is a Wronskian which is a non-zero constant [9]:
W ¼ y1 ðtÞy_ 2 ðtÞ y_ 1 ðtÞy2 ðtÞ ¼ const:
sinðxt1 Þ; cosðxt1 Þ : sinðxt2 Þ; cosðxt2 Þ
For convenience, Eq. (10) can be rewritten as the operator equation in a compact way
ð27Þ
This immediately leads to the inequality with regard to the kernel G
Z 0
T
Z
T
G2 ðt; s : kp Þdsdt < 1
ð28Þ
0
The amplitudes, a and b, so calculated, in turn, enable the phase shift a and the amplitude r in Eqs. (3) and (4) to be completely determined.
for an arbitrary (large) constant T > 0. We recall that a kernel satisfying the inequality in Eq. (28) is called a Hilbert–Schmidt kernel and the integral operator with a Hilbert–Schmidt kernel is known as a Hilbert–Schmidt integral operator [13]. Thus, it is found that the kernel G in Eq. (11) and the operator K in Eq. (26) are a Hilbert–Schmidt kernel and a Hilbert–Schmidt integral operator, respectively.
4.2. Nonlinear damping
5.2. Unbounded operator K1: discontinuity
Since the forcing f(t) in Eq. (6) is shown to be determined through Eq. (21) in Section 4.1, it is possible to estimate the following ratio
According to functional analysis, a Hilbert–Schmidt integral operator is shown to be a compact operator, that is, it transforms bounded sets into compact sets. Furthermore, a compact operator, defined on infinite dimensional spaces such as function spaces, is continuous, while its inverse is not continuous. Therefore, the operator K in Eq. (26) is compact and thus K1 is a discontinuous operator (or unbounded operator). Accordingly, F(t) in Eq. (10), represented in terms of the displacement y(t) as the below equation, does not depend continuously on the y(t) in Eq. (7), lacking stability properties of solutions [14,15]:
a b
¼ M 1
uðt 1 Þ
uðt 2 Þ
f ðtÞ uðtÞ ; _ yðtÞ
ð21Þ
ð22Þ
_ provided that yðtÞ – 0. In addition, we note that, from Eq. (15), the ratio is identical to the damping function D(y) that we are looking for. Thus, the estimation of the ratio implies the measurement of D(y), which is accomplished by using a proper time with regard to system responses. In contrast to a zero-crossing time t = t1, t2,. . . as in Section 4.1, we shall select a non-zero-crossing time in terms of the velocity y_ of the system response during oscillations, denoted as t = t1, t2, . . ., such that
_ i Þ – 0 for 0 < t 1 < t 2 < . . . : yðt
ð23Þ
F ¼ K 1 y
ð29Þ
Due to the discontinuity of K1, a small perturbation in y(t) or a small noisy error can be extremely amplified so that unreliable solutions may result from a small amount of noisy data. In fact, although Eq. (10) may be directly discretized into a numerical system of simultaneous linear equations (for example, by using the
80
T.S. Jang / Computers and Structures 120 (2013) 77–85
Fig. 1. Two system responses of the dynamic oscillatory system with the initial conditions of Eq. (35): c = 4.
Fig. 3. A typical graphical illustration of the L-curve: c = 4 and d = 0.01. Fig. 2. Phase diagram corresponding to the dynamic responses in Fig. 1: c = 4.
Table 1 Recovered a and b (or, r and a) of excitations for different noise level d: c = 4. Recovered
Exact d = 0.01 d = 0.05 d = 0.10 d = 0.15
a
b
4 3.7667 4.4980 3.5749 4.5311
8 8.0082 7.8518 7.5785 8.8896
r
a
Average error (%)
8.9443 8.8442 9.0489 8.3791 9.9778
1.1071 1.1308 1.0506 1.1300 1.0994
– 2.9 6.8 7.9 12.2
trapezoidal numerical integration rule), the directly discretized numerical system turns out to be ill-conditioned [14].
we thus would like an approximation for the solution to depend continuously on y(t). Tikhonov’s regularization is the most commonly used method of regularization [15]. In this study, however, we suggest the use of the Landweber’s regularization, which is known as a simple but effective iterative method. Thereby, the solution F in Eq. (10) is obtained as a limit of the sequence {Fj(t)} in the following iterative form of operator equation
F j ðtÞ ¼ F j1 ðtÞ kL fLðF j1 ðtÞÞg þ kL ðyðtÞÞ;
In order to prevent the lack of stability from affecting the performance of the present study of simultaneous measurement, we will suppress the instability through a stabilization (or regularization) technique. The main idea of the method of regularization is to replace an (unstable) integral equation of the ‘‘first’’ kind by a nearby (stable) integral equation of the ‘‘second’’ kind. Mathematically,
ð30Þ
for a real positive constant, k > 0. In Eq. (30), the operator L is the Volterra integral transform in the form in Eq. (10)
LðzÞ 5.3. Landweber’s regularization: a stabilization technique
j ¼ 1; 2 . . .
Z
t
Gðt; s : kp ÞzðsÞds
ð31Þ
0
for an arbitrary function z(t). The symbol L⁄ stands for the adjoint operator of L, and I the identity operator [15]. Based on the fact that Landweber’s regularization converges to the solution for an arbitrary initial guess F0(t), we start in this study with the initial guess of zero function, i.e., F0(t) 0, without loss of generality [15]. According to the theory of regularizations, there exist an optimal number of iterations in Eq. (30) for finding an accurate
T.S. Jang / Computers and Structures 120 (2013) 77–85
81
Fig. 4. Convergence behaviors of Fj(t): c = 4 and d = 0.01 (a) j = 104, (b) j = 105, (c) j = 106 and (d) j = 108.
and stable solution in the presence of noisy data. We shall apply Lcurve criterion [16] to determine the optimal number of iterations in Eq. (30), which will be discussed in more detail in Section 6. 6. Numerical experiment In this section, the workability of the proposed method will be demonstrated for measuring simultaneously the full nonlinear damping and the phase shift and amplitude of the external harmonic excitation during forced nonlinear oscillations. To this end, a model equation is chosen for the numerical experiment. 6.1. Forced Van der Pol equation
Fig. 5. Recovered f(t) asin(xt) + bcos(xt) compared with the exact one: c = 4 and d = 0.01.
We introduce a forced Van der Pol equation having a cubic nonlinear restoring to the numerical experiment, which evolves in time according to the second order differential equation.
€ þ cðy2 1Þy_ þ k1 y þ k2 y3 ¼ a sinðxtÞ þ b cosðxtÞ: my
ð32Þ
In the case of k2 = a = b = 0, Eq. (32) reduce to the Van der Pol equation which has a long history of being used in physical sciences, as non-conservative systems with nonlinear damping. The immediate comparison of equations (1) and (32) reveals the functional forms of D(y) and R(y) in Eq. (1), represented as
DðyÞ ¼ cðy2 1Þ and RðyÞ ¼ k1 y þ k2 y3 :
ð33Þ
For convenience and without loss of generality, the mass m in Eq. (32) is normalized as unit. The values for the spring constants k1 and k2 are taken as 50 and 5, respectively, and the amplitudes of excitation are set to be
a ¼ 4 and b ¼ 8:
Fig. 6. Recovered D(y): c = 4 and d = 0.01.
ð34Þ
We further select the interval of time domain as 0 < t < 1.2. Figs. 1 and 2 depict the time evolution solutions for Eq. (32) and their phase diagrams, respectively: Runge–Kutta type integration
82
T.S. Jang / Computers and Structures 120 (2013) 77–85
Fig. 7. Recovered f(t) asin(xt) + bcos(xt) compared with the exact ones for four different various frequencies: c = 4 and d = 0.01.
Fig. 8. Recovered D(y) compared with the exact one for four different various frequencies: c = 4 and d = 0.01.
schemes are used as numerical integration methods, in which the value for the coefficient c is taken as 4 and the initial condition is imposed by
yð0Þ ¼ 0;
_ yð0Þ ¼ 0:
ð35Þ
With data depicted in Figs. 1 and 2, we are now ready to follow the (inverse) measuring procedure proposed in the previous sections. By using the data of Figs. 1 and 2, we will attempt to identify the damping and forcing of a forced Van der Pol equation in Eq. (32) and the results are to be compared with the corresponding exact expressions, D(y) = c(y2 1), in Eq. (33) and f(t) in Eq. (6) with the coefficients of Eq. (34), respectively.
Fig. 9. Recovered f(t) asin(xt) + bcos(xt) and D(y) for 3 cases of pseudo-spring constant kp = 1, 3 and 5.
T.S. Jang / Computers and Structures 120 (2013) 77–85
83
Fig. 10. Recovered f(t) asin(xt) + bcos(xt) and D(y) for four different noise levels when jopt = 106: c = 4.
6.2. Noisy data and L-curve criterion As mentioned earlier, for the present simultaneous measurement, we are required to measure the displacement and velocity data of motion responses beforehand. In practice, however, the measured data are always deteriorated to some extent by noisy error. This implies that the left-hand side in the integral Eq. (10) is never exactly known but only up to an error of, say, noise level
d > 0. We thus assume that we know, from measuring devices, d > 0 and noisy data yd(t) such that
kyðtÞ yd ðtÞk2 6 d;
ð36Þ
in which the symbol kk2 denotes L2 norm [9–11,17]. In this study, four cases of noise levels of noisy data are randomly generated, which are tabulated in Table 1. Generally, in
84
T.S. Jang / Computers and Structures 120 (2013) 77–85
Table 2 Recovered a and b (or, r and a) of excitations for different damping coefficients c when d = 0.01. c
Case 1 Case 2 Case 3
2 3 4
Recovered a
b
4.426 4.037 3.767
7.623 8.292 8.002
r
a
Average error (%)
8.814 9.222 8.442
1.045 1.118 1.130
7.7 2.3 2.9
the presence of noisy data, Landweber’s regularization approaches the true solution in the initial stage of iteration, but it potentially deviate far away at a large number of iterations. As a result, the iterated solution of Landweber’s regularization goes to a useless solution because of the amplification of hidden noise. Depending on a given noise level d > 0, the accuracy of the Landweber’s regularization can be influenced by the number of iterations. One of the most effective ways to find the optimal number of iterations may be of an L-curve criterion [16], which mainly concerns how to select an appropriate iterated solution. A brief sketch is as follows. If we consider the two quantities
rj ¼ LðF j ðtÞÞ yd ðtÞ2
and sj ¼ F j ðtÞ2
j ¼ 1; 2 . . .
ð37Þ
which correspond to the norms of the residual and the iterated solution, respectively, then, the residual rj tends to decrease, while the value of sj increases as the number of iterations increase. This can be best illustrated by introducing a plot of the pair of (rj, sj), which may yield a curve exhibiting a typical ‘‘L’’ shape. According to L-curve criterion [16], an optimal number of iterations are considered to be the one corresponding to the corner of an L-shaped curve. Fig. 3 shows a typical shape of the L-curve in the present study, where the convergence constant k > 0 is chosen as k ¼ 1: here, c = 4 and d = 0.01. In relating the calculation of Eq. (31), essential for the Landweber’s iteration of Eq. (30), the pseudo-spring constant is set to be unit in this section, i.e., kp = 1. The optimal iteration number jopt appears to be clear at the corner of the L-curves as shown in Fig. 3, which is jopt = 106. The convergence behaviors of iterative solutions Fj(t) in the iteration Eq. (30) are illustrated in Fig. 4. Here, the Fj(t) are compared with the exact one, F(t) in Eq. (8), corresponding to Eq. (33). We see that not-enough iteration
Fig. 11. Recovered f(t) asin(xt) + bcos(xt) and D(y) for three different coefficients of c when jopt = 106.
T.S. Jang / Computers and Structures 120 (2013) 77–85
produces a poor but stable approximation, while too many iterations result in a worse (unstable) solution. 6.3. Measuring D(y), a and r Following the simultaneous measuring procedure as proposed in Section 4, we try to calculate the first two zero-crossing times t1, t2 in terms of the velocity y_ from data in Fig. 1. The results are as follows
t 1 ¼ 0:46 s; t 2 ¼ 0:91 s:
ð38Þ
Using the two times t1, t2, we can determine the amplitudes of the two unknowns a and b in Eq. (6) by solving the matrix equation of Eq. (21) and thus find out the phase shift a and the amplitude r from Eqs. (3) and (4). The recovered ones are a = 3.7667, b = 8.0082, equivalently, a = 1.1308rad and r = 8.8442, which have the average error about 3%. Fig. 5 depicts the recovered external harmonic excitation f(t) asin(xt) + bcos(xt) in Eq. (6) compared to the exact one, which is f(t) in Eq. (6) by using Eq. (34). Next, we select 25 non-zero-crossing times of Eq. (23) in terms _ of the velocity yðtÞ of the system response during oscillation. With the selected non-zero-crossing times, we apply Eq. (24) to recover the nonlinear damping function D(y). The results are shown in Fig. 6, where it is relatively accurate when compared with the exact results D(y) = c(y2 1) of Eq. (33). Figs. 7 and 8 also depict the identified f(t) in Eq. (6) and D(y) for the optimal iteration number jopt = 106, respectively. However, four different various frequencies are used for the results, which are in a good agreement with the exact ones of f(t) in Eq. (6) with Eq. (34) and D(y) = c(y2 1) of Eq. (33), respectively. We notice that it is clear from Eq. (5) that D(y) and f(t) are invariant with respect to the pseudo-spring constant kp, that is, they do not depend on a choice of kp. This is also numerically confirmed as shown in Fig. 9, where D(y) and f(t) are calculated based on the choices of the three different kp = 1, 3 and 5. It may be hard to see any effect of a pseudo-spring constant kp on D(y) and f(t), as expected. In Fig. 10, we examine f(t) and D(y) according to the four different noise levels d as listed in Table 1. Finally, we pick up three cases of different coefficients for c in Table 2, to check their dependency on f(t) and D(y). As exhibited Fig. 11, the results are compared with exact values of f(t) in Eq. (6) with Eq. (34), which show reasonable solutions. 7. Conclusion In the study reported here, an inverse problem is formulated for the simultaneous identification of the functional form of nonlinear damping as well as the phase shift and amplitude of the external harmonic forcing from measured data of a forced oscillatory motion. However, mathematically, it is involved in solving the firstkind integral equation. Thus, the identifying process is ill-posed in the sense of stability, which leads to a numerical instability that will influence the performance of the simultaneous identification. This implies that some small amount of noisy response data can be arbitrarily amplified, resulting in unreliable solutions. The difficulty is remedied by the use of Landweber’s regularization method, where the appropriate number of iterations for the regularization is determined by applying the L-curve criterion in the presence of noisy data. As a model equation, a forced van der Pol equation is used to determine the workability of the identification. Even though, the example provided herein is a highly nonlinear system, the method proposed in this paper for the simultaneous identification yields a good approximation of the exact expression in mea-
85
suring the full nonlinear damping and the phase shift and amplitude of the harmonic excitation in a stable but accurate manner. That is, we have shown that through the numerical experiment the proposed identification procedure works well, however, in this study, a real experiment was not carried out with the purpose of proving the workability, because it may be at times extremely difficult to do. For future work, it is necessary to test the workability of the proposed identification with a real experiment or an alternative from a different source. It is also required to extend the identification, limited to a forced nonlinear single degree of freedom system, to more complex systems, e.g., multiple degree of freedom systems. Acknowledgement The present paper is supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Grant no.: 2011-0010090). In addition, this research was supported by the Lloyd’s Register Educational Trust Research Centre of Excellence at Pusan National University and also by the National Research Foundation (NRF) grant funded by the Korea government (MEST) (Grant no.: K20901000005-09E0100-00510). Lloyd’s Register Educational Trust is an independent charity working to achieve advances in transportation, science, engineering and technology education, and training and research worldwide for the benefit of all. Finally, the author would like to extend special thanks to the student, Mr. Jinsoo Park in the Department of Naval Architecture and Ocean Engineering, Pusan National University for his help for the manuscript.
References [1] Iourtchenko DV, Duval L, Dimentberg MF, The damping identification for certain SDOF systems, developments in theoretical and applied mechanics, In: Proceedings of SECT AM XX conference, Auburn university, 2000. [2] Iourtchenko DV, Dimentberg MF. In-service identification of non-liner damping from measured random vibration. J Sound Vib 2002;255:549–54. [3] Kang JS, Park SK, Shin SB, Lee HS. Structural system identification in time domain using measured acceleration. J Sound Vib 2005;288:215–34. [4] Lee XY, Law SS. Identification of structural damping in time domain. J Sound Vib 2009;328:71–84. [5] Mohammad KS, Worden K, Tomlinson GR. Direct parameter estimation of linear and non-linear structures. J Sound Vib 1992;152:471–99. [6] Tomme JV, Evaluation of damping measurements in materials and structures, In: 13th International seminar on modal analysis, 1988. [7] Lazan BJ. Damping of materials and members in structural mechanics. Oxford Pergamon Press; 1968. [8] Naprstek J, Identification of linear and non-linear dynamic systems with parametric noises, in: Proceedings of the 4th international conference on stochastic structural dynamics, 1999. [9] Jang TS, Choi HangS, Han SL. A new method for detecting nonlinear damping and restoring forces in nonlinear oscillation systems from transient data. Int J Non-Linear Mech 2009;44:801–8. [10] Jang TS. Non-parametric simultaneous identification of both the nonlinear damping and restoring characteristics of nonlinear systems whose dampings depend on velocity alone. Mech Syst Signal Process 2011;25:1159–73. [11] Jang TS, Baek HS, Choi HS, Lee SG. A new method for measuring nonharmonic periodic excitation forces in nonlinear damped systems. Mech Syst Signal Process 2011;25–6:2219–28. [12] Ugol’kov VN. Methods of measuring the phase shift and amplitude of harmonic signals using integral samples. Meas Tech 2003;46:495–501. [13] Stakgold I. Boundary value problems of mathematical physics, vol. 1. SIAM; 1967. [14] Jang TS, Han SL. Application of Tikhonov’s regularization to an unstable two dimensional water waves: spectrum with compact support. Ships Offshore Struct 2008;3:41–7. [15] Groetsch CW. Inverse problems in the mathematical sciences. Vieweg; 1993. [16] Hansen PC. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev 1992;34:561–80. [17] Kirsch A. An introduction to the mathematical theory of inverse problems. Berlin: Springer; 1996.