Existing solutions for continuous systems with localized, non-smooth nonlinearities (such as impacts) focus on exact methods for satisfying the nonlinear constitutive equations. Exact methods often require that the non-smooth nonlinearities be expressed as piecewise-linear functions, which results in a series of mapping equations between each linear regime of the nonlinearities. This necessitates exact transition times between each linear regime of the nonlinearities, significantly increasing computational time, and limits the analysis to only considering a small number of nonlinearities. A new method is proposed in which the exact, nonlinear constitutive equations are satisfied by augmenting the system's primary basis functions with a set of non-smooth basis functions. Two consequences are that precise contact times are not needed, enabling greater computational efficiency than exact methods, and localized nonlinearities are not limited to piecewise-linear functions. Since each nonlinearity requires only a few non-smooth basis functions, this method is easily expanded to handle large numbers of nonlinearities throughout the domain. To illustrate the application of this method, a pinned–pinned beam example is presented. Results demonstrate that this method requires significantly fewer basis functions to achieve convergence, compared with linear and exact methods, and that this method is orders of magnitude faster than exact methods.
Nonlinearities are pervasive in the engineering problems that are common today; however, either the suite of analytical tools to efficiently study them often require significant approximations or simplifications or the resulting simulations are very computationally expensive. Systems with localized, non-smooth nonlinearities include joints in large structures , gears with backlash [2,3], steam generator tubes in nuclear power plants , rotating shafts with cracks , rotors and other rotating systems with finite clearances [6,7] and fluid conveying pipes with supports  among many other applications [9–11]. Typically, these nonlinearities are approximated by piecewise-linear models; however, the simplification of the non-smooth nonlinearity as a piecewise-linear model can significantly affect the dynamics of the system. Other types of nonlinearities that are of interest in engineering systems include Hertzian contact and coefficient of restitution models [12–14], cubic springs [8,15], Coulomb friction [16,17] and Iwan friction models [18,19], none of which can be accurately modelled with piecewise-linear constitutive models.
The existing analytical studies of nonlinearities in continuous systems is limited for several reasons. If the nonlinearity is expressed as a discontinuity in the domain, a non-smooth spatial and temporal transformation is shown to be effective in removing singular terms from the equation of motion , but the resulting equations still require numerical solution. When the nonlinearity is limited to being expressed as a piecewise-linear model, a number of modal mapping methods exist that calculate the mode shapes of the system in each regime of the piecewise-linear nonlinearity and then use the orthogonality of the mode shapes to map the displacement of the system across the non-smooth transitions of the nonlinearity [4,21–23]. These methods, though, are often constrained in the order of the system [21,22] or the permissible parameter spaces [4,23]. A recent approach eliminates these constraints by proposing a mapping method based on the L2 norm of the system ; however, in all cases, the limiting factor for computational time is finding the exact time that the system transitions from one regime of the piecewise-linear model to another. As long as the nonlinearity can be specified as a piecewise-linear function, these methods are able to exactly satisfy the nonlinearity's constitutive equations (and are thus referred to as exact methods in what follows). However, non-smooth nonlinearities such as impacts, joints and other typical constraints in engineering systems are often not amenable to being modelled in a piecewise-linear manner. Thus, even though exact methods are able to precisely satisfy the piecewise-linear equations, the nonlinearity itself is an approximation of the actual physical phenomenon of interest. Therefore, a solution method that approximately satisfies a realistic model of the nonlinearity may be more accurate than a solution method that exactly satisfies a simplified model of the nonlinearity . In the literature, however, only numerical studies (such as finite difference, finite-element or other discretized methods) are presented in order to study arbitrary, non-smooth, localized nonlinearities in continuous systems.
The method proposed in the present work satisfies the nonlinear constitutive relationships by using a set of non-smooth basis functions. This approach has several advantages over exact and numerical methods: it is neither restricted in the order of the system nor the permissible parameter space, the nonlinearities are not limited to only piecewise-linear constitutive models, and it is computationally more efficient. Part of the increased efficiency is due to the use of continuous mode shapes (when compared with discretized numerical methods) and to not needing to find the exact time of state changes (when compared with exact methods). As a result, this method is well suited for the predictive design and optimization of structures with strong, non-smooth, localized nonlinearities such as contact or dry friction, particularly when the design variable is associated with the nonlinearity's constitutive model. This method originates in the method of discontinuous basis functions, which is developed for the study of large, discrete (finite-element) systems with dampers or joints [26,27] and has non-smooth basis functions that are discontinuous. This method has proved to be effective in analysing discrete systems with arbitrary localized nonlinearities; however, prior to this research no obvious approach for applying this method to continuous systems existed.
The method of augmentation by non-smooth basis functions is first presented from a theoretical perspective in §2. This method is applied to a pinned–pinned beam in §3, and the method is compared with existing techniques for modelling these types of systems. Parameter studies are used to gain more insight into the nature of the nonlinear phenomena modelled, and the extension of the method to multiple nonlinearities and piecewise-nonlinear models is discussed.
2. The method of augmentation by non-smooth basis function
Given a continuous system with position vector , time t, displacement , and locations of the NNL localized nonlinearities in the domain, the equation of motion for the system is expressed as 2.1where is the equation of motion of the system without the nonlinearities, is the constitutive force due to the jth nonlinearity, and f comprises the additional forces acting on the system. Note that in what follows, indicates a vector-valued function, while the absence of an undertilde indicates a scalar-valued function.
(a) Ordinary basis functions
The first step in the method of augmentation by non-smooth basis functions is to find the mode shapes (referred to as the ordinary basis functions in subsequent sections) of a reference system defined as 2.2with the same boundary conditions as the original system. For most systems, this is a linear equation that can be cast in matrix-operator form ; that is, 2.3The operators , and are the mass, gyroscopic and stiffness operators, respectively, and the comma subscript notation is used to denote differentiation with respect to the subsequent variables. (In non-gyroscopic and undamped systems (), the following state–space operator notation is not strictly necessary, and more traditional modal analysis approaches may be used.) Equation (2.3) is cast as 2.4which has known mode shapes 2.5 2.6 and 2.7with natural frequencies ωj. For systems with a vector-valued displacement function , as is the case for a Timoshenko beam or multiple continuous systems that are coupled together, the matrix of stiffness operators must be nonsingular, and the matrix of gyroscopic operators must be antisymmetric. A general damping matrix that violates the conditions necessary for the orthogonality relationships may be absorbed on the right-hand side of equation (2.1) in this formulation in order to guarantee the properties of the inner products. This is illustrated for the case of modal damping in §3.
Defining the inner product of two vectors and over the domain Γ 2.8with denoting the complex conjugate of , and the system matrices 2.9and 2.10orthogonality yields 2.11with the Kronecker delta function denoted by δmn. A further condition placed on the mode shapes is that they must be normalized with respect to (as is evident in equation (2.11)). Since is skew-symmetric 2.12
In the case where is a nonlinear operator, such as the von Karman plate equations, alternative methods such as the method of quadratic components [29,30] are readily available to find the mode shapes and natural frequencies of the system. The restrictions placed on the mode shapes found by alternative methods are that they are orthonormal with respect to the system's mass matrix as above (note that this is not always automatically achieved ). Additionally, the number N of mode shapes needed is determined either by the highest frequency of interest for the analysis of the original system or by a convergence study.
(b) Non-smooth basis functions
A set of systems are defined that are the same as the reference system but with a linear spring at the location of the jth nonlinearity and homogeneous boundary conditions; that is, 2.13where δ(a−b) is the Dirac delta function. Only a small number of mode shapes from equation (2.13) are needed for this system since the non-smooth basis functions are used to augment a set of basis functions that already spans the frequency range of interest, and the solution can be found using the same procedure as for the reference system in §2a as the additional term is linear. The resulting set of mode shapes 2.14 2.15 and 2.16and corresponding natural frequencies ϖj(k) vary as a function of the spring stiffness k. Above a critical stiffness, the ϖj(k) approach an asymptote as the linear spring acts like a pinned connection for the mode shape (at least for fourth-order systems). Likewise, below the same critical stiffness the ϖj(k) approach a second asymptote as the contribution from the linear spring to the forces acting on the system becomes negligible. To select an appropriate k, a lower bound is given as the equivalent stiffness of the system for a load applied to the same location . That is, a static calculation for the deformation of the system parametrized in terms of a prescribed displacement d due to a load P can be expressed as P=kd, where k is the lower bound for the spring stiffness. To determine the non-smooth basis function, the derivative 2.17is needed. Note that the terms ϖj,kφj are neglected as they are not necessary for deriving the non-smooth basis functions. The purpose of taking the derivative of the mode shapes in equation (2.17) is to remove the dependence on k that the mode shapes have. In the instance where this derivative is numerically determined about some nominal value of k, the above guidance can be used to select an appropriate, non-arbitrary k. Because the non-smooth basis functions are derived independent of the discrete nonlinearity's constitutive model, the resulting basis functions do not need to be recalculated even when the nonlinearity is changed. The set of functions must be orthonormalized with respect to the ordinary basis functions, which can be efficiently accomplished using a Gram–Schmidt procedure resulting in the set of non-smooth basis functions 2.18The projection of vector on to vector is defined as 2.19After each non-smooth basis function is found, is normalized via equation (2.11) and the associated natural frequency is given as 2.20This complete set of basis functions spans the entire response space. While each is at least continuous, the derivative (third in fourth-order problems, first in second-order problems) contains a discontinuity at , which allows for the nonlinearity to be accounted for. With the response of the system 2.21the solution of equation (2.1) is readily available, where NNS is the number of non-smooth basis functions used in the expansion solution. The advantage of the present approach is that the exact time of transition from one regime of the nonlinearity to another is not needed as there is no change in basis functions. Consequently, time-stepping methods [31–33] and other temporal discretization methods  may be used. In what follows, an implicit-explicit Runge–Kutta backward-Euler method  is used. While analytical expressions are developed for most quantities (i.e. the mode shapes, natural frequencies, evaluations of the integrals in the force equations, etc.), they are not necessary for analysis by this method; numerical solutions can suffice instead of analytical expressions. For instance, this is the case when modes are coupled by a general damping matrix. Larger systems must be solved, but the model reduction remains efficient.
3. Pinned–pinned beam
This example demonstrates the application of the method of augmentation by non-smooth basis functions to a beam with pinned boundary conditions (figure 1), having length L, position x along the beam, total displacement w(x,t), boundary excitation at x=0, and a piecewise-linear nonlinearity that contains a deadband region located at x=ℓ. This example is chosen since an exact solution is available from  for comparison. Because this example deals with a non-gyroscopic system, modal operator notation is not necessary, and the standard, non-operator based, modal analysis is applied instead.
The equation of motion in terms of w(x,t) is 3.1with the nonlinearity, which represents a stiff flange that bounds the beam's displacement, specified by 3.2and boundary conditions 3.3In the first part of what follows, a symmetric constraint with k(−1)=k(1)=106 N m−1 is used.
(a) Ordinary basis functions
The linear reference system for this problem is defined as 3.4having boundary conditions described by equation (3.3). In order to satisfy the inhomogeneous boundary condition, the superposition function 3.5is used. Defining the elastic displacement u(x,t)=w(x,t)−y(x,t), the equation of motion is rewritten as 3.6Assuming the separable solution 3.7the normalized mode shapes that constitute the ordinary basis functions are 3.8with wavelengths 3.9and natural frequencies 3.10Substitution of (3.7) into (3.6), multiplying by ϕn(x) and integrating with respect to x yields the modal equation of motion 3.11The modal forces due to the superposition are 3.12
(b) Non-smooth basis functions
The non-smooth basis functions are determined by first considering the system 3.13subject to pinned boundary conditions at x=0 and L, with defined in equation (3.4). In order to find the mode shapes of this system, u(x,t), which is defined over x∈[0,L], is conceptually divided into two contiguous regions (denoted by m=1, 2) about x=ℓ, and the problem is recast as 3.14subject to 3.15with u1(x,t) defined over x∈[0,ℓ] and u2(x,t) over x∈[ℓ,L]. Solutions of equations (3.14) and (3.15) are found using the separable solution 3.16and 3.17where λn is found through application of the boundary conditions and 3.18An appropriate stiffness k for calculating the non-smooth basis functions is estimated from 3.19in which the deflection d of the beam subject to a static load P at the location of the nonlinearity is calculated. Rearranging, 3.20
This value of k is verified by calculating the natural frequency ϖn as a function of k, shown in figure 2. Values of k that are too low or too high will yield non-smooth basis functions that cannot adequately account for the nonlinearity in the system (either by not being able to adequately resolve the nonlinearity's restoring forces for low values of k or by not having any displacement at the location of the nonlinearity for high values of k). From equation (3.20), k=0.03 N m−1, which is in the transitionary region for ϖ in figure 2. The sensitivity to k is shown in figure 3 for the contact force and contact duration of the first impact. In order to emphasize the sensitivity, a smaller set of modes than is needed for convergence is used. For the piecewise-linear constraint, convergence is achieved for values of k≥0.02 N m−1 and k≤20 N m−1. As the constitutive model for the discrete nonlinearity is changed, such as to the piecewise-nonlinear nonlinearity of §3d, the converged region is found to extend over 10−3≤k≤8 N m−1. Thus, the estimate for an appropriate k from equation (3.20) is valid for both models. Note that the magnitude of the contact impulse and contact duration for the piecewise-nonlinear nonlinearity are much higher than for the piecewise-linear nonlinearity due to the piecewise-nonlinear nonlinearity being significantly more compliant. The derivative of the mode shapes φn,k(x) are numerically evaluated at the appropriate k, then orthogonalized with respect to the ordinary basis functions via a Gram–Schmidt procedure (2.18) to yield a set of non-smooth basis functions that are non-normalized. The final set of non-smooth basis functions are found by normalizing the orthogonal mode shapes from the Gram–Schmidt procedure with respect to the linear reference system's operators 3.21In a similar manner, the corresponding natural frequency for each non-smooth basis function is found by 3.22Both and φ1 are shown in figure 4. Note that while the mode shape (figure 4a) is continuous, there is a discontinuity in the third derivative (figure 4b). Substitution of the normalized mode shapes into (3.1) and integrating with respect to the normalized mode shapes yields the modal equations of motion for the ordinary and non-smooth modes, respectively 3.23and 3.24with n=1, 2, …, N and j=1, 2, …, NNS. Light, modal damping with damping coefficient ζ is introduced into the system via the procedure of Meirovitch & Ryland .
(c) Nonlinear model comparison
For a system such as that shown in figure 1, it is unclear what metric should be used to test for convergence or model comparison. Because the analysis of this system is focused on the response at the location of the nonlinearity, a plausible metric is the displacement or velocity at the nonlinearity's location; however, for this system, the maximum displacement is bounded and the maximum velocity is approximately bounded. Studying a quantity such as the steady-state impact velocities across a range of frequencies is an effective method to demonstrate convergence , but it can be computationally intensive. A criteria of interest when studying systems with non-smooth nonlinearities is whether the model is able to accurately capture the existing nonlinear resonance. Several metrics exist to classify nonlinear systems, such as Lyapunov exponents and Floquet multipliers, that are useful for determining the stability of the system; however, like other metrics, such as the total strain energy in the system, that may provide a good comparison between two different methods from a system-level perspective, the assessment of how well two different methods account for the nonlinearity is lost.
In studying the convergence of the present method, a feasible approach is to compare the contact time and contact force across the first impact of the system. This metric, thus, assesses the nonlinear transient dynamics of the system rather than the steady-state response. Subsequent analysis should focus on the presence of nonlinear resonances and the ability of the model to adequately capture them; the present work, though, is focused on investigating the convergence characteristics. Figure 5a shows the convergence of the contact impulse, which is calculated by integrating the contact force over the entire impact event, as the number of ordinary basis functions is varied from N=1–40 while NNS is varied from 1 to 3. Above N=24 modes, the contact impulse does not change significantly, and above N=38, the contact duration is approximately constant. As the number of non-smooth basis functions is increased, it is apparent that NNS≥2 is needed for adequate convergence. In figure 6, the convergence of the system with respect to NNS is shown for N=40. Above NNS=4, no change is observed in either the contact impulse or the contact duration.
For comparing the exact solution method of Brake  with the present method, the material and geometric properties for the system are chosen for consistency and are given in table 1. Since the system of figure 1 contains a non-smooth nonlinearity, and is sensitive to even extremely small amounts of numerical error, it is reasonable to ask if any similarities between two different methods is expected. The two systems should be identical prior to the engagement of the nonlinearity; after that, small differences can accumulate over time leading to results that are 180° out of phase with one another. Figure 7 shows the time histories for both  (figure 7a) and the present method (figure 7b). For the present method, N=40 with NNS=4, and for the exact method, N=40 in each piecewise-linear state of the nonlinearity. Similar results are observed for the first 0.5 s; however, as with many systems containing non-smooth nonlinearities, small differences become significant over time. Consequently, time histories are poor metrics for the comparison of two solution methods for the same nonlinear problem, and comparisons are made by studying the contact impulse and duration of the first impact. It is worth noting that the exact method predicts a response in which the beam engages both the upper and lower sides of the constraint three times over five periods of the excitation while the present method predicts a response in which the beam engages both sides of the constraint five times over five periods of the excitation. A reasonable question to ask is: what are the convergence requirements for the exact method?
The advantage of augmenting the set of ordinary basis functions with a set of non-smooth basis functions is shown in figure 8, which compares the convergence of the present method, with NNS=3, and the convergence of the exact method . The set of only ordinary basis functions  requires significantly more basis functions than when the ordinary basis functions are augmented with a set of non-smooth basis functions. When augmented with a set of NNS=3 non-smooth basis functions, approximately N=40 ordinary basis functions are needed. Without augmentation by a set of non-smooth basis functions, even using N=200 ordinary basis functions is insufficient for convergence. In strictly comparing the two methods, it is not surprising that different responses are calculated when the number of basis functions is restricted to N=40. Owing to the computational intensity of calculating the response for a system with several hundred basis functions, it is not practical to find and use a converged set of basis functions for the exact method with this system.
One advantage of using this solution method instead of an exact method is the significant reduction in computation time. In comparing the two methods and using the same N, the present method is 92 times faster than  on average for the simulations used to generate figures 7 and 8 when the same time step is used away from impacts. When converged solutions are used for both methods, this factor increases by orders of magnitude. Another advantage of the present method is that the nonlinearity is not restricted to piecewise-linear functions, as detailed in what follows.
(d) Example of a piecewise-nonlinear model
For the piecewise-linear nonlinearity (equation (3.2)), the displacement and phase diagram of the system are shown in figure 9. The high rebound velocities and no apparent signs of steady-state behaviour seen in figure 9b are indicative of the low level of damping in the system resulting in the displacement of the beam away from the nonlinearity increasing in magnitude with every period of the excitation. In figure 9c, the nonlinear forces acting on the beam as a function of time t is shown. Thick bands of response, such as between 0 and 0.5 s, indicate a chatter response in which the beam is engaging and disengaging the nonlinearity multiple times over a short period of time.
The piecewise-linear nonlinearity of the previous section is replaced with a piecewise-nonlinear model that is defined with the force deflection profile described by 3.25in place of equation (3.2). The restoring force for negative displacements is sufficiently stiff to keep ; however, because the restoring force for positive displacements is significantly less stiff, displacements of w(ℓ,t)>dx are observed (figure 10). Even though the nonlinearity is changed from a piecewise-linear to a piecewise-nonlinear function, no new ordinary basis functions or non-smooth basis functions need to be calculated; the same basis functions as derived for the piecewise-linear nonlinearity are used for the piecewise-nonlinear model. The convergence of the system (with respect to the first impact) for the piecewise-nonlinear model is shown in figure 11. Compared with the piecewise-linear nonlinearity, fewer non-smooth basis functions are needed for convergence with respect to the first impact; however, the use of only one non-smooth basis function yields results that are significantly different from that of models that include higher numbers of non-smooth basis functions. Because the piecewise-nonlinear model includes a restoring force for negative displacements similar to the piecewise-linear nonlinearity, the true convergence of this system (as opposed to convergence across the first impact only) is the maximum number of modes needed to demonstrate convergence for the first impact against each side of the constraint.
The restoring force of the nonlinearity as a function of time is shown in figure 12. Because the restoring force for negative displacements is much stiffer than for positive displacements, the forces must be plotted on separate scales for clarity. Despite a relatively coarse time step, the restoring force for negative displacements is still adequately modelled. This is significant because the main constraint on computational time for this method is time-step size. Since this method does not map from one state of the nonlinearity to another (such as in the exact methods of the earlier studies [4,21,24]), the exact time of contact between the beam and the nonlinearity is not required.
(e) Extension to multiple nonlinearities
In extending the method of augmentation by non-smooth basis functions to multiple nonlinearities, only a small number of additional non-smooth basis functions (in some cases as few as one) are needed, corresponding to each additional nonlinearity. These new non-smooth basis functions must be normalized with respect to the other basis functions; beyond this, though, no additional calculations are needed. Figure 13 shows the convergence for a system with two nonlinearities located at ℓ1=0.33 m and ℓ2=0.67 m. Both nonlinearities are modelled with the constitutive relationship described in equation (3.2), with all other properties given in table 1. The convergence study shows that only two non-smooth basis functions per nonlinearity are needed for convergence. Similar to the systems with only one nonlinearity, approximately 40 ordinary basis functions are needed in order to demonstrate convergence of the contact forces.
As the number of nonlinearities is further increased to four (with ℓ1=0.2 m, ℓ2=0.4 m, ℓ3=0.6 m and ℓ4=0.8 m, and all other parameters are given in table 1), the convergence study (figure 14) shows that only one non-smooth basis function is needed per nonlinearity, which implies that as the number of nonlinearities in this system increases, the total number of basis functions needed for convergence will increase only by the number of additional nonlinearities. Unlike the system with only one nonlinearity (figure 5), the contact duration converges before the contact impulse for this system with four nonlinearities. Figure 15 shows the phase diagram (figure 15a) and the displacement of the beam as a function of x at 50 evenly spaced time increments, each 0.15 s apart (figure 15b). The additional nonlinearities significantly attenuate the magnitude of the displacement along the entire length of the beam, not just at the location of the nonlinearities. Additionally, every time a nonlinearity is engaged, a sharp change in the velocity is seen along the beam. This is particularly evident in the phase diagram for the displacement and velocity of the beam at x=ℓ1 (figure 15a). As a result, high-frequency motion is observed as waves travelling along the beam due to contact with the nonlinearities.
A new method for analysing localized nonlinearities in continuous systems is developed. This method solves for the state of the system by using non-smooth basis functions to account for the nonlinearities. These non-smooth basis functions are used to augment the ordinary basis functions (eigenmodes) of the system. Unlike exact methods that are limited to piecewise-linear nonlinearities, this new method of augmentation by non-smooth basis functions neither divides the system into multiple regimes based on the nonlinearities, nor requires the exact time of engagement of a nonlinearity. Additionally, the only assumption placed on the nonlinearities is that they can be modelled in a discrete manner. An example of a pinned–pinned beam with rigid and compliant amplitude constraints is provided to demonstrate the application of this method. This system is chosen because of the availability of an exact solution  for comparison. Both piecewise-linear and piecewise-nonlinear models are used, and the number of nonlinearities used in the illustrative results are varied. The primary conclusions of this work are
— the method of augmentation by non-smooth basis functions is a robust and efficient technique for modelling the effects of nonlinearities on continuous systems. The advantage of this method lies in not needing separate sets of basis functions for each regime of the nonlinearities, not needing to find the precise times that different regimes of the nonlinearities are engaged and being able to consider any arbitrary nonlinearity;
— changing the constitutive properties of the nonlinearity, such as from a linear spring to a piecewise-nonlinear model, requires no additional calculations as this method is derived independent of the nonlinearity's model;
— only a small number of non-smooth basis functions are needed for convergence. Extension of the method to multiple nonlinearities requires only several additional non-smooth basis functions for each additional nonlinearity. For large numbers of nonlinearities, convergence is achieved with only one non-smooth basis function per nonlinearity. This is in contrast to a mapping method, which requires a geometric increase in the number of basis functions for each additional nonlinearity; and
— when compared with an exact method, the total number of modes needed for convergence is significantly less than that using the present method. In studying the same problem using the present method and an exact method, the exact method required hundreds of additional basis functions for convergence.
We would like to thank our colleagues Michael J. Starr, Jamey T. Bond and Kurtis Ford for their feedback and support.
↵† Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporations, for the US Department of Energy's National Nuclear Security Administration under Contract DE-AC04-94AL85000.
- Received April 26, 2013.
- Accepted July 17, 2013.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.