Efficient Analysis of Stationary Dynamics of Piecewise-Linear Nonlinear Systems Modeled Using General State-Space Representations

Tien, Meng-Hsuan; Lu, Ming-Fu; D'Souza, Kiran

doi:10.1115/1.4054152

Abstract

In this paper, a new technique is presented for parametrically studying the steady-state dynamics of piecewise-linear nonsmooth oscillators. This new method can be used as an efficient computational tool for analyzing the nonlinear behavior of dynamic systems with piecewise-linear nonlinearity. The new technique modifies and generalizes the bilinear amplitude approximation method, which was created for analyzing proportionally damped structural systems, to more general systems governed by state-space models; thus, the applicability of the method is expanded to many engineering disciplines. The new method utilizes the analytical solutions of the linear subsystems of the nonsmooth oscillators and uses a numerical optimization tool to construct the nonlinear periodic response of the oscillators. The method is validated both numerically and experimentally in this work. The proposed computational framework is demonstrated on a mechanical oscillator with contacting elements and an analog circuit with nonlinear resistance to show its broad applicability.

1 Introduction

Piecewise-linear (PWL) systems are often utilized to model the dynamic behavior of many physical systems with nonsmooth properties. The analysis of the dynamics of these systems is of importance in many fields including health monitoring of engineering structures [1–3], design of electrical circuits [4,5], and gene regulatory networks modeling in biological systems [6,7]. The analysis of the nonlinear behavior of these nonsmooth systems is usually computationally challenging, since the techniques that have been developed for analyzing nonlinear dynamical systems require much more computational cost than traditional linear methods [8,9]. Thus, it is crucial to develop efficient computational tools for predicting the dynamics of nonsmooth systems so that efficient parametric studies can be performed in a reasonable time frame.

The nonsmooth properties of engineering systems are usually triggered by discrete events in their dynamic responses. For instance, the intermittent contact observed in engineering structures due to cracks or delamination [10–13], or interface among components [14] often results in nonsmooth features. Also, the nonsmooth dynamics and mechanics are of great importance for modeling jointed structures [15]. A great amount of research has been carried out on the nonlinear modeling of mechanical joints over the last decade [16–18]. Furthermore, the nonsmooth property also occurs in many electrical systems when current-limit elements such as diodes are utilized to realize particular functions [4]. Since the nonlinearity induced by the discrete event usually dominates the dynamics of these nonsmooth systems, they are often modeled using PWL systems in which the dynamics are assumed to behave linearly in each of the systems' subregimes. Therefore, an effective computational tool that can capture the dynamic features of PWL systems is required for analyzing these nonsmooth systems.

Although PWL systems usually consist of multiple linear subsystems, traditional liner analysis techniques such as modal analysis [19] are not able to capture the event-triggered nonlinearity in these systems since they usually exhibit strong nonlinear features [5,20]. Therefore, numerical integration (NI) is usually required to compute the time response of these nonlinear systems. Common NI techniques include Runge–Kutta methods [21] and Newmark methods [22]. Although these methods are commonly employed in nonlinear system analyses, NI usually requires a fairly small time-step to predict the nonsmooth evolution of nonlinear variables in nonsmooth systems and hence reduces computational efficiency significantly. Recently, a new class of transient dynamic analysis tools has been developed to speed up time marching in the simulation process [23,24]. These techniques are referred to as the hybrid symbolic numeric computational (HSNC) method and have been shown to be more efficient than NI for simulating the dynamic response of PWL systems. The HSNC method employs the analytical solution of the system in each of its linear regimes and stitches these PWL responses together using a numerical incremental search process. Since the method only requires numerical computation at transition points where the system switches its state, the computational cost can be reduced significantly compared to traditional NI methods. The HSNC method is particularly useful for analyzing transient and chaotic responses of complex PWL systems. The method has been recently adapted to analyze the nonlinear rattling behavior of gear systems [25,26].

Although efficient transient dynamic analyses for PWL systems can be achieved using HSNC, a method that can directly capture the steady-state dynamics of these systems is preferred when parametric studies are needed. Thus, another branch of techniques based on the harmonic balance (HB) method has been developed to capture the steady-state vibrational response of PWL nonsmooth systems [27–31]. These methods have been applied to a broad array of nonlinear systems and are not limited to PWL applications. These HB-based methods approximate the periodic response of nonlinear systems using a truncated Fourier series with the unknown coefficients being solved using numerical methods. Although these methods are more efficient than transient dynamic analysis tools, numerous harmonics are usually required to capture the nonsmoothness in PWL systems in particular. Thus, these methods still require considerable computational cost and often incur convergence issues when dealing with PWL systems.

Recently, a new efficient method has been created to approximate the steady-state dynamic response of PWL systems. This method is referred to as the bilinear amplitude approximation (BAA) method [32] and is applicable for capturing the vibrational response of a subset of PWL nonsmooth systems whose dynamic behavior is dominated by two distinct linear states. The BAA method is based on a similar idea as that of HSNC; namely, the nonlinear dynamic response of PWL systems can be obtained by stitching together the analytical responses of the system in each of its linear regimes. The method derives a set of constraint equations to enforce the continuity condition of the PWL responses within a vibrational cycle. A numerical optimization solver is then used to solve for the unknown parameters in the analytical responses when combining each individual linear response to obtain the nonlinear vibrational response. It has been shown that the BAA method is more efficient than the transient dynamic analysis tools and HB-based techniques [32]. Furthermore, the BAA method does not have similar convergence issues as that of HB-based methods, since this method only requires a few linear modes to obtain an accurate approximated response. The BAA method has been successfully applied to investigate the nonlinear dynamics of complex structures with cracks [33]. Also, the method has been extended to develop an efficient technique that is able to capture the vibrational response of systems with Coulomb friction [34]. Note that the idea of utilizing a semi-analytic approach for computing the dynamic response of PWL systems was also used in some previous studies [35,36]. These studies are focused on investigating the nonlinear behaviors of single degree-of-freedom (DOF) oscillators in contrast to the BAA method, which was developed for enabling the analysis of high-dimensional structural systems that contain stationary dynamics.

Although the BAA method is an efficient and effective tool, it is developed based on the standard form of equations of motion of structural systems. It also requires the system to be proportionally damped since its formulation to obtain the analytic solution of the modal response is built off of second-order differential equations. Hence, the goal of this paper is to extend the original BAA method to more general PWL systems so that physical systems in a variety of disciplines (e.g., electrical engineering) modeled using state-space models can also be efficiently investigated. In the extended method, the complex analytical expression of PWL responses of bilinear state-space models is derived, and a numerical tool is used to find the parameters in the analytical expressions by enforcing appropriate compatibility conditions. Nonlinear vibrational responses of the system can then be obtained by coupling the analytical solutions of the system in each linear regime. The extended method is as efficient as the original BAA method since they both utilize efficient analytical approaches in determining the PWL responses. However, the extended method is also able to analyze state-space models whose modal properties are quantified using complex numbers. This capability significantly extends the applicability of the method. Also, the proposed method is validated by both numerical simulation and experimental measurement in this work.

The remainder of the paper is organized as follows. First, the proposed method is introduced. Next, the results of applying the proposed method for a spring-mass mechanical oscillator and a nonlinear electrical circuit are presented. Finally, conclusions are discussed.

2 Methodology

In this section, an extended BAA method is introduced for analyzing general PWL systems represented by a state-space model. In general, the state-space representation of engineering systems can be expressed as

\dot{u} = A (u) u (t) + B (u) q (t)

(1)

where

A (u)

and

B (u)

represent the state matrix and input matrix, respectively. Note that

A (u)

has the size of

n \times n

and

B (u)

has the size of

n \times p

⁠. Furthermore,

u (t) \in {IR}^{n}

and

q (t) \in I R^{p}

represent the state vector and input vector, respectively. For systems whose dynamics can be modeled using two distinct linear subsystems, the state matrix and input matrix can be formulated as

\begin{matrix} A (u) = {\begin{matrix} A_{1} & if u_{nl} \geq a \\ A_{2} & if u_{nl} < a \end{matrix} \\ B (u) = {\begin{matrix} B_{1} & if u_{nl} \geq a \\ B_{2} & if u_{nl} < a \end{matrix} \end{matrix}

(2)

where $u_{nl}$ represents the nonlinear state variable that determines the system's state, $A_{1}$ and $A_{2}$ represent the state matrices in the first and second linear states, and $B_{1}$ and $B_{2}$ represent the input matrices in the first and second linear states. These matrices are all constant matrices since they represent the system matrices of the linear subsystems. Note that a in Eq. (2) represents a threshold that distinguishes the first and second linear states.

Since the dynamic response of the system in each of its linear regimes can be expressed analytically, eigendecompositions

A_{1}

and

A_{2}

are performed to obtain the closed-form solution of the linear subsystems. However,

A_{1}

and

A_{2}

in state-space models are generally nonsymmetric matrices and the eigenvalues of these matrices could be complex. Also, the eigenvectors of these matrices are not mutually orthogonal. This leads to the following adjoint eigenvalue problem. The eigenvalue problem of the matrices A's can be expressed as

A x_{j} = λ_{j} x_{j}, where j = 1, \dots, n

(3)

λ_{j}

and

x_{j}

represent the eigenvalues and right eigenvectors of the state matrices A's. Note that if

λ_{j}

and

x_{j}

are a set of complex eigenvalue and eigenvector, then the complex conjugates

{\bar{λ}}_{j}

and

{\bar{x}}_{j}

are also a set of eigenvalue and eigenvector. Next, the adjoint eigenvalue problem of the transposed matrices

A^{T}

can be formulated as

A^{T} y_{r} = λ_{r} y_{r}, where r = 1, \dots, n

(4)

Note that

A^{T}

and A share the same eigenvalues, and the right eigenvectors

x_{j}

's and the left eigenvectors

y_{r}

's satisfy the following biorthogonality:

y_{r}^{T} A x_{j} = 0, if λ_{j} \neq λ_{r}

(5)

These eigenvectors can be normalized, by enforcing

y_{r}^{T} x_{j} = δ_{j r}

where

δ_{j r}

is the Kronecker delta, to obtain

Y^{T} A X = Λ

(6)

where the matrix

Y = [y_{1}, \dots, y_{n}]

includes the left eigenvectors, the matrix

X = [x_{1}, \dots, x_{n}]

includes the right eigenvectors, and

Λ = diag [λ_{1}, \dots, λ_{n}]

is the corresponding eigenvalue matrix. Equation can be reformulated as

A = X Λ Y^{T}

(7)

Note that Eq. represents the modal decompositions of the matrices A's. Next, the state vectors

u (t)

of the linear subsystems in Eq. can be expressed as a linear combination of the right eigenvectors by applying the state space expansion theorem [19]

u (t) = X ξ (t)

(8)

where

ξ (t)

include the modal coordinates. Then, substituting Eq. into the linear subsystems in Eq. and premultiplying by

Y^{T}

⁠, the equations of motion in the modal space are obtained

\dot{ξ} (t) = Λ ξ (t) + n (t)

(9)

where

n (t) = Y^{T} B q (t)

⁠. Note that Eq. is a set of independent modal equations

{\dot{ξ}}_{j} (t) = λ_{j} ξ_{j} (t) + n_{j} (t), where j = 1, \dots, n

(10)

Note that $n_{j} (t) = y_{j}^{T} B q (t)$ ⁠.

Next, a hybrid analytic-numerical process is employed to obtain the nonlinear forced response of the PWL systems. If the system is excited by a harmonic force

q (t) = q_{0} e^{i ω t}

⁠, the modal coordinates of the system in its linear states can be represented as combinations of the transient dynamic motion and the steady-state dynamic motion

\begin{matrix} ξ_{j, 1} (t) = c_{j, 1} e^{λ_{j, 1} t} + \frac{y_{j, 1}^{T} B_{1} q_{0}}{i ω - λ_{j, 1}} e^{i (ω t + ψ)} \\ ξ_{r, 2} (t) = c_{r, 2} e^{λ_{r, 2} t} + \frac{y_{r, 2}^{T} B_{2} q_{0}}{i ω - λ_{r, 2}} e^{i (ω t + ψ)} \end{matrix}

(11)

where $ξ_{j, 1} (t)$ (⁠ $j = 1, \dots, n$ ⁠) represents the modal response of the PWL system in the first state; $ξ_{r, 2} (t)$ (⁠ $r = 1, \dots, n$ ⁠) represents the modal response of the system in the second state; $c_{j, 1}, c_{r, 2} \in ℂ$ are scalar complex amplitude coefficients of the transient responses; $ω \in I R$ represents the frequency of excitation; and $ψ \in I R$ represents the phase difference between the excitation and the steady-state responses. Note that ψ is included as an additional phase shift so that the time origin can be reset when coupling these PWL responses. The coupling process will be explained shortly. Also note that the phrases “transient response” and “steady-state response” represent the responses of the system in each of its linear states and should not be confused with the total nonlinear response of the nonsmooth system. Furthermore, if ξ_j associated with the eigenvector $x_{j}$ is a complex coordinate, its complex conjugate $\bar{ξ_{j}}$ is also a modal coordinate and the eigenvector associated with $\bar{ξ_{j}}$ is $\bar{x_{j}}$ ⁠, the complex conjugate of $x_{j}$ ⁠. Thus, complex coefficients c_j and ${\bar{c}}_{j}$ exist in conjugate pairs. The physical coordinates of the system in its linear regimes can be obtained by applying the transformation $u_{1} (t) = X_{1} ξ_{1} (t)$ and $u_{2} (t) = X_{2} ξ_{2} (t)$ ⁠, where the subscripts 1 and 2 represent the two states in Eq. (8). The key idea of the method is to couple $u_{1} (t)$ and $u_{2} (t)$ to construct the nonlinear vibrational response of the PWL system by applying appropriate compatibility conditions. Note that the modal expression Eq. (11) used in this work is derived for complete eigenproblems. For systems whose state matrices are defective, the generalized eigenvectors are required to complete the eigenbasis and Eq. (11) needs to be modified for those cases [37].

A steady-state vibrational cycle of the PWL system is schematically shown in Fig. 1. Note that when the nonlinear state variable

u_{n l} \geq a

⁠, the system's dynamics is described by the first linear subsystem; when

u_{n l} < a

⁠, the system's dynamics is governed by the second linear subsystem. T₁ represents the period of time the system stays in the first linear state; T₂ represents the period of time the system stays in the second linear state; and T represents the full period of the nonlinear vibrational cycle. In order to construct this vibrational cycle, a set of compatibility conditions is derived and applied to obtain the unknowns

c_{j, 1}, c_{r, 2}

⁠, and

ψ

in Eq.

\begin{matrix} u_{1} (T_{1}) = u_{2} (T_{1}) \\ u_{1} (0) = u_{2} (T) \\ u_{n l, 1} (0) = a \\ u_{n l, 1} (T_{1}) = a \end{matrix}

(12)

Fig. 1

View large Download slide

One vibrational cycle for the nonlinear state variable u_nl

The first equation and the second equation in Eq. represent the continuity condition for all of the state variables. The variables must satisfy the first condition at the moment the system switches from the first state to the second state; similarly, the variables must satisfy the second condition at the moment the system switches from the second state to the first state. The third and fourth equations represent the threshold value of the nonlinear variable u_nl at the moments of transition. Note that the period of vibration T can be determined by the user. If the system responds at the excitation frequency, then the period of vibration can be calculated by

T = \frac{2 π}{ω}

⁠. If sub- or super-harmonic motions are of interest, then, the period is assumed to be a multiple or a factor of

T = \frac{2 π}{ω}

⁠. Furthermore, T₁ in Eq. is an additional unknown that needs to be obtained to construct the full nonlinear vibrational cycle and T₂ can be calculated by

T - T_{1}

once T₁ is obtained. Next, Eq. is employed in Eq. to obtain the modal representation of the compatibility conditions

\begin{matrix} X_{1} ξ_{1} (T_{1}) - X_{2} ξ_{2} (T_{1}) = 0 \\ X_{1} ξ_{1} (0) - X_{2} ξ_{2} (T) = 0 \\ X_{nl, 1} ξ_{1} (0) - a = 0 \\ X_{nl, 1} ξ_{1} (T_{1}) - a = 0 \end{matrix}

(13)

where $X_{nl, 1}$ is the portion of the mode shapes corresponding to the nonlinear variable $u_{nl, 1}$ ⁠. In other words, $X_{n l, 1}$ represents the kth row of $X_{1}$ ⁠, where the kth variable in $u (t)$ is the nonlinear state variable u_nl(t). The unknowns (⁠ $c_{j, 1}, c_{r, 2}, ψ$ ⁠, T₁) in Eqs. (11) and (13) can be solved using numerical optimization tools. The optimization solver “lsqnonlin” in matlab [38] was used in this work to find these parameters by minimizing the residual in Eq. (13). The solver uses the trust-region-reflective algorithm to iteratively search for a local minimum of the objective function until a tolerance is achieved. Note that if the system is driven by $q_{0} cos ω t$ ⁠, only the real part of Eq. (13) needs to be solved; if the system is excited by $q_{0} sin ω t$ ⁠, only the imaginary part of Eq. (13) needs to be solved. In practice, one can arbitrarily choose either the real or imaginary part to solve since the phase delay in the excitation does not affect the steady-state vibration properties of the system. A steady-state vibrational cycle of the PWL system can then be constructed by coupling the solved linear responses of each linear state.

Since the proposed method is able to obtain the steady-state response of PWL systems, it is particularly useful for conducting parametric studies. A parametric sweep can be efficiently performed by starting from an initial parameter value. At this initial parameter value, multiple random initial guesses for (⁠ $c_{j, 1}, c_{r, 2}, ψ$ ⁠, T₁) are provided to the nonlinear optimization solver and the solution with the minimum residual of Eq. (13) is chosen as the initial values for the second parameter point. In the remaining sweep process, the solution from the previous parameter value is used as initial values for computing the variables at the next parameter point. This sweep strategy ensures convergence while speeding up the analysis. Note that although the method is proposed for PWL systems with two distinct linear subsystems in this work, Eq. (13) can be expanded to include additional compatibility conditions for systems with more linear states (i.e., more general PWL systems that are not simply bilinear). Also note that the stability of the obtained solutions can be determined by the modified Floquet theory proposed by Leine [39].

3 Results

The proposed method is demonstrated on both mechanical and electrical systems to verify its broad capability of analyzing PWL nonsmooth systems in general fields. Results of applying the methodology to a mechanical oscillator with contacting elements and a nonlinear analog circuit with a diode are discussed.

3.1 Mechanical Oscillator With Contacting Masses.

A three DOF spring-mass oscillator with contacting elements as shown in Fig. 2(a) is analyzed using the proposed method. The masses, spring constants, and damping coefficients are set to be

m_{1} = m_{2} = 2.0

kg,

m_{3} = 10

kg,

k_{1} = 1.20

N/m,

k_{2} = 1.68

N/m,

k_{3} = 8.00

N/m,

c_{1} = 0.060

kg/s,

c_{2} = 0.0168

kg/s, and

c_{3} = 0.400

kg/s. Also, the initial gap size is g = 0 m. Note that the damping coefficients are chosen such that the system is not proportionally damped and must be modeled using a state-space representation; thus, the system studied in this work represents a fairly general mechanical oscillator that the original BAA method [32] would not be able to handle. In this case study, a harmonic excitation

F (t) = f_{0} cos (ω t)

with amplitude

f_{0} = 0.01

N is applied on mass m₃. In order to model the intermittent contact phenomenon between the masses m₁ and m₂, two linear subsystems are used to model the dynamics of the nonsmooth oscillator in each of its linear states. The first subsystem shown in Fig. 2(b) represents the linear state of the oscillator when the gap is open; the second subsystem shown in Fig. 2(c) represents the linear state of the oscillator when the gap is closed. In the second subsystem, a set of contact stiffness

k^{*}

and damping

c^{*}

is added to the model for minimizing the penetration between the contacting masses. In this case study, the values

k^{*} = 1000

N/m and

c^{*} = 50

kg/s are used to model the contact stiffness and damping, which are at least two orders of magnitude larger than the structural stiffness and damping. The equations of motion of the nonsmooth oscillator can be expressed as

\begin{matrix} M_{1} {\ddot{x}}_{1} + C_{1} {\dot{x}}_{1} + K_{1} x_{1} = F (t) if x_{1} - x_{2} < 0 \\ M_{2} {\ddot{x}}_{2} + C_{2} {\dot{x}}_{2} + K_{2} x_{2} = F (t) if x_{1} - x_{2} \geq 0 \end{matrix}

(14)

where the subscripts 1 and 2 of the system matrices and the coordinate vectors

x

's represent the linear states of the system when the gap is open and closed, respectively. The coordinate vectors

x = {[x_{1}, x_{2}, x_{3}]}^{T}

represent the displacements of m₁, m₂, and m₃, respectively. The matrices in Eq. can be expressed as

\begin{matrix} M_{1} = M_{2} = [\begin{matrix} m_{1} & 0 & 0 \\ 0 & m_{2} & 0 \\ 0 & 0 & m_{3} \end{matrix}] \\ C_{1} = [\begin{matrix} c_{1} & 0 & - c_{1} \\ 0 & c_{2} & - c_{2} \\ - c_{1} & - c_{2} & c_{1} + c_{2} + c_{3} \end{matrix}] \\ C_{2} = [\begin{matrix} c_{1} + c^{*} & - c^{*} & - c_{1} \\ - c^{*} & c_{2} + c^{*} & - c_{2} \\ - c_{1} & - c_{2} & c_{1} + c_{2} + c_{3} \end{matrix}] \\ K_{1} = [\begin{matrix} k_{1} & 0 & - k_{1} \\ 0 & k_{2} & - k_{2} \\ - k_{1} & - k_{2} & k_{1} + k_{2} + k_{3} \end{matrix}] \\ K_{2} = [\begin{matrix} k_{1} + k^{*} & - k^{*} & - k_{1} \\ - k^{*} & k_{2} + k^{*} & - k_{2} \\ - k_{1} & - k_{2} & k_{1} + k_{2} + k_{3} \end{matrix}] \end{matrix}

(15)

Fig. 2

View large Download slide

(a) Three DOF spring-mass oscillator with contacting elements. (b) The linear subsystem associated with the open state. (c) The linear subsystem associated with the closed state.

Since the linear subsystems are not proportionally damped, the equation of motion Eq. needs to be transformed to its state-space form as shown in Eq. . In the state-space representation the coordinate system is expanded to

u = {[x_{1}, x_{2}, x_{3}, {\dot{x}}_{1}, {\dot{x}}_{2}, {\dot{x}}_{3}]}^{T}

and the corresponding state matrices and the input matrices are

\begin{matrix} A_{1} = [\begin{matrix} 0 & I \\ - M_{1}^{- 1} K_{1} & - M_{1}^{- 1} C_{1} \end{matrix}] \\ A_{2} = [\begin{matrix} 0 & I \\ - M_{2}^{- 1} K_{2} & - M_{2}^{- 1} C_{2} \end{matrix}] \\ B_{1} = [\begin{matrix} 0 \\ M_{1}^{- 1} \end{matrix}] \\ B_{2} = [\begin{matrix} 0 \\ M_{2}^{- 1} \end{matrix}] \end{matrix}

(16)

where I represents the identity matrix.

Next, the nonlinear vibrational response of the system is computed using the proposed method. Since the state of the system depends on the relative displacement of m₁ and m₂, the nonlinear variable for this mechanical oscillator is

u_{n l} = x_{1} - x_{2}

and the threshold value is set to a = 0. In order to apply the proposed algorithm, the original state vector u is converted to an appropriate coordinate system

\bar{u} = {[u_{n l}, x_{2}, x_{3}, {\dot{x}}_{1}, {\dot{x}}_{2}, {\dot{x}}_{3}]}^{T}

by applying the transformation

u = α \bar{u}

⁠, where

α = [\begin{matrix} 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}]

(17)

The system matrices associated with the converted state vector can also be obtained by applying the transformation ${\bar{A}}_{1} = α^{T} A_{1} α, {\bar{A}}_{2} = α^{T} A_{2} α, {\bar{B}}_{1} = α^{T} B_{1} α$ ⁠, and ${\bar{B}}_{2} = α^{T} B_{2} α$ ⁠. Note that u_nl is the only variable that determines the state of the system in the converted coordinate system. The proposed method is then applied to the transformed system, and the computed forced response is compared with the one obtained using NI in Fig. 3. For performing NI, the explicit Runge–Kutta method [21] and the event function in matlab [38] are employed in the computation. Figure 3 shows that the steady-state responses obtained using the proposed method and NI have an excellent agreement over the plotted frequency range in which multiple resonances are accurately captured by the extended BAA method. Furthermore, the proposed method only requires 0.039 s to obtain the steady-state vibrational response for a specific frequency. By contrast, the average CPU time required by NI is 2.557 s. The computations are conducted using a Dell XPS 15 laptop (Round Rock, TX) (2.60 GHz), and the extended BAA method only requires 1.53% of the CPU time of NI. The proposed method is expected to have increasing computational savings with respect to NI in the same way as the original BAA method as the system investigated becomes more and more complex [33,40].

Fig. 3

Steady-state forced response of the three DOF spring-mass oscillator. [(−−),(−),(···)] represents the responses of [m1,m2,m3] computed using the proposed method, and [(□),(◯),(△)] represent the responses of [m1,m2,m3] computed using NI.

View large Download slide

Steady-state forced response of the three DOF spring-mass oscillator. $[(- -), (-), (\cdot \cdot \cdot)]$ represents the responses of $[m_{1}, m_{2}, m_{3}]$ computed using the proposed method, and $[(□), (◯), (△)]$ represent the responses of $[m_{1}, m_{2}, m_{3}]$ computed using NI.

3.2 Nonlinear Analog Circuit With a Nonsmooth Resistance.

In this section, the proposed method is used to analyze the response of the analog circuit shown in Fig. 4. This circuit is composed of a capacitor (C), an inductor (L), several resistors (R, R₁, R₂), and a diode. Note that the resistor in the circuit is PWL since the diode restricts the current flow to one direction. If the circuit system is excited by an AC voltage source with frequency

ω

and amplitude

f_{0}

⁠, the system can be modeled as

\dot{u} = A (v) u (t) + B sin (ω t)

(18)

where

u (t) = {[v (t), i (t)]}^{T}

⁠. Note that the state variable v(t) represents the voltage across the capacitor C, and the state variable i(t) is the current across the inductor L and the resistor R. The system matrices in Eq. can be expressed as

\begin{matrix} A (v) = {\begin{matrix} [\begin{matrix} \frac{- (R_{1} + R_{2})}{C R_{1} R_{2}} & \frac{1}{C} \\ \frac{- 1}{L} & \frac{- R}{L} \end{matrix}] if v \geq 0 \\ [\begin{matrix} \frac{- 1}{C R_{1}} & \frac{1}{C} \\ \frac{- 1}{L} & \frac{- R}{L} \end{matrix}] if v < 0 \end{matrix} \\ B = [\begin{matrix} 0 \\ \frac{f_{0}}{L} \end{matrix}] \end{matrix}

(19)

Fig. 4

View large Download slide

A nonlinear circuit with PWL resistance: (a) circuit diagram and (b) physical realization of the circuit on a breadboard

The steady-state response of the circuit is analyzed using the proposed method. In this study, the parameter values are chosen as $R = 50 Ω$ ⁠, $L$ = 10 mH, $R_{1} = 2500 Ω, R_{2} = 500 Ω$ ⁠, and $f_{0} = 9$ V. The phase portraits of the system for a few capacitor values (⁠ $C = 0.01 μ$ F, $1 μ$ F, and $100 μ$ F) and driving frequencies (⁠ $ω = 10$ Hz, $100$ Hz, $1$ kHz, and $10$ kHz) computed using the proposed technique and NI are compared in Fig. 5. In each plot, $(-)$ and $(◯)$ represent the solutions computed using the proposed technique and NI, respectively. Note that the plots in the same row represent phase portraits with the same driving frequency, and the plots in the same column represent phase portraits with the same capacitor values. Figure 5 shows an excellent agreement between the solutions computed using the proposed method and NI solutions. Furthermore, the system exhibits stronger nonlinearity as the capacitor value and/or the driving frequency is decreased since the trajectories in these phase portraits are more twisted.

Fig. 5

View large Download slide

Phase portraits of the nonlinear circuit for a few capacitor values and driving frequencies. In each plot, $(-)$ and $(◯)$ represent the solutions computed using the proposed method and NI, respectively. The horizontal axis represents the voltage v[V] and the vertical axis represents the current i[A] in each plot.

Next, the extended BAA method is used to conduct a frequency sweep to validate its capability of capturing the dynamic response in a broad parameter range. The peak-to-peak voltage of the system for a couple of capacitor values (⁠ $C = 1 μ$ F and $C = 100 μ$ F) is shown in Figs. 6(a) and 6(b), respectively. The excitation frequency range in these case studies is chosen as $[10 Hz, 10^{4.5} Hz]$ ⁠. The frequency sweep curves computed using the proposed method (–) are compared with the ones obtained by NI (⁠ $°$ ⁠) and experimental measurement (+). It is shown that the responses obtained using the proposed method, NI, and experimental measurement exhibits excellent agreement in the plotted frequency. The average CPU time required by the new method to compute the vibrational response for a specific frequency is 0.08 s. By contrast, NI requires 3.01 s to simulate to a steady-state response at each frequency. The new method only requires 2.7% of the CPU time of NI. Thus, the proposed method can be used as an efficient computational framework for designing circuit systems when current-limiting elements are required. Furthermore, this work provides the first experimental validation of the proposed hybrid analytic-numeric method and shows its excellent agreement for analyzing an electrical system. The original BAA method was only validated on a structural system [41].

Fig. 6

View large Download slide

Peak-to-peak amplitude of voltage obtained by the proposed method (–), NI (⁠ $°$ ⁠), and experimental measurement (+): (a) $C = 1 μ$ F and (b) $C = 100 μ$ F

Next, a multivariable parametric study is conducted using the proposed method to analyze the dynamic properties of the circuit. The driving frequency ω and the capacitor C are chosen as the control parameters in this study. Note that this type of parametric study can be efficiently performed using the proposed method and is generally computationally challenging using traditional NI. The peak-to-peak amplitude of voltage and current for the frequency range $[10 Hz, 10^{4.5} Hz]$ and capacitor range $[10^{- 2} μ F, 10^{2} μ F]$ are computed and plotted in Figs. 7(a) and 7(b), respectively. Figure 7 shows that the circuit exhibits a low-pass behavior since the amplitude of voltage drops quickly after a cutoff frequency. Also, the bandwidth shrinks as the capacitor value increases. Furthermore, it is observed that the amplitude of current rises around the cutoff frequency. Although this nonlinear circuit exhibits similar low-pass behavior as a simple linear resistor-inductor-capacitor circuit, the system is different due to how the diode element affects the electrical load and the nonlinearity must be accounted for when analyzing it. The investigation of these properties is critical since the system exhibits a stronger nonlinear behavior in the lower frequency range as discussed above.

Fig. 7

View large Download slide

(a) Parametric study of voltage using the proposed method and (b) parametric study of current using the proposed method

4 Conclusions

In this paper, an efficient methodology for analyzing the steady-state periodic response of piecewise-linear (PWL) nonsmooth systems in general fields is introduced. The technique extends the bilinear amplitude approximation method, which was developed previously for capturing the dynamic response of proportionally damped structural systems, to more general systems so that PWL oscillators modeled using state-space representations can also be efficiently investigated. The new method employs the eigen-decomposablility of a state-space model in distinct linear regimes to obtain the closed-form solution for each linear subsystem. The steady-state vibrational response of the PWL oscillator is then constructed by coupling the closed-form solutions with appropriate compatibility conditions. An efficient computation of the nonlinear dynamics can be achieved using the proposed hybrid analytical-numeric approach.

The method is successfully demonstrated on a mechanical oscillator with contacting elements and an electrical circuit with a diode to show its broad applicability. Also, the method is validated both numerically and experimentally. The proposed method shows the capability of capturing the nonlinear vibrational response of PWL nonsmooth oscillators for a wide parameter range and is an efficient computational framework for investigating the effects of multiple parameters on the system dynamics.

Acknowledgment

This paper is based on work partially supported by the National Science Foundation (United States) under Grant No. 1902408, program manager Dr. Harry Dankowicz, and the Ministry of Science and Technology (Taiwan, R.O.C) under Grant No. MOST 109-2222-E-007-006-MY3. Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the National Science Foundation and the Ministry of Science and Technology.

Funding Data

Ministry of Science and Technology, Taiwan (Grant No. MOST 109-2222-E-007-006-MY3 Funder ID: 10.13039/501100004663).
National Science Foundation (Grant No. 1902408; Funder ID: 10.13039/100000001).

References

1.

Brownjohn

,

J. M. W.

,

De Stefano

,

A.

,

Xu

,

Y.-L.

,

Wenzel

,

H.

, and

Aktan

,

A. E.

,

2011

, “

Vibration-Based Monitoring of Civil Infrastructure: Challenges and Successes

,”

J. Civ. Struct. Health Monit.

,

1

(

3–4

), pp.

79

–

95

.10.1007/s13349-011-0009-5

2.

Doebling

,

S. W.

,

Farrar

,

C. R.

, and

Prime

,

M. B.

,

1998

, “

A Summary Review of Vibration-Based Damage Identification Methods

,”

Shock Vib. Dig.

,

30

(

2

), pp.

91

–

105

.10.1177/058310249803000201

Google Scholar

Crossref

3.

D'Souza

,

K.

, and

Epureanu

,

B. I.

,

2008

, “

Multiple Augmentations of Nonlinear Systems and Generalized Minimum Rank Perturbations for Damage Detection

,”

J. Sound Vib.

,

316

(

1–5

), pp.

101

–

121

.10.1016/j.jsv.2008.02.018

4.

Bajaj

,

N.

,

Chiu

,

G. T.-C.

, and

Rhoads

,

J. F.

,

2018

, “

A Megahertz-Frequency Tunable Piecewise-Linear Electromechanical Resonator Realized Via Nonlinear Feedback

,”

J. Sound Vib.

,

425

, pp.

257

–

274

.10.1016/j.jsv.2018.02.053

Google Scholar

Crossref

5.

Matsumoto

,

T.

,

Chua

,

L.

, and

Komuro

,

M.

,

1985

, “

The Double Scroll

,”

IEEE Trans. Circuits Syst.

,

32

(

8

), pp.

797

–

818

.10.1109/TCS.1985.1085791

Google Scholar

Crossref

6.

Gebert

,

J.

,

Radde

,

N.

, and

Weber

,

G.-W.

,

2007

, “

Modeling Gene Regulatory Networks With Piecewise Linear Differential Equations

,”

Eur. J. Oper. Res.

,

181

(

3

), pp.

1148

–

1165

.10.1016/j.ejor.2005.11.044

Google Scholar

Crossref

7.

Shahrear

,

P.

, and

Glass

,

L.

, and Del Buono, N.,

2015

, “

Analysis of Piecewise Linear Equations with Bizarre Dynamics

,” Ph.D. thesis,

University of Bari Aldo Moro

,

Italy

.

8.

Allemang

,

R.

,

1980

, “

Investigation of Some Multiple Input/Output Frequency Response Experimental Modal Analysis Techniques

,” Ph.D. thesis,

Mechanical Engineering Department, University of Cincinnati

, Cincinnati, OH.

9.

Ewins

,

D. J.

,

1984

,

Modal Testing: Theory and Practice

,

Research Studies Press

,

Taunton, UK

.

10.

Dimarogonas

,

A. D.

,

1996

, “

Vibration of Cracked Structures: A State of the Art Review

,”

Eng. Fract. Mech.

,

55

(

5

), pp.

831

–

857

.10.1016/0013-7944(94)00175-8

Google Scholar

Crossref

11.

Shiryayev

,

O.

, and

Slater

,

J.

,

2010

, “

Detection of Fatigue Cracks Using Random Decrement Signatures

,”

Struct. Health Monit.

,

9

(

4

), pp.

347

–

360

.10.1177/1475921710361324

Google Scholar

Crossref

12.

Burlayenko

,

V. N.

, and

Sadowski

,

T.

,

2010

, “

Influence of Skin/Core Debonding on Free Vibration Behavior of Foam and Honeycomb Cored Sandwich Plates

,”

Int. J. Non-Linear Mech.

,

45

(

10

), pp.

959

–

968

.10.1016/j.ijnonlinmec.2009.07.002

Google Scholar

Crossref

13.

Hein

,

H.

, and

Feklistova

,

L.

,

2011

, “

Computationally Efficient Delamination Detection in Composite Beams Using Haar Wavelets

,”

Mech. Syst. Signal Process.

,

25

(

6

), pp.

2257

–

2270

.10.1016/j.ymssp.2011.02.003

Google Scholar

Crossref

14.

Bilotta

,

A.

,

Faella

,

C.

,

Martinelli

,

E.

, and

Nigro

,

E.

,

2012

, “

Indirect Identification Method of Bilinear Interface Laws for FRP Bonded on a Concrete Substrate

,”

J. Compos. Const.

,

16

(

2

), pp.

171

–

84

.10.1061/(ASCE)CC.1943-5614.0000253

Google Scholar

Crossref

15.

Brake

,

M.

,

2018

,

The Mechanics of Jointed Structures: Recent Research and Open Challenges for Developing Predictive Models for Structural Dynamics

,

Springer

, Cham, Switzerland.

Google Scholar

Crossref

16.

Porter

,

J. H.

,

Balaji

,

N. N.

, and

Brake

,

M. R. W.

,

2022

, “

A Non-Masing Microslip Rough Contact Modeling Framework for Spatially and Cyclically Varying Normal Pressure

,”

Nonlinear Structures & Systems

, Vol.

1

,

G.

Kerschen

,

M. R.

Brake

, and

L.

Renson

, eds.,

Springer International Publishing

, Cham, Switzerland, pp.

53

–

59

.

Google Scholar

Crossref

17.

Gimpl

,

V.

,

Fantetti

,

A.

,

Klaassen

,

S. W.

,

Schwingshackl

,

C. W.

, and

Rixen

,

D. J.

,

2022

, “

Contact Stiffness of Jointed Interfaces: A Comparison of Dynamic Substructuring Techniques With Frictional Hysteresis Measurements

,”

Mech. Syst. Signal Process.

,

171

, p.

108896

.10.1016/j.ymssp.2022.108896

Google Scholar

Crossref

18.

Firrone

,

C. M.

,

Battiato

,

G.

, and

Epureanu

,

B. I.

,

2018

, “

Modeling the Microslip in the Flange Joint and Its Effect on the Dynamics of a Multistage Bladed Disk Assembly

,”

ASME J. Comput. Nonlinear Dyn.

,

13

(

1

), p.

011011

.10.1115/1.4037796

Google Scholar

Crossref

19.

Meirovitch

,

L.

,

2003

,

Fundamentals of Vibration

,

Mcgraw-Hill

, New York.

20.

Thompson

,

J. M. T.

,

Bokaian

,

A. R.

, and

Ghaffari

,

R.

,

1983

, “

Subharmonic Resonances and Chaotic Motions of a Bilinear Oscillator

,”

IMA J. Appl. Math.

,

31

(

3

), pp.

207

–

234

.10.1093/imamat/31.3.207

Google Scholar

Crossref

21.

Dormand

,

J.

, and

Prince

,

P.

,

1980

, “

A Family of Embedded Runge-Kutta Formulae

,”

J. Comput. Appl. Math.

,

6

(

1

), pp.

19

–

26

.10.1016/0771-050X(80)90013-3

Google Scholar

Crossref

22.

Newmark

,

N. M.

,

1959

, “

A Method of Computation for Structural Dynamics

,”

J. Eng. Mech.

,

85

(

3

), pp.

67

–

94

.10.1061/JMCEA3.0000098

23.

Tien

,

M.-H.

, and

D'Souza

,

K.

,

2019

, “

Analyzing Bilinear Systems Using a New Hybrid Symbolic-Numeric Computational Method

,”

ASME J. Vib. Acoust.

,

141

(

3

), p.

031008

.10.1115/1.4042520

Google Scholar

Crossref

24.

Tien

,

M.-H.

, and

D'Souza

,

K.

,

2019

, “

Transient Dynamic Analysis of Cracked Structures With Multiple Contact Pairs Using Generalized HSNC

,”

Nonlinear Dyn.

,

96

(

2

), pp.

1115

–

1131

.10.1007/s11071-019-04844-7

Google Scholar

Crossref

25.

Donmez

,

A.

, and

Kahraman

,

A.

,

2022

, “

Experimental and Theoretical Investigation of Vibro-Impact Motions of a Gear Pair Subjected to Torque Fluctuations to Define a Rattle Noise Severity Index

,”

ASME J. Vib. Acoust.

,

144

(

4

), p.

041001

.10.1115/1.4053264

Google Scholar

Crossref

26.

Donmez

,

A.

, and

Kahraman

,

A.

,

2022

, “

Characterization of Nonlinear Rattling Behavior of a Gear Pair Through a Validated Torsional Model

,”

ASME J. Comput. Nonlinear Dyn.

,

17

(

4

), p.

041006

.10.1115/1.4053367

Google Scholar

Crossref

27.

Kim

,

T. C.

,

Rook

,

T. E.

, and

Singh

,

R.

,

2005

, “

Super- and Sub-Harmonic Response Calculations for a Torsional System With Clearance Nonlinearity Using the Harmonic Balance Method

,”

J. Sound Vib.

,

281

(

3–5

), pp.

965

–

993

.10.1016/j.jsv.2004.02.039

28.

Lau

,

S. L.

, and

Zhang

,

W.-S.

,

1992

, “

Nonlinear Vibrations of Piecewise-Linear Systems by Incremental Harmonic Balance Method

,”

ASME J. Appl. Mech.

,

59

(

1

), pp.

153

–

160

.10.1115/1.2899421

Google Scholar

Crossref

29.

Poudou

,

O.

,

2007

, “

Modeling and Analysis of the Dynamics of Dry-Friction-Damped Structural Systems

,” Ph.D. thesis,

The University of Michigan

, Ann Arbor, MI.

30.

Saito

,

A.

,

Castanier

,

M. P.

,

Pierre

,

C.

, and

Poudou

,

O.

,

2009

, “

Efficient Nonlinear Vibration Analysis of the Forced Response of Rotating Cracked Blades

,”

ASME J. Comput. Nonlinear Dyn.

,

4

(

1

), p.

011005

.10.1115/1.3007908

Google Scholar

Crossref

31.

Sheng

,

W.

,

Hua

,

L.

,

Yang

,

C.

,

Zhang

,

Y.

, and

Tan

,

X.

,

2018

, “

Nonlinear Vibrations of a Piecewise-Linear Quarter-Car Truck Model by Incremental Harmonic Balance Method

,”

Nonlinear Dyn.

,

92

(

4

)0, pp.

1719

–

1732

.10.1007/s11071-018-4157-6

32.

Tien

,

M.-H.

, and

D'Souza

,

K.

,

2017

, “

A Generalized Bilinear Amplitude and Frequency Approximation for Piecewise-Linear Nonlinear Systems With Gaps or Prestress

,”

Nonlinear Dyn.

,

88

(

4

), pp.

2403

–

2416

.10.1007/s11071-017-3385-5

Google Scholar

Crossref

33.

Tien

,

M.-H.

,

Hu

,

T.

, and

D'Souza

,

K.

,

2019

, “

Statistical Analysis of the Nonlinear Response of Bladed Disks With Mistuning and Cracks

,”

AIAA J.

,

57

(

11

), pp.

4966

–

4977

.10.2514/1.J058190

Google Scholar

Crossref

34.

Altamirano

,

G.

,

Tien

,

M.-H.

, and

D'Souza

,

K.

,

2021

, “

A New Method to Find the Forced Response of Nonlinear Systems With Dry Friction

,”

ASME J. Comput. Nonlinear Dyn.

,

16

(

6

), p.

061002

.10.1115/1.4050686

Google Scholar

Crossref

35.

Natsiavas

,

S.

,

1990

, “

On the Dynamics of Oscillators With Bi-Linear Damping and Stiffness

,”

Int. J. Non-Linear Mech.

,

25

(

5

), pp.

535

–

554

.10.1016/0020-7462(90)90017-4

Google Scholar

Crossref

36.

Natsiavas

,

S.

, and

Gonzalez

,

H.

,

1992

, “

Vibration of Harmonically Excited Oscillators With Asymmetric Constraints

,”

ASME J. Appl. Mech.

,

59

(

2S

), pp.

S284

–

S290

.10.1115/1.2899502

Google Scholar

Crossref

37.

Golub

,

G.

,

Golub

,

P.

,

Van Loan

,

C.

,

Van Loan

,

C.

, and

Van Loan

,

P.

,

1996

,

Matrix Computations

,

Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press

, Baltimore, MD.

38.

MATLAB,

2019

,

Version: R2019b

,

The MathWorks

,

Natick, MA

.

39.

Leine

,

R.

,

2000

, “

Bifurcations in Discontinuous Mechanical Systems of the Fillippov-Type

,” Ph.D. thesis,

Mechanical Engineering

, Technische Universiteit Eindhoven, Eindhoven, North Brabant, The Netherlands.

40.

Tien

,

M.-H.

,

Hu

,

T.

, and

D'Souza

,

K.

,

2018

, “

Generalized Bilinear Amplitude Approximation and X-Xr for Modeling Cyclically Symmetric Structures With Cracks

,”

ASME J. Vib. Acoust.

,

140

(

4

), p.

041012

.10.1115/1.4039296

Google Scholar

Crossref

41.

Noguchi

,

K.

,

Saito

,

A.

,

Tien

,

M.-H.

, and

D'Souza

,

K.

,

2022

, “

Bilinear Systems With Initial Gaps Involving Inelastic Collision: Forced Response Experiments and Simulations

,”

ASME J. Vib. Acoust.

,

144

(

2

), p.

021001

.10.1115/1.4051493

Google Scholar

Crossref

Efficient Analysis of Stationary Dynamics of Piecewise-Linear Nonlinear Systems Modeled Using General State-Space Representations

Abstract

1 Introduction

2 Methodology

3 Results

3.1 Mechanical Oscillator With Contacting Masses.

3.2 Nonlinear Analog Circuit With a Nonsmooth Resistance.

4 Conclusions

Acknowledgment

Funding Data

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Efficient Analysis of Stationary Dynamics of Piecewise-Linear Nonlinear Systems Modeled Using General State-Space Representations

Abstract

1 Introduction

2 Methodology

3 Results

3.1 Mechanical Oscillator With Contacting Masses.

3.2 Nonlinear Analog Circuit With a Nonsmooth Resistance.

4 Conclusions

Acknowledgment

Funding Data

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only