main.dvi

Transcription

1 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 Force State Maps Using Reproducing Kernel Particle Method and Kriging Based Functional Representations Vikas amdeo 1, and C S Manohar 1,3 Abstract: The problem of identification of nonlinear system parameters from measured time histories of response under known excitations is considered. Solutions to this problem are obtained by using the force state mapping technique with two alternative functional representation schemes. These schemes are based on the application of reproducing kernel particle method (RKPM) and kriging techniques to fit the force state map. The RKPM has the capability to reproduce exactly polynomials of specified order at any point in a given domain. The kriging based methods represent the function under study as a random field and the parameters describing this field are optimally determined based on available observations. The present study investigates the performance of RKPM and kriging based fits to the force state maps for a variety of nonlinear dynamical systems. The study also examines the application of force state maps in determining the fixed points limit cycles of the system and their stability, determining the properties of the linear system which would result if nonlinearity were to be absent, and (c) dealing with nonlinearities that are continuous but not differentiable and nonlinearities that are discontinuous at a set of points within the domain of interest. Illustrative examples on single and multi-degrees of freedom nonlinear systems are presented to demonstrate the scope of the proposed procedures. Keyword: Force state map. System identification. RKPM. Kriging. 1 Department of Civil Engineering, Indian Institute of Science, Bangalore India Research Student (vikasnamdeo@yahoo.co.in). 3 Corresponding author. Professor (manohar@civil.iisc.ernet.in) 1 Introduction Problems of nonlinear dynamical system identification are presently one of the most challenging problems in structural mechanics. Currently there exists wide ranging methods that work in time or frequency domains to tackle this problem (Kerschen et al., 006). Force state mapping is one of the time domain methods for nonlinear system identification that can be viewed as a nonparametric method for system identification and this method constitutes the focus of the present study. The essence of the method can be explained by considering the dynamics of a single degree of freedom (SDOF) nonlinear system mẍ+g(x,ẋ)=g(t); x(0)=x 0, ẋ(0)=ẋ 0 (1) Assuming that the mass of the system is known and that we are able to measure the applied force g(t) and the response acceleration ẍ(t), an expression for the unknown restoring force G(x,ẋ) can be obtained as G[x(t),ẋ(t)] = g(t) mẍ(t) () Furthermore, if x(t) and ẋ(t) are also measured, the function G[x(t),ẋ(t)] can be plotted as a surface over the phase plane plot of x(t) versus ẋ(t) and this surface is called the restoring force surface. It is of interest to note that the nature of restoring force surface G(x, ẋ) remains unchanged if excitation g(t) or the initial conditions x(0) and ẋ(0) are changed. Figure 1 shows a typical restoring force surface for a nonlinear system along with the associated phase plane plot. The idea of force state map method is to fit the above surface using suitable set of known functions which, in turn, lead to the characterization of the system. Figure 1 also shows the restoring force surface that has been fitted numerically based on

2 14 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 Figure 1: Illustration of a typical force state map; phase plane plot; the associated force state map. observed x(t), ẋ(t) and computed G[x(t),ẋ(t)]. The idea of force state map was originally proposed by Masri and Caughey (1979), Crawley and Aubert (1986), and Crawley and O Donnell (1986,1987a,b). Masri and Caughey (1979) proposed that the restoring force surface be fitted using Chebychev polynomials as G[x(t),ẋ(t)] = m n i=0 j=0 C ij T i (x)t j (ẋ) (3) Crawley and Aubert (1986) advocated the use of measured force state maps, either as plain look-up tables, or, as a mathematically fitted surface, to characterize nonlinear behavior of joints. These authors also noted the versatility of the approach to handle arbitrary strong nonlinearity. Crawley and O Donnell (1987) considered the force state mapping method as a three step procedure involving selection and application of the excitation force, processing of the system response to construct the force state map, and (c) extraction of the pertinent system parameters from the map. These authors note the need to measure simultaneously the quantities x(t), ẋ(t), ẍ(t) and g(t) and have discussed the possible use of state estimation methods to consistently estimate all the states when measurements are fewer than what is needed. Furthermore, assuming that the form of the nonlinearity is known, the authors used the least squares algorithm to fit the restoring force surface. Masri et al., (198, 1987a, b) considered the problem of identification of parameters of multidegree of freedom (MDOF) systems using force state mapping techniques. Here the normal modes of the linear system are used to expand the response of the nonlinear system. The resulting equations in the generalized coordinates are clearly expected to be mutually coupled. The authors, however, make an ad hoc simplifying assumption that the coupling terms in the generalized equations can be ignored, which, in turn, enables the construction of a set of force state maps for the resulting uncoupled dofs. The difficulties associated with the simultaneous measurement of displacement, velocity and acceleration, as is required in the implementation of the force state mapping method, has been discussed by Worden (1990a, b). The author discusses numerical strategies for estimating two of the required quantities in terms of measurement on the remaining one and also examines the issues related to the choice of the test signals in achieving this. In a series of papers, Al-Hadid and Wright (1989,1990,199) have examined the following issues: the performance of the method when ordinary polynomials are used to for the force state map instead of the orthogonal polynomials such as the Chebychev polynomials, identification of location of discrete nonlinearities, (c) examination of the performance of the method when applied to the study of a T-beam, and (d) the performance of the method vis-à-vis

3 Force State Maps 15 the need for the knowledge of the mass or the modal mass. The possibility of developing the force state map method in frequency domain for identification of nonlinearities at joints has been explored by Kim and Park (1994). The study by Duym and Schoukens (1996) consider the errors in force state map method arising out of measurement noise and incompleteness of model and propose a method for optimal points of grid points in fitting the force state map. The application of the force state mapping method to identify a system with double well potential has been investigated by Shin and Hammond (1998a). Experimental data from a beam oscillating under the influence of two symmetrically placed magnets are employed in this study and the authors investigate the performance of the method in its linear and chaotic response regimes. The same authors (Shin and Hammond 1998b), in a separate study, note the possibility of construction of phase portrait form a single measured state variable by using delay coordinates and the resulting pseudophase portrait is subsequently used in construction of force state map. An analysis of the errors in the parameter estimates obtained using force state mapping has been discussed by Meskell and Fitzpatrick (00). These authors employ phasor notations to investigate parameter identification problems associated with linear SDOF systems. The applications of force state mapping method to specific applications have been discussed by a few authors. Thus Meskell et al., (001) discuss the application of the method for identification of damping in fluid-structure interaction problems. The parameters associated with fluid-elastic instabilities in a triangular tube array have been investigated by Meskell and Fitzpatrick (003). Phani and Venkatraman (003) have used the method in their study on damping behavior of sandwich beam with electro-rheological fluid. One of the issues in the implementation of the force state mapping method has been the choice of functional representation of the restoring surface. The Chebychev polynomials proposed by Masri and Caughey (1979), and, used subsequently by several authors, are a set of orthogonal polynomials. Consequently, one need not re-estimate the coefficients if a lower order model is acceptable (Worden and Tomlinson, 001). However they obscure the meaning of the system parameters and their use is computationally intensive. The polynomial models discussed by Al-Hadid and Wright (1989, 1990, 199) are not best suited if the restoring surface is non-differentiableor discontinuous at a set of points. In fact, Worden and Tomlinson (001) note:...inthefamilyofpolynomial approximations to a given function over a given interval, there will be one which has the smallest maximum deviation from the function over the interval. This approximating polynomialthe minimax polynomial- has so far eluded discovery.... In this context it is of interest to note that in the area of computational mechanics that employs meshfree and particle methods, several functional representation schemes have been developed. Overviews on these developments are available in the works of Li and Liu (00) and Atluri (004,005). These methods include reproducing kernel particle method (RKPM), moving least squares (MLS) method, radial basis function (RBF), and multi-quadratic radial basis function (MQ-RBF) method. There exists a wide range of literature on the application of these methods for representing field variables in meshfree computational methods: see, for instance, the works of Liu 1995, Li and Liu 00, Liu 003, Atluri (004, 005), Choi and Marcozzi 004, Han et al., 005, Hon et al., 005, La Rocca et al., 005, Atluri et al., 006a,b, ie et al., 006, Sladek et al., 006, Shaw and Roy 007, Yuan et al., 007, Long et al., 008, Libre et al., 008, Ma 008, Orsini et al., 008, and Le et al., 008, A discussion on the possible mutual equivalence of these methods is also available (Li and Liu 00).Yet another class of methods for functional representation, known as kriging, is widely used in design analysis of computer experiments (Santner et al., 003). The present authors believe that these developments have important bearing on development of force state maps in the context of nonlinear system identification. With this in view, we consider, in the present paper, the problem of exploring the usefulness of RKPM and kriging

4 16 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 based methods for force state map construction. The RKPM methods have been developed in the context of representing field variables in meshfree computational methods (see, for instance, Liu 1995, Li and Liu 00, Liu 003, and Shaw and Roy 007). The RKPM based functions have the capability to represent polynomials of a given order exactly at any point within a specified domain. The method also possesses attributes such as capabilities for space or time localization, hp-like adaptivity, and multi-resolution analysis. More importantly, the methods have the capabilities to satisfactorily represent functions that have points of non-differentiability and (or) discontinuities within the domain of interest. To the best of the authors knowledge, the usefulness of RKPM based functional representations for the restoring surface in the context of the force state mapping method has not been explored in the existing literature. It is proposed in the present study to address this issue. Specifically the present study aims to answer the following questions:. How do the RKPM based methods perform in the context of force state mapping method? Can these methods be used to model nonlinearities that are either not differentiable or continuous at a set of points within the domain of interest?. Is qualitative analysis of nonlinear system possible with FSM? Can we determine fixed points and limit cycles of the system and investigate their stability using force state mapping method? (c). Using the force state mapping method, under what conditions would it be possible to determine the properties of the linear system that would result if nonlinearity were to be absent? (d). How do restoring force surface fit methods perform if nonlinearities are not known and are not of polynomial type? Another functional representation scheme which has met with success in engineering applications (Sacks et al., 1989, Santner et al., 003, Lophaven et al., 003, Kaymaz 005 and Panda and Manohar 008) but not explored in the context of force state map construction is the kriging based approximation. These representations adopt a probabilistic approach and the essential idea is to treat the function to be represented as nonstationary Gaussian random field with a postulated form of the nonstationary mean term and a homogeneous random field term having a specified covariance functional form. The functional form of the covariance function is taken to be known while its parameters are considered to be unknowns. Based on the available data these unknown parameters are estimated using the method of maximum likelihood estimation. In a subsequent step, the values of the function at points where its values are unknown are estimated using the criterion of minimization of mean square error of prediction. In the context of force state map construction, the advantage of application of kriging based approximation is that these methods are capable of taking into account the effect of the measurement noise in the observations. The present study explores the performance of kriging based force state map construction in nonlinear structural system identification. As a prelude to developing answers to these questions we first present a brief overview of functional representations using RKPM and kriging technique. Reproducing Kernel Particle Method (RKPM).1 Representation of smooth functions The exposition of the basics of RKPM has been provided by several authors: see for example, the papers by Liu et al., (1995a,b), Belytschko et al., 1996, Chen et al., (1996), Aluru, (000), Chen et al., (003), and Shaw and Roy (007). We briefly summarize here the essence of the procedure. For this purpose we consider a one-dimensional function u(x) which is sufficiently smooth and defined over x Ω R 1. Let us consider the transformation ν(x)=tu(x)= φ(x, s)u(s)ds (4) Ω

5 Force State Maps 17 Here T is a linear transformation, φ(x,s) is known as the kernel function, and ν(x) is the kernel approximation to the function u(x). If φ(x,s) = δ(s x), it follows that ν(x)=u(x) x Ω, that is, ν(x)= δ (x s)u(s)ds = u(x) (5) Ω When the condition φ(x s)=δ(x s) is not satisfied, the representation ν(x) of u(x) could suffer from phase and amplitude distortion and distortions near the boundary. Liu et al., (1995a, b) proposed a correction and introduced the kernel approximation as u a (x)= φ a (x;x s)u(s)ds (6) Ω with the modified kernel function φ a (x;x s) given by φ a (x;x s)=c(x;x s)φ a (x s) (7) where C(x;x s) is the correction function, to be determined; and φ a (x s) is the kernel function. Furthermore, the kernel function φ a (x s) is taken to be known and, typically, assumed to be given by a cubic spline function, that is, 0 z < (z+) 3 6 z 1 φ(z)= 3 z (+z) 1 z 0 3 (8) z ( z) 0 z 1 (z )3 6 1 z 0 z > with z = x s a. The parameter a appearing in this equation is termed as the dilational parameter. Its choice ensures that the support of φ a (x s) is compact. Figure illustrates typical plots of φ a (x s) for s = and for different values of a. The determination of the correction function C(x;x s) is based on the condition that u a (x) reproduces the function u(x) exactly at all x Ω whenever u(x) is a polynomial. If u(x) is not a polynomial, u a (x) is conditioned to reproduce Taylor sexpansion ofu(x) up to a specified order. Figure : Plots of kernel function φ a (x s) at s = for different values of the dilational parameter a. For computational purposes it is assumed that the domain Ω is discretized by a set of grid points {x i } with the associated nodal values of function u(x) given by u i = u(x i ); i = 1,,...,. The integral representation of equation (6) can now be replaced by a summation of the form (Liu et al., 1995a,b) u a (x)= C (x x i )φ a (x x i )u i (9) where { φ a (x x i )= 1 a φ ( )} x x i a is the known kernel function centered at x = x i. According to Liu et al., (1995a,b) the correction function C (x x i ) is taken as C (x x i )=b 0 (x)+(x x i )b 1 (x)+... +(x x i ) p b p (x) =H T (x x i )b(x) with H T (x x i )= { (x xi ) 0 (x x i ) 1 (x x i ) } p 1 ( p +1) being the set of monomial basis functions and b(x)= { b 0 (x) b 1 (x) b p (x) } T (10) being the ( p +1) 1 set of unknown functions to be determined. Equation (9) may now be written

6 18 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 as, u a (x)= ψ i (x)u i (11) The functions ψ i (x); i = 1,,..., are the RKPM shape functions. In terms of the correction and kernel functions, the shape functions are given by ψ i (x)=h T (x x i )b(x)φ a (x x i ) (1) To determine the unknown vector b(x), wenow impose the condition that u a (x) reproduces u(x) exactly up to a desired order as mentioned before. To develop this notion, we begin by consider the problem of exact reproduction of a constant. The requisite condition is clearly given by ψ i (x)x 0 i = 1 (13) It also follows that ψ i (x)(x x i ) 0 = 1 (14) Similarly, the reproduction of function x is given by ψ i (x)x i = x (15) and this can be re-written as x ψ i (x)(x x i )=0 (16) oting that ψ i (x)=1, the above equation can also be written as x ψ i (x) from which it follows ψ i (x)x i = 0 (17) ψ i (x)(x x i )=0 (18) Extending this logic to the reproduction of x, one gets ψ i (x)x i = x (19) This can be re-written as = = = = ψ i (x) ( x i x) ψ i (x) ( x i x +x ) ψ i (x)x i x ψ i (x)x i x ψ i (x)(x x i ) ψ i (x)+x ψ i (x) ψ i (x)x i +x ψ i (x) (0) Similarly, it can be shown that the reproduction condition for the function x p is given by ψ i (x)x p i = x p (1) leading to ψ i (x)(x x i ) p = 0 () Equations (14), (18), (0) and () can be written in a compact form as, ψ i (x)h (x x i )=H(0) (3) where H(0)= { } ( p +1) 1 ow using equations (1) and (3), we obtain H T (x x i )b(x)φ a (x x i )H (x x i )=H(0) (4)

7 Force State Maps 19 oting that the product H T (x x i )b(x)φ a (x x i ) is a scalar and b(x) is independent of the indexi, the above equation can be re-written as [ H (x x i )H T (x x i )φ a (x x i ) ] b(x)=h(0) (5) basis function H(x). It follows H 0 (x)=h(x)= { 1 x x x } p H 1 (x)= { 0 1 x p x p 1 }. H p (x)= { p! } (9) This can further be written in a compact form as M(x)b(x)=H(0) (6) where the ( p + 1) ( p + 1) matrix M(x) = H (x x i )H T (x x i )φ a (x x i ) is called the moment matrix. Consequently, the unknown vector b(x) is obtained as b(x)=m 1 (x)h(0) (7) Using equations (1) and (7), the RKPM shape function, ψ i (x) can be written as ψ i (x)=h T (x x i )M 1 (x)h(0)φ a (x x i ) (8) The determination of b(x) requires the inversion of the moment matrix M(x) for every value of x that is of interest.. RKPM with reproduction of derivatives The RKPM formulation can be extended to achieve reproduction of derivatives up to a desired order and these issues have been studied by a few authors (Chen et al., 1996, Jun et al., 1998, Aluru 000, and Shaw and Roy 007). The scheme for derivative reproduction of a given function u(x) is based on the assumption that the β th derivative of the RKPM representation u a (x) will exactly reproduce the β th derivative of u(x) whenever u(x) is a polynomial of degree p. Moreover, the RKPM representation for the derivative is obtained in terms of the nodal values u i = u(x i ); i = 1,,..., of the function u(x). The method proposed by Shaw and Roy (007) has been demonstrated to be computationally advantageous and we follow this procedure in the present study. Let H β (x) represent the β th derivative of monomial For the purpose of illustration we restrict the discussion to the case of β = 1. We derive the moment equation for this case and then subsequently extend the formulation to the β th order derivative reproduction of a given function. For this purpose we introduce a new family of RKPM shape functions denoted by θ β i (x) which reproduce the β th order derivative of the function u(x). We consider the case of the function to be a constant and clearly the condition for the reproduction of the first derivative by the RKPM representation is given by θi 1 (x)=0 (30) where θ 1 i (x) represents the derivative of θ i(x) with respect to x. The above equation can also be written as θ 1 i (x)(x x i ) 0 = 0. (31) Similarly, the derivative reproduction condition for a function x, isgivenby θ 1 i (x)x i = 1 (3) oting that be re-written as x θ 1 i (x) from which it follows θi 1 (x)=0, the above equation can θ 1 i (x)x i = 1 (33) θ 1 i (x)(x x i)= 1 (34)

8 130 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 Similarly, the condition for reproduction of first derivative of the function x leads θ 1 i (x)x i = x (35) The above equation can be re-written as θ 1 i (x)x i x = x θ 1 i (x) x θ 1 = = θ 1 i (x)x θ 1 θ 1 i (x)(x x i ) i (x)x i + i (x)xx i + θ 1 i (x)x i θ 1 i (x)x i (36) The above formulation can be generalized to show that the derivative reproduction condition for the function x p is given by θi 1 (x)x p i = p x p 1 leading to (37) θ 1 i (x)(x x i ) p = 0 (38) Equations (31), (34), (36) and (38) can be written in a compact form as θ 1 i (x)h (x x i )=( 1)H 1 (0) (39) where H 1 (0)= { } ( p +1) 1 Clearly, in equation (39) θi 1 (x) can be written in terms of corrected kernel function as follows, θi 1(x)= C 1 (x x i )φ a (x x i ) (40) where C 1 (x x i )=H T (x x i ) b(x) is the correction function for the first-order derivative reproduction. ow, using equations (39) and (40), one gets, H T (x x i ) b(x)φ a (x x i )H (x x i ) =( 1)H 1 (0) (41) where b(x) = { b 0 (x) b 1 (x) b p (x) } T is the vector of unknown coefficients. Equation (41) may be written as, [ H (x x i )H T (x x i )φ a (x x i ) ] b(x) =( 1)H 1 (0) (4) This can be written in a compact form as M(x) b(x)=( 1)H 1 (0) (43) from which it follows b(x)=( 1)M 1 (x)h 1 (0) (44) Using equations (40) and (44), the shape function for the first-order derivative reproduction of an arbitrary function can be written as, θi 1 (x)= ( 1)H T (x x i )M 1 (x)h 1 (0)φ a (x x i ) (45) The shape function for β th order derivative reproduction of an arbitrary function using similar procedure as we described above, one can write as θ β i (x)= ( 1) β H T (x x i )M 1 (x)h β (0)φ a (x x i ) (46) Clearly, the above formulation does not require the differentiation of the kernel function φ a (x) and requires the inversion of moment matrix M(x) only once..3 Representation of non-smooth functions In the context of developing force state maps using RKPM representations, it is crucial to be able to obtain satisfactory representations for maps with points of discontinuities and (or) points of non-differentiabilities. Thus, for example, in modeling energy dissipation due to friction, representation for signum function, that is, functions of the form sgn(x a), need to be obtained; similarly, in modeling piecewise linear systems, representation for modulus functions of the form x a

9 Force State Maps 131 becomes necessary. The RKPM provides powerful tools to achieve these difficult goals. Thus, the problem of RKPM models for functions (in more than one variables) which possess either points of discontinuities and (or) points at which the function is non-differentiable in one or more independent variables, has been studied by a few authors (Krongauz et al., 1998, Zhang et al., 004, Lu et al., 005, Joyot et al., 006 and Shaw and Roy, 007). The essence of the idea here is to enrich the RKPM representation by deliberately introducing discontinuities so that discontinuity either in the function or its derivative(s) could be modeled. We illustrate the idea by considering a domain Ω and a function u(x) with x Ω such that the domain Ω is divided in to two sub-domains Ω 1 and Ω such that Ω 1 Ω = Ω and Ω 1 Ω = Φ where Φ = null set (Figure 3). For the purpose of simplicity we assume Ω R 1 and define {u i } to be the nodal values of u(x) at x = x i ; i = 1,,...,. Furthermore, we group the nodal values into two groups, namely, { u I i (x)} I with x Ω 1 and { u II i (x) } II i =1 with x Ω such that I + II =. The nodal co-ordinates are also grouped as x I i and x II i such that x I i Ω 1 and x II i Ω with i = 1,,..., I and i = 1,,..., II. We wish to obtain a RKPM based model for u(x) which allows for the presence of discontinuity in u(x). To achieve this, we write the kernel approximation u a (x) as u a (x)= = I II i =1 ψ I i (x)ui i if x Ω 1 ψ i II (x)u II i if x Ω (47) Furthermore, using the indicator function I Ωi (x), given by I Ωi (x)=1 ifx Ω i (48) = 0 ifx / Ω i and, also noting that, ψi I (x)=h T ( x x I ) i b I ( ) (x)φ a x x I i ψ II i (x)=h T ( x x II i ) b II ( ) (49) (x)φ a x x II i we obtain 1 L d Figure 3: The division of the domain Ω into Ω 1 and Ω which are separated by a line of discontinuity. u a (x)= [ I H T ( ] x x I ) i b I ( ) (x)φ a x x I i u I i I Ω1 (x) + [ II i =1 H T ( x x II i ) b II ( ) (x)φ a x x II i u II i ] I Ω (x) (50) It is clear that, the RKPM representation has been enriched with functions that allow for possible existence of discontinuity in u(x) within x Ω. The determination of the functions b I (x) and b II (x) now follows exactly the same procedure as has been described in the previous section with the provision that b I (x) are determined based on the knowledge of u(x) at x i = Ω 1 ; i = 1,,..., I and b II (x) based on u(x) at x i = Ω ; i = 1,,..., II. The details of this procedure are best illustrated by considering a few examples. Example 1: Consider f (x)=0.5sgn(x 6) with 4 x 7. Here f (x) has a discontinuity at x = 6. Accordingly, we divide 4 x 7 into two intervals Ω 1 :4 x < 6andΩ :6 x 7. ote that, x = 6 is included in Ω and excluded in Ω 1.Figure 4 shows the RKPM shape functions valid for Ω 1 and Ω respectively. Figure 5 compares the RKPM representations of f (x) using enrichment

10 13 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 Figure 4: Shows RKPM shape functions ψ i (x) for the function f (x) at and near the point of discontinuity. Figure 5: RKPM representation of a function f (x); shows enlarged view of a fig. 5 at the point of discontinuity x = 6. and without using enrichment. In the latter case the RKPM representation can be observed to be unsatisfactory. Example : Here we consider f (x)=0.8 x 6 with 4 x 7. Here f (x) is continuous every where in 4 x 7 but its derivative has a discontinuity at x = 6. To implement RKPM we define Ω 1 :4 x < 6andΩ :6 x 7. Within Ω 1 and Ω we construct RKPM approximation to f (x) with provisions for exact reproduction of the first derivative. Figure 6 shows the RKPM shape functions ψ i (x)&θi 1 (x) and the resulting representations are shown in Figure 7. The superior performance of the enriched RKPM is again evident from the Figures. Example 3: Here we consider f (x)=0.5sgn(x 5)+0.8 x 6 with 4 x 7. This function has discontinuity at x = 5 and a point of nondifferentiability at x = 6. We define three regions, Ω 1 :4 x < 5, Ω :5 x < 6andΩ 3 :6 x 7 and represent f (x) in three regions as f (x)= (x 6); 4 x < 5 = (x 6); 5 x < 6 = (x 6); 6 x 7 (51) The RKPM representations obtained are shown in Figure 8. The RKPM with enrichment is shown to perform better especially around x = 5andx =

11 Force State Maps 133 (c) (d) Figure 6: - Shows RKPM shape functions ψ i (x) for the function f (x) at and near the point of nondifferentiability; (c)-(d) shows RKPM derivative shape functions θi 1 for the function f (x) at and near the point of discontinuity. 6, where the RKPM (without enrichment) show relatively large errors. In obtaining these results it is assumed that u(x) is available at an interval of Δx = The dilational parameter a = 6is considered for the above case (explanation for this choice is provided in later section). Remark: The generalization of the formulations in the preceding subsections (equations 8, 46 and 50) for the case of vector of multi-dimensional functions is conceptually straightforward and we do not present these details here. The details can be found in the thesis by amdeo (007). 3 Functional representation using kriging Here we begin by considering a function g(x) defined with respect to a l-dimensionalinput vector x. Let { x 1 x x n } T be such that xi R l (i = 1,, n) denote the n l matrix of input points at which the output g(x i ), i = 1,, n are available. In the context of force state map construction, x represents the state vector consisting of displacement and velocity components, while g(x) represents the restoring force (as in equation ). We model the function using the representa-

12 134 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 (c) (d) Figure 7: RKPM representation of a function f (x); shows enlarged view of fig. 7, having nondifferentiability at x = 6; (c) shows RKPM derivative reproduction of f (x); (d) shows enlarged view of fig. 7(c). tion g(x i )= p f j (x i )β j +Z(x i )= f T (x i )β +Z(x i ); j=1 i = 1,,...n (5) Here the superscript T denotes matrix transposition, { f j (x i ) } p is a set of known functions, { } j=1 p β j a set of unknown constants to be determined and Z(x) is a homogenous, zero mean, j=1 Gaussian random field evolving in x with an unknown covariance function. Thus, for example, with l=3 and quadratic basis functions, at the i th sampling point we get [ ] g x (1) i,x () i,x (3) i = β 1 +β x (1) i +β 3 x () i +β 4 x (3) i ( ) +β 5 x (1) ( ) ( ) i +β6 x (1) i x () i +β 7 x (1) i x (3) i ( ) +β 8 x () ( ) ( ) i +β9 x () i x (3) i +β 10 x (3) i [ ( ) ] + Z x (1) i,x () i,x (3) T i (53) so that, in this approximation, p = 1+l + l(l+1) = 10. The covariance function is taken to be of the form Cov[Z (x r )Z (x s )] = σ R(x r x s,θ) (54)

13 Force State Maps 135 Figure 8: RKPM representation of a function f (x); shows enlarged view of a fig. 8 at the point of discontinuity x = 5; (c) shows enlarged view of a fig. 8 at the point of non-differentiability x = 6. (c)

14 136 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 Here σ is the unknown variance of the random field Z (x i ) and R(x r x s,θ) is the (r,s) th element of the n n matrix of correlation coefficients with θ being a q 1 vector of unknown parameters in the covariance model. The functional form of the correlation matrix R(x r x s,θ) is taken to be known but the function itself is considered to be dependent on the unknown parameters θ. Some of the choices for the functional form of R(x r x s,θ) employed in the existing literature include the following (Sacks, et al., 1989, Lophaven et al., 003) R(θ rs,x r x s )=exp [ θ rs x r x s ] (55) [ ] R(θ rs,x r x s )=exp θ rs x r x s θ n+1 ; (56) 0 < θ n+1 < [ R(θ rs,x r x s )=exp θ rs (x r x s ) ] (57) R(θ rs,x r x s )=max {0,1 θ rs x r x s } (58) R(θ rs,x r x s )=1 1.5ξ rs +0.5ξrs; 3 (59) ξ rs = min(1,θ rs x r x s ) R(θ rs,x r x s )=1 3ξrs +ξrs; 3 (60) ξ rs = min(1,θ rs x r x s ) R(θ rs,x r x s )=ς(ξ rs ); (61) 1 15ξrs +30ξrs 3 for 0 ξ rs 0. ς (ξ rs )= 1.5(1 ξ rs ) 3 for 0. ξ rs < 1 0 forξ rs 1 The work of Santner et al., (003) discusses the issues related to smoothness of the functional representation in equation (5) vis-à-vis the choice of the form of R(x r x s,θ) and values of parameters θ. In the present study we employ the third of the above correlation coefficient model. It may be noted that the function g( ), as in equation (5), constitutes a non-homogeneous Gaussian random field with mean E [g(x)] = f T (x)β and covariance E [{ g(x r ) f T (x r )β }{ g(x s ) f T (x s )β }] = σ R(x r x s,θ). The non-homogeneity of the random field here arises due to the dependence of mean value of g(x) on the parameter x. Thus, the set of unknowns to be determined in equation (5) are { } p β j j=1, σ and {θ i } q. Corresponding to the points x = { } T, x 1 x x n we can re-write equation (5) as g(x)= F(x)β + Z(x) (6) Here g(x) is a n 1 vector of random variables, F(x)is a n p matrix of deterministic functions with (i, j) th element given by F ij (x)= f j (x i ) and Z(x) is a n 1 vector of zero mean Gaussian random variables. Thus the probability distribution function (PDF) of g(x) is given by g(x)= [ F(x)β,σ R(x r x s,θ) ] (63) where denotes the normal PDF with mean vector F(x)β and covariance matrix σ R(x r x s,θ). In order to determine the unknown model parameters { } p β j j=1, σ and θ we invoke the method of maximum likelihood estimation. Accordingly, we construct the negative log-likelihood function L ( β,σ,θ ) = 1 { nlnσ +ln R } +(g F(x)β ) t 1 σ R 1 (g F(x)β)+nlnπ (64) and minimize this function with respect to the parameters { } p β j j=1, σ and {θ i } q. The minimization with respect to { } p β j can be shown to lead j=1 to the condition ˆβ = ( F T R 1 F ) 1 F T R 1 g (65) Similarly, the minimization with respect to σ leads to ˆσ = 1 n (g Fβ )T R 1 (g Fβ ) (66) ow by substituting equations (65) and (66) into equation (64) we get the objective function only in terms of parameter vector θ as ) L( ˆβ, ˆσ,θ = 1 { nln ˆσ (θ)+ln R(θ) +nlnπ } (67) The minimization of this function, with respect to the remaining variables {θ i } q, can now be

15 Force State Maps 137 carried out using numerical optimization methods. The conditions θ i > 0; i = 1,,,q are imposed as constraints in this optimization problem. These constraints are needed given that the correlation function that we are using is as given by the third of the models listed in equation (57). A penalty function to artificially enforce very high value of the given function is used to avoid negative values of θ. The resulting predictor g(x) with β = ˆβ ( ˆθ ), σ = ˆσ ( ˆθ ) and θ = ˆθ is designated as the maximum likelihood empirical best linear unbiased predictor for the force state map. At this stage we have determined the values of all the unknowns in the model in equation (63). ow we consider the question of predicting the output at a point x 0 at which the function g(x) has not been measured. 3.1 Prediction of g(x) at an unmeasured point Let x 0 be a point in the input space at which we now aim to estimate g(x 0 ) and from equation (5), we get g(x 0 )= f T (x 0 )β + Z(x 0 ). It is clear that, g(x 0 ) is a Gaussian random variable and the vector { g(x 0 ) g(x) } T forms a (n +1) 1 vector of Gaussian random variables with PDF given by { g(x0 ) g(x) } = ([ f T (x 0 ) F(x) ] ]) ˆβ ( ˆθ), ˆσ ( ˆθ)[ 1 r T 0 ( ˆθ) r 0 ( ˆθ) Rr 0 ( ˆθ) (68) Here, r 0 = [ R(x 0 x 1 ) R(x 0 x n ) ] T.The estimate ĝ(x 0 ){ of g(x 0 ) that minimizes the mean square error E [g(x 0 ) ĝ(x 0 )] } is well known to be given by (Papoulis and Pillai, 001) ĝ(x 0 )=E [g(x 0 ) g(x)] = f T (x 0 ) ˆβ ( ˆθ ) +r T 0 ( ) ˆθ R 1 ( ˆθ )[ g(x) F ˆβ ( (69) )] ˆθ This expression constitutes the kriging model and this forms the model for the force state map in characterizing the nonlinear dynamical system. In implementing this model we still need to choose the functions { f j (x) } p and, in the present study, j=1 we take quadratic functions (as illustrated in equation 53 for l=3) so that p = 1 +l + l(l+1).ifcross terms are ignored, one getsp = l RKPM for force state map construction As has been mentioned already, the objective of the present investigation is to explore the performance of RKPM based models for force state maps. We consider the following cases which are sequenced in the order of increasing complexity (in terms of number of dofs and treatment of discontinuities and non-differentiable restoring force surface). Case 1: Here we consider SDOF systems with restoring force surface being a smooth function of displacement and velocity variables. This would mean that f [x(t),ẋ(t)] in equation () is continuous and differentiable in x and ẋ over the domain of interest. The construction of RKPM based force state map involves the following steps: 1. Estimate mass m.. Measure g(t) and ẍ(t) at t i = iδt; i = 1,,..., and derive f [x(t),ẋ(t)] = g(t) mẍ(t) at t i = iδt and i = 1,,...,. 3. Measure x(t) and ẋ(t) at t i = iδt; i = 1,,...,. It may be noted that the evaluation of f [x(t), ẋ(t)] in the preceding step requires already the measurement of ẍ(t). This would mean that at every dof of interest it is required to measure x(t), ẋ(t) and ẍ(t). This may not be often possible. In case only one of the quantities x(t) and ẋ(t) is measured, one could estimate the unmeasured quantity using Kalman filter Crawley and O Donnell (1987). Thus for the purpose of illustration, if we assume that it has been possible to measure only ẋ(t) and ẍ(t), an estimate of x(t) can be obtained by considering the process and measurement equations given respectively by [ẋ ] [ ][ ] [ ] [ ] 0 1 x 0 0 = + ẍ ẍ 0 0 ẋ 1 m + ξ 1 y = ẋ m = [ 0 1 ][ ] (70) x + χ ẋ

16 138 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 where ẍ m denote measured acceleration, ẋ m denote the measured velocity, ξ represents noise associated with the measured acceleration, and χ represents noise associated with the measured velocity. The above equation can be written in the standard state space form as, Ẋ = AX +Bu d +Bξ y = HX + χ where X = { x (71) ẋ } T, ud = ẍ m (driving input), H = [ 0 1 ], B = [ 0 1 ] T [ ] T, B = 0 1, [ ] 0 1 and A =. One could now apply the 0 0 Kalman filter to carry out the state estimation problem leading to the estimate of x(t). 4. The above measurements lead to the discretized values x i,ẋ i and f i (x i,ẋ i ) for i = 1,,...,. Plotf i (x i,ẋ i ) versus x i &ẋ i. 5. The points (x i,ẋ i ) define the nodes in the RKPM representation and with the corresponding f (x i,ẋ i ) being the nodal-values of the restoring force. It should be noted that the points (x i,ẋ i ); i = 1,,..., would not be equi-spaced in the x ẋ plane. In the implementation of RKPM we divide the range of values of x and ẋ into a set of equi-spaced grid points and determine the approximation f a (x,ẋ) at these grid points. The value of the dilational parameter a needs to be selected at every grid point so that an acceptable approximation to f (x i,ẋ i ) is obtained locally at every grid point. To achieve this, at any grid point, n number of particles in the neighborhood are selected where n = (n d + p )! (7) n d! p! Here n d = dimension of the problem (= for the SDOF systems) and p = order of the polynomial (to be chosen). Thus, if n d = and p = 3, one gets n = 10. At any given grid point, the parameter a is selected such that a = the distance from the grid to the farthest among n particles in the neighborhood of the grid point. Once the parameter a is chosen the RKPM shape functions are evaluated at the specified grid point by inverting the moment matrix M(x) (equation 8). This exercise needs to be repeated at all the chosen grid points. In the numerical work it was observed that the inversion of M (x) at certain grid points became erroneous due to poor conditioning of M(x). A remedy to this difficulty lied in increasing n in the evaluation of a, and, hence, the size of matrix M(x) increases. 6. The representation in the previous step provides the complete solution in the sense that the RKPM based force state map f a [x(t),ẋ(t)] could be used to predict numerically the system response to any excitation different from g(t) for any prescribed initial conditions x(0) and ẋ(0). 7. The fixed points of the dynamical system (equation 1) given by ẋ = 0and f (x,0)=0. Thus by constructing f a (x,ẋ = 0), one could determine the fixed points by determining the zeros ] given by f a (x,ẋ = 0) =0. Let [x j,0 ; j = 1,,...denote there fixed points. The stability of the fixed points could be examined by constructing the Jacobian matrix [ ] 0 1 J j = f x 1 f dx (x j,0) (73) and determining the eigenvalues of this matrix. As is well known, if the real points of the eigenvalues are nonzero, the linear stability analysis is conclusive (Simmons, 197). It is of interest to note that the gradients f x and f ẋ required in the above analysis can be accurately established by employing RKPM with derivative reproduction at ẋ = 0and x = x j ; j = 1,, The knowledge of f a (x,ẋ) could also be employed to study the possible existence of limit cycles. This could be achieved by employing some of the existing theorems for verification of existence (or otherwise) of

17 Force State Maps 139 limit cycles or by numerical simulations of free vibration trajectories. 9. If a parametric form of f (x, ẋ) isknown,and is taken to be a polynomial, the coefficients of this polynomial could be established by examining the gradients of f (x,ẋ) typically at origin. Here again, the gradient reproduction capabilities of RKPM representations offer powerful means to arrive at accurate estimates of unknown parameters. 10. By linearizing the function f a (x,ẋ) around the origin, it is possible to establish properties of the linear system that would result if nonlinearities were to be absent. The procedure described in this sub-section is illustrated by considering the following oscillators. 1. Linear SDOF system: mẍ +cẋ +kx = g(t); x(0)=0, ẋ(0)=0 (74). Duffing s oscillator: mẍ +cẋ +kx+αx 3 = g(t); x(0)=0, ẋ(0)=0 (75) 3. Duffing s oscillator with negative linear stiffness: mẍ +cẋ kx+αx 3 = g(t); x(0)=0, ẋ(0)=0 (76) 4. Van der Pol s oscillator: mẍ εẋ ( 1 x ) +kx = g(t); x(0)=0, ẋ(0)=0 (77) It may be noted that in all these examples, the restoring force is a polynomial in x and ẋ. In all these examples it is assumed that m = 1kg, k = 9/m and g(t)=p s=1 (a s cosω s t +b s sinω s t) (78) (c) Figure 9: Studies on linear SDOF system; phase plane plot; restoring force surface based on RKPM; (c) plot of f [x,ẋ = 0]. with P = 10, ω s = sω 0 ; s = 1,,...,, = 15, ω 0 =.4rad/s. The parameters a s and b s,

18 140 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 (c) (e) (d) (f) Figure 10: Studies on Duffing s oscillator; phase plane plot; restoring force surface based on RKPM with p = 3; (c) plot of f [x,ẋ = 0]; (d) restoring force surface based on RKPM with p = 1; (e) restoring force surface based on RKPM p = ; (f) restoring force surface based on RKPM with p = 4.

19 Force State Maps 141 s = 1,,..., are taken to be realizations of zero mean, Gaussian, mutually independent distributed random variables with unit standard deviations. The measurements data are generated synthetically by using the ODE 45 routine of Matlab with relative tolerance reltol = 1.0e-1 and absolute tolerance abstol = 1.0e-1. The program outputs the results at a uniform time interval of Δt = π 80ω nat s, where ω nat = k/mrad/s. The results on the force map models obtained using RKPM are shown in Figures 9-1. Figure 9 shows the results for the linear SDOF system with c = 0.1s/m. The measured response is depicted in Figure 9 in the form of a phase plane plot and the RKPM based force state map is shown in Figure 9. The coefficients f a x and f a ẋ are evaluated at x = 0&ẋ = 0, which represents the values for k and c respectively. The values of k and c are obtained as k = 9.00/m and c = 0.1s/m and this matches very well with the corresponding reference values. The plot of f a [x,ẋ = 0] versus x is shown in Figure 9(c) from which the linear dependence for spring force on displacement is inferred. Also, from Figure 9(c) it is concluded that x = 0 and ẋ = 0 is a fixed point and it is verified that this point is stable. Results shown in Figure 10-(f) correspond to the case of Duffing s oscillator with α = 10 /m 3 and m,c,k as in the previous example. Based on the measurements shown in Figure 10, the force state map as in 10 is obtained which, in turn, leads to the coefficients c = f a ẋ, x=0,ẋ=0 k = f a x,and6α = 3 f a x=0,ẋ=0.thevalues of k = 9.00 /m, c = 0.1s/m, and α = x=0,ẋ=0 x /m 3 so found, compare very well with the corresponding reference values of k = 9/m, c = 0.1s/m, and α = 10/m 3. Figure 10(c) show the plot of f a [x,ẋ = 0] versus x with p = 3, from which it emerges that x = 0 is the fixed point which is further verified to be stable. It must be noted that with p = 3, the RKPM representation, in principle, permits the presence of quadratic nonlinear terms and also cubic terms involving products of x(t) and ẋ(t) and the term ẋ 3 (t) in the restoring force surface representation. It has been verified in the numerical work that the coefficients (c) Figure 11: Studies on Duffing s oscillator with negative linear stiffness; phase plane plot; restoring force surface based on RKPM; (c) plot of f [x,ẋ = 0].

20 14 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 (c) Figure 1: Studies on Van der Pol s oscillator; phase plane plot; restoring force surface based on RKPM; (c) free vibration characteristic obtained from the original system and from the identified system; full line: results from the reference system; cycle line: results from identified system. of these terms turn out to be close to zero. It is also to be noted that the choice of p = 3 would not be apparent at the outset. One could begin with p = 1 and increase its value progressively. This, in fact has been done in the present case: Figures 10(d)-(f) show the RKPM force state maps for the cases of p = 1,, and 4 respectively. Based on the examination of errors associated with each of these representations it may be concluded that the choice of p = 3 is the most appropriate for this example. Results from studies on the Duffing s oscillator with negative linear stiffness (equation 76) are depicted in Figure 11 -(c). The fact that this system has an unstable fixed point at the origin, and two stable points at x = 0.88&0.89m and ẋ = 0, is borne out by the RKPM based force state map (Figure 11c). Here, again, the model for the force state map leads to the estimates k = 9.00/m, c = 0.1s/m, and α = 10.00/m 3, which agree well with the corresponding reference values of k = 9/m, c = 0.1s/m, and α = 10/m 3. Figure 1 shows the results of studies on the Van der Pol s oscillator with ε = 10. Based on the phase plane plot shown in Figure 1 and the RKPM based force state map shown in Figure 1, it is not obvious that the system possesses a limit cycle. However, a stability analysis of the fixed point reveals the presence of an unstable fixed point at the origin. A subsequent simulation of free oscillations (Figure 1c) indeed reveals the presence of a stable limit cycle. It may also be noted that once the model f a (x,ẋ) for the force state map is obtained, it is possible to apply criterion such as Poincare-Bendixson s criterion to investigate the possible existence of limit cycle. The values of k = 9.00/m, and ε = so found, compare very well with the corresponding reference values of k = 9/m and ε = Case : Linear MDOF systems Here we assume that the vibrating system under study can be modeled as a n-dof linear system governed by the equation MẌ +CẊ +KX = g(t); X (0)=X 0, Ẋ(0)=Ẋ 0 (79)

21 Force State Maps 143 Based on the assumption that the n n mass matrix M is known, it becomes possible to express the n 1 vector of force state maps as f [ X(t),Ẋ(t) ] = g(t) MẌ(t) (80) Each element of is expected to be functions of the n 1 variables X(t) and Ẋ(t). The implementation of force state map thus requires RKPM models for in n-dimensions. This obviously is a computationally demanding task. Also, the study would need elaborate measurements to be made on elements of vectors X(t) and Ẋ(t). Clearly, the elements C ij and K ij of damping matrix C and stiffness matrix K are respectively given by C ij = f [ ] i X,Ẋ & K Ẋ j ij = f [ ] i X,Ẋ X=0 X j X=0 Ẋ=0 Ẋ=0 (81) If the focus of study is limited to determination of C ij andk ij, it is enough to obtain RKPM fit to the force state map only near the origin. This would reduce the computational effort involved in parameter identification. For the purpose of illustration we consider the frame structure shown in Figure 13. This frame structure consists of a rigid slab mounted on three aluminum columns and a steel column. This would mean that the structure has planar asymmetry and under the action of horizontal base motion, the frame not only translates along the horizontal axes but also twists about the vertical axis. This structure has been studied experimentally by Manohar and Venkatesha (006). Based on the free body diagram shown in Figure 14, the equation of motion for this structure can be obtained as m z +k 1 (z v 1 θ)+c 1 (ż v 1 θ)+k 4 (z v 1 θ) +c 4 (ż v 1 θ)+k 8 (z +v θ)+c 8 (ż v θ) +k 5 (z +v θ)+c 5 (ż+v θ)= mẍ s (t) mÿ +k (y +u 1 θ)+c (ẏ +u 1 θ)+k 3 (y u θ) +c 3 (ẏ u θ)+k 7 (y +u 1 θ)+c 7 (ẏ +u 1 θ) +k 6 (y u θ)+c 6 (ẏ u θ)=0 [ I θ +u 1 (k +k 7 )(y +u 1 θ)+(c +c 7 )(ẏ+u 1 θ) ] [ u (k3 +k 6 )(y u θ)+(c 3 +c 6 )(ẏ u θ) ] [ +v (k5 +k 8 )(z +v θ)+(c 5 +c 8 )(ż+v θ) ] [ v 1 (k1 +k 4 )(z v 1 θ)+(c 1 +c 4 )(ż v 1 θ) ] = 0 z(0)=0; y(0)=0; θ(0)=0; ż(0)=0; ẏ(0)=0; θ(0)=0 (8) where k 1 = k = k 3 = k 4 = k 7 = k 8 = K a represents the stiffness coefficients for the aluminum columns and k 5 = k 6 = K s represents the stiffness coefficients for the steel column. The stiffness coefficients for the aluminum and steel columns are given by K a = 1E ai a h 3 ; K s = 1E si s c h 3 (83) c where E a I a = flexural rigidity of aluminum column and E s I s = flexural rigidity of steel column. The location of mass center of the frame structure (see Figure 15) is evaluated using the relation as given by z = m S b s +(m 1 +m 4 ) y = m S d s +(m 4 +m 3 ) ( ) bs b s1 ) +(m +m 3 )( bs +b s1 ρ s tb s d s +m 1 +m +m 3 +m 4 ( ) ds d s1 ) +(m 1 +m )( ds +d s1 ρ s tb s d s +m 1 +m +m 3 +m 4 m 1 = m = m 3 = A a ρ a (h c /) m 4 = A s ρ s (h c /), m S = ρ s tb s d s (84) where ρ s = mass density of steel, ρ a = mass density of aluminum, A a,a s = areas of cross-section of aluminum and steel columns respectively, and h c = height of the columns. The torsional stiffness of the frame structure due to shear rigidity of individual columns can be written as k = 3G aj a h c + G sj s h c ; J a = πd4 a 3 ; J s = πd4 s 3 (85)

22 144 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 Figure 13: Details of the frame structure made up of a steel slab supported on one steel column and three aluminum columns. where G a J a = torsional rigidity of aluminum column, G s J s = torsional rigidity of steel column, and D a, D s = diameter of the aluminum and steel columns respectively. The mass m of the rigid slab and mass moment of inertia I of the frame structure is given by m =m slab +m columns =ρ s tb s d s +3 ρ aa a h c I = ρ stb s d s ( b 1 s +ds ) +ρ s tb s d s [ ( z b s + ρ sa s h c ) +( y d ) ] s + ρ aa a h c (u 1 +v 1 +u +v 1 +u 1 +v ) + ρ sa s h c ( u +v ) +3 ρ aa a h c (r a ) + ρ sa s h c (r s ) (86) Here r a,r s = radii of the aluminum and steel columns respectively. It may be noted that the mass moment of inertia I is being computed with respect to the mass center of the frame structure. For the numerical analysis, we compute the damping matrix C using the modal damping ratios as given by the relation, C = [ ϑ t] 1 [ ] 1 [κ] ϑ where, κ is a diagonal matrix with entry on the n th row being η n ω n and ϑ is the modal matrix for the system. The forcing function x s (t) is taken to be a single frequency sinusoidal excitation given by x s (t)=δcosλt (87) where Δ is taken to be an amplitude of support motion and λ is a frequency of excitation. For the reference system, we take λ = π f with f = 5.86 Hz, Δ = m, E a = /m, E s = /m, G a = /m, G s = /m, ρ a = kg/m 3, and ρ s = kg/m 3. The measurements data are generated synthetically for 0 t f 15s with a time increment of Δt = π 80ω n3 s, where ω n3 =

23 Force State Maps 145 k 1 c k c k 3 3 k 4 x (t) s c 1 c 4 x (t) s v υ z x (t) s k 8 v 1 k 5 x (t) s c u 8 1 c 7 k 7 y u c 6 k 6 c 5 k (y + u υ) c (y + u υ) 1 k (y - u υ) c 3 (y - u υ) k (z - v υ - x ) s... c (z - v υ - x ) 1 1 s k (z - v υ - x ) s... c (z - v υ - x ) 4 1 s k (z + v υ - x ) + 8 s... c (z + v υ - x ) 8 s.. Iυ.. my.. mz k (z + v υ - x ) + 5 s... c (z + v υ - x ) 5 s k (y + u υ) + k (y - u υ) c 7 (y + u 1 υ) c 6 (y - u υ) Figure 14: Rigid mass-damper-spring model representation of the frame structure subjected to harmonic base motion x s (t); free body diagram showing the forces acting on the frame.

24 146 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 Aluminum bs Aluminum A B m 1 m ds / 1 ds z ds / 1 D y m m 3 4 Aluminum (0,0) bs / 1 bs / 1 C Steel Figure 15: Location of the mass center. the third natural frequency of the system. ow, for the identification of system parameters we use equation (81), to estimate the elements of damping matrix C and stiffness matrixk. In doing so, we utilize the derivative reproduction capabilities of RKPM based force state maps. It is observed that the identified damping matrix C and stiffness matrix K compare very well with the corresponding reference matrices as shown below: (units of [K] is /m and [ C ] is s/m) K re f. = ; K est. = C re f. = ; C est. = Case 3: MDOF systems with smooth nonlinearities The procedure described in the preceding subsection can be extended in a straight forward manner. To illustrate this we consider a -dof system with cubic nonlinearity (see Figure 16 assuming μ = 0): m 1 ẍ 1 +c 1 ẋ 1 +c (ẋ 1 ẋ )+k 1 x 1 +k (x 1 x ) +αx 3 1 = g 1 (t) m ẍ +c (ẋ ẋ 1 )+k (x x 1 )=g (t) x 1 (0)=0; x (0)=0; ẋ 1 (0)=0; ẋ (0)=0 (88) k, α 1 1 μ c 1 m 1 x (t) k g (t) 1 c m x (t) g (t) Figure 16: -dof system having cubic nonlinearity α and frictional damping μ. The system parameters considered are m 1 = m = 1kg, k 1 = 10/m, k = 40/m; the damping

25 Force State Maps 147 parameters are selected in such a way that, when α = 0, the two modes of oscillations have damping ratios η 1 = η = 0.0. The excitations g 1 (t) and g (t) are taken to be similar to g(t) as in examples considered in section 4.0 with ω 0 = 0.1ω n1,whereω n1 = the first natural frequency of the system when α = 0. The measurements data are synthetically generated for 0 t f 15s with Δt = π 80ω n s, where ω n = the second natural frequency of the system when α = 0. If the problem of system identification is now viewed within the framework of parametric identification methods, the task is now to determine c 1,c,k 1,k and α from the RKPM based force state maps fi a [x,ẋ]; i = 1,. Thus, one gets c 11 = c 1 +c = f 1 a ẋ 1 c 1 = c = f 1 a ẋ c 1 = c = f a ẋ 1 c = c = f a ẋ k 11 = k 1 +k = f 1 a x 1 k 1 = k = f 1 a x k 1 = k = f a x 1 k = k = f a x α = 1 3 f a 1 6 x 3 1 x i =0 ẋ i =0;, x i =0 ẋ i =0;, x i =0 ẋ i =0;, x i =0 ẋ i =0;,,,, x i =0 ẋ i =0;, x i =0 ẋ i =0;, x i =0 ẋ i =0;, x i =0 ẋ i =0;, x i =0 ẋ i =0;, ;,,,, (89) The results so obtained by evaluating the gradients of the RKPM based force state maps in the neighborhood of x = 0andẋ = 0areshownbelow: [ ] K est. = /m; C est. = [ ] s/m. It may be noted that in evaluating these gradients the derivative reproduction capabilities of RKPM are used. It was noted that the accuracy of identification is influenced by choice of n in determining the RKPM shape functions. A choice of n = 80 was found to provide satisfactory results. It is observed that the identified values of system parameters compare very well with the corresponding reference values. The nonlinearity parameter obtained using RKPM approximation is α = /m 3, which coincides with the reference value of α. 4.3 Case 4: SDOF systems with non-smooth nonlinearities We exemplify this class of problems by considering the system mẍ +cẋ +kx+k s xu (x 0)+μsgn(ẋ)=g(t); x(0)=0, ẋ(0)=0 (90) This system has a piecewise linear restoring force characteristics resulting from a tension-only spring k s and the Coulomb s damper with coefficient μ. Specifically, the restoring force surface is discontinuous at the origin in both x and ẋ. The starting point in developing models for restoring force is the measured restoring force surface data such as those shown in Figures Based on a careful inspection of these plots, it is possible to anticipate the possible existence of non-smooth characteristics in the restoring force surface or its derivatives. It is apparent that the dilemma of deciding if the restoring surface, or its derivative, have discontinuities would not exist if the form of the nonlinearity is taken to be known a priori and we limit the identification problem to the determination of the assumed model. For purpose of illustration, the system parameters considered are m = 1.77Kg, k = 500/m, k s = 1.5k, μ = 1.5, g(t) =Pcosλt, P = 10, λ = π f with f = 5.86Hz, and the damping parameter c is selected in such a way that, when k s = 0and μ = 0, the system has damping ratioη = 0.0. The response is determined for 0 t f 10s with Δt = 00ω π n s, whereω n = k/m. A visual inspection of this force state map reveals the possible

26 148 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 (c) (d) Figure 17: SDOF system with only bilinear stiffness: k s = 1.5k/m (Case 4a); phase plane plot; measured restoring force surface; (c) projection of the surface along ẋ = 0; (d) projection of the surface along x = 0. discontinuities in the restoring force surface at the origin. Following this finding, the domain Ω : < x,ẋ < could be divided into sub domains Ω i ; i = 1,,... and an enriched RKPM representation (section.3) is used subsequently to construct the refined force state maps. We consider, for the purpose of illustration, three cases: Case 4: k s 0, μ = 0; Ω 1 : x x < 0; ẋ ẋ ẋ Ω :0 x x; ẋ ẋ ẋ Case 4: k s = 0, μ 0; Ω 1 : x x x; ẋ ẋ < 0 Ω : x x x; 0 ẋ ẋ Case 4(c): k s 0, μ 0; Ω 1 :0< x x; 0< ẋ ẋ Ω : x x 0; 0 < ẋ ẋ Ω 3 : x x 0; ẋ ẋ 0 Ω 4 : 0< x x; ẋ ẋ 0 For each of these cases we develop the enriched RKPM based force state map and there results are shown in Figures 0-(c). In fitting the enriched RKPM based force state map, it is assumed that n = 50. The determination of system parameters could now be achieved by examining the restoring force surface. Thus, from Figure 0 it can be observed that the restoring force surface

27 Force State Maps 149 (c) (d) Figure 18: SDOF system with only Coulomb s friction damper: μ = 1.5 (Case 4b); phase plane plot; measured restoring force surface; (c) projection of the surface along ẋ = 0; (d) projection of the surface along x = 0. is continuous and differentiable in ẋ for all values of x. Conversely, the surface is not differentiable at x = 0 although, it is continuous for all values of x and ẋ. This clearly, points towards the bilinear characteristics of the system and a linear damping characteristic in its velocity. This is borne out by results shown in Figures 1(a,b) in which the sections of the restoring force surface along x = 0 (Figure 1a) and along ẋ = 0(Figure 1b) are shown. Based on these Figures the value of damping coefficient c and stiffness coefficients k&k s are determined, and are found to be matched the corresponding reference values very well. Similar results for the case 4 and 4(c) are shown in Figures and 3. Here again, the system parameters identified from the sections of restoring force surface agree very well with the corresponding values. 4.4 Case 5: MDOF systems with non-smooth nonlinearities This class of problems perhaps offers the greatest challenge in nonlinear structural system identification. Here again we discuss the RKPM based strategy within the framework of a specific example. Accordingly, we consider a modified version of example considered in section 4. (see Figure

28 150 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 (c) (d) Figure 19: SDOF system with bilinear stiffness: k s = 1.5k/m and Coulomb s friction damper: μ = 1.5 (Case 4c); phase plane plot; measured restoring force surface; (c) projection of the surface along ẋ = 0; (d) projection of the surface along x = 0. 16) given by m 1 ẍ 1 +c 1 ẋ 1 +c (ẋ 1 ẋ )+k 1 x 1 +k (x 1 x ) +αx μsgn(ẋ 1 )=g 1 (t) m ẍ +c (ẋ ẋ 1 )+k (x x 1 )=g (t) x 1 (0)=0; x (0)=0; ẋ 1 (0)=0; ẋ (0)=0 (91) The system parameters, excitation, and initial conditions are selected as in example considered in section 4. with μ = 5. A value of n = 10 was found to be necessary to obtain satisfactory RKPM representations. To simplify the analysis we assume that the functional form of the restoring force surfaces is known a priori. The phase space is divided in to two regions Ω 1 and Ω which are separated by a line of ẋ 1 = 0. The RKPM representation is accordingly enriched to allow for this discontinuity. The system parameters are again estimated by examining the sections of restoring force surface and its derivatives and following estimates of system parameters are obtained: [ ] K est. = /m; C est. (1, )= c = s/m; C est. (, 1)= c = s/m;

29 Force State Maps 151 Figure 1: Sections of enriched RKPM based force state map (Case 4a); section along x = 0; section along ẋ = 0. duced from the sections of restoring force surface (see Figure 4). (c) Figure 0: Enriched RKPM based force state maps for Case 4; case 4a: k s 0;μ = 0; case 4b: k s = 0;μ 0; (c) case 4c: k s 0;μ 0. C est. (, )=c = s/m; α est. = /m 3 The parameters μ and C (1, 1) =c 1 + c are de- 5 Kriging based models for force state maps Here we re-consider the examples considered in section 4.0 (equations 74-77) and obtain the kriging based force state maps (based on equations 5-69). The numerical values of the system parameters and excitation are retained as before. Steps 1-4 for force state map construction remain as in section 4.0. Based on the knowledge of x(t), ẋ(t) and g(t) and, in conjunction with covariance model (third of the models in equation 57) we determine the unknown model parameters

30 15 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , 008 Figure : Sections of enriched RKPM based force state map (Case 4b); section along x = 0; section along ẋ = 0. Figure 3: Sections of enriched RKPM based force state map (Case 4c); section along x = 0; section along ẋ = 0. { } p β j j=1, σ and {θ i } q by optimizing the maximum likelihood function L (equation 67). The constrained nonlinear optimization here is solved numerically using the genetic algorithm subroutine ga solver available on the Matlab platform. In the numerical work it is assumed that p=10 and attention focused on cases in which q=. The ga solver was run for 100 generations. Once the model parameters are thus determined, we divide the x ẋ space into a set of equi-spaced grid points, and at each of these points the values of force state map are estimated using equation 69. The procedure for investigation of fixed points and limit cycles follows the approach similar to what has been used in section 4.0. Fig- ures 5-7 respectively show the numerical results obtained for the cases of Duffing s oscillator with positive linear stiffness, Duffing s oscillator with negative linear stiffness and for Van der Pol s oscillator. In each of these cases the measured x(t),ẋ(t),ẍ(t) and g(t) are seeded by synthetically generate Gaussian noise (with flat broad band spectra) with zero mean and specified standard deviation. The standard deviations are chosen to be % and 5% of maximum values of the respective signals in the absence of noise. A selection of the numerical results for the three cases is as follows: 1. Duffing s oscillator with positive linear stiffness (equation 75):

31 Force State Maps 153 Figure 4: Studies on -dof system having discontinuous and polynomial-type nonlinearities; section of enriched RKPM based force state map along x = 0. a. Standard deviation of noise (at 5% of deterministic maxima): displacement: m, velocity: m/s, acceleration: m/s, and applied force: b. Estimates { of model } parameters: θ = 51.03, , ˆσ = , and { ˆβ =.4, 10.50, 0.39, 0.37, 0.10, 0.00, 9.78, } Duffing s oscillator with negative linear stiffness (equation 76): a. Standard deviation of noise (at 5% of deterministic maxima): displacement: m, velocity: m/s, acceleration: m/s, and applied force: b. Estimates { of model } parameters: θ = 51.03, , ˆσ = , and ˆβ = { 5.90, 3.75, 0.01, 0.90, 0.10, 0.01, 9.57, 0.00, 0.05, 0.00 } (c) Figure 5: Force state maps using kriging model; Example of Duffing s oscillator with positive k; observed noisy data on system response; force state map based on kriging model; (c) restoring force at ẋ = 0.

32 154 Copyright c 008 Tech Science Press CMES, vol.3, no.3, pp , Van der Pol s oscillator (equation 77): a. Standard deviation of noise (at 5% of deterministic maxima): displacement: m, velocity: m/s, acceleration: m/s, and applied force: b. Estimates of model parameters: θ = { 3.789, }, ˆσ = , and { ˆβ = 3.93, 9.40, , 0., 0.04, } 0.00 (c) Figure 6: Force state maps using kriging model; Example of Duffing s oscillator with negative k; observed noisy data on system response; force state map based on kriging model; (c) restoring force at ẋ = 0; the inset shows the details of the curve near the origin illustrating the influence of noise. Figures 5a-c show the results for the case of Duffing s oscillator with positive stiffness (equation 75). Based on the data shown in figure 5a, the kriging based force state map is obtained in figure 5b and the cross section of this surface along ẋ = 0 (that is, the section Ĝ(x,ẋ = 0)) is shown in figure 5c for noise levels of 0%, % and5%.basedontheresultsshowninfigure5c it may be concluded that the kriging based force state map construction works satisfactorily when measurements are corrupted by reasonably large levels of noise. Again, based on Ĝ(x,ẋ = 0) (figure 5c) it may be inferred that the origin is a stable fixed point and it is verified that it is stable. It may be noted that in the representation of the kriging surface (equation 53) terms up to order three have been retained in the study. Similar results for the case of Duffing s oscillator with negative linear stiffness term (equation 76) are shown in figures 6a-c. Again, as may be inferred from figure 6c, the force state map construction is shown to perform satisfactorily in this case also for noise levels of 0%, % and 5%. In the numerical work it was observed that at higher levels of noise, for some realizations of force state map, there were minor fluctuations in the section Ĝ(x,ẋ = 0) and this is depicted in the inset in figure 5c. Such fluctuations could potentially lead to misleading conclusions about the nature of the nonlinearity but however could be correctly interpreted correctly by repeated fits to the force state maps.

33 Force State Maps 155 (c) (d) Figure 7: Force state maps using kriging model; Example of Van der Pol s oscillator; observed noisy data on system response; force state map based on kriging model; (c) restoring force at ẋ = 0; (d) Limit cycle oscillation constructed using kriging based force state map and the observed data. Results on the models for force state maps for Van der Pol s oscillator (equation 77) are shown in figure 7a-d. The prediction on linear stiffness term is observed to be satisfactory from the results in figure 7c for the three noise levels considered. To examine the nonlinear nature of energy dissipation, the free vibration characteristic of the system was studied using the kriging surface fitted using forced vibration data. This results has been compared with the known limit cycle behavior of the system in figure 7d. The departure of the shape of the limit cycle from a circle clearly point towards the strong nature of nonlinearity present and the results presented also point towards successful performance of kriging based force state maps for strongly nonlinear systems even in the presence of significant measurement noise levels. 6 Closingremarks Force state maps offer powerful solutions to problems of nonlinear system identification. The construction of these maps constitutes a nonparametric method for system identification. The method could also be used to determine the parameters of models when the functional form of the restoring force surface is postulated. This typically involves evaluation of gradients of the force state map and (or) exploration of the force state map in specific regions of phase space. The present study has investigated the performance of force state map technique when the restoring force surface is modeled using two hitherto unex-