Abstract

The application of servocontrolled mechanical-bearing-based precision motion stages (MBMS) is well-established in advanced manufacturing, semiconductor industries, and metrological applications. Nevertheless, the performance of the motion stage is plagued by self-excited friction-induced vibrations. Recently, a passive mechanical friction isolator (FI) has been introduced to reduce the adverse impact of friction in MBMS, and accordingly, the dynamics of MBMS with FI were analyzed in the previous works. However, in the previous works, the nonlinear dynamics components of FI were not considered for the dynamical analysis of MBMS. This work presents a comprehensive, thorough analysis of an MBMS with a nonlinear FI. A servocontrolled MBMS with a nonlinear FI is modeled as a two DOF spring-mass-damper lumped parameter system. The linear stability analysis in the parametric space of reference velocity signal and differential gain reveals that including nonlinearity in FI significantly increases the local stability of the system's fixed-points. This further allows the implementation of larger differential gains in the servocontrolled motion stage. Furthermore, we perform a nonlinear analysis of the system and observe the existence of sub and supercritical Hopf bifurcation with or without any nonlinearity in the friction isolator. However, the region of sub and supercritical Hopf bifurcation on stability curves depends on the nonlinearity in FI. These observations are further verified by a detailed numerical bifurcation, which reveals the existence of nonlinear attractors in the system.

1 Introduction

Motion stages provide high-speed and high-precision positioning in manufacturing and metrology-related processes, which include machining, additive manufacturing, and semiconductor fabrications [14]. Among magnetic-based, flexural-based, fluidic-based, and mechanical-bearing-based motion stages (MBMS), MBMS are more prominent in the industries due to their cost-effectiveness, high off-axis stiffness, wide ranges of motions, and easy-to-install feature [5]. The motion control of MBMS in different applications is commonly realized by using different servocontrollers [68]. However, it can lead to the problem of self-excited limit cycles in the tracking error caused by the sliding friction between contact surfaces, which is also known as friction-induced vibrations (FIV) [6,914]. The adverse effects of friction on the control performance include oscillations of stick-slip phenomena tracking errors and long settling times. Therefore, understanding the dynamics of self-excited friction-induced vibrations under different conditions is essential to mitigate tracking error oscillation, which leads to better motion stage performances.

Different controllers have been proposed to eliminate friction-induced vibrations. The fundamental concept of these controllers to suppress FIV is to counteract the friction force by providing an equal and opposite force and suppress these FIV. These controllers can be divided into three categories as (1) high-gain controllers, (2) model-based controllers, and (3) advanced controllers. The use of these controllers in compensating for the effect of friction is well-established in the literature [11,13,1518]. However, the performance of these controllers can be limited by environmental noise in the case of high-gain controllers, model inaccuracy in the case of model-based controllers, and low-performance computational/actuator hardware in the case of advanced controllers.

Recently, a mechanical device known as the friction isolator (FI) was proposed as a more robust and efficient approach to mitigating self-excited friction-induced vibrations in the MBMS system [19,20]. Unlike conventional motion stages in which the bearings are rigidly installed to the table, the compliant motion stages adopt FI as motion-compliant joints between the bearings and the table, isolating the table from nonlinear frictional effects. For more design details of FI, readers are referred to Refs. [19,20]. Furthermore, parametric analyses on the motion stage with FI showed that with optimum design parameter selections, FI could increase the fixed-point stability region to allow the use of higher control gains and reduce the amplitudes of the control error limit cycles [21,22]. However, the above-mentioned studies analyzed the dynamics of an MBMS with linear FI only. Nevertheless, the mechanical design of FI can introduce nonlinear elements, which are significantly compared to the linear elements [19] and hence, need to be considered in the dynamics of compliant motion stages for optimum selection of control parameters. The present paper examines this problem and provides a framework for further analysis. For this, we extend the model proposed in Refs. [22] by introducing a nonlinear element in the system to take into account nonlinear FI. The preliminary linear stability analysis reveals that linear FI underestimates the stability regime, and hence, it is desirable to consider nonlinearity in FI for a better selection of control parameters for steady operations. Furthermore, the nonlinear analysis using the analytical method, more specifically the method of multiple scales, reveals the existence of super and subcritical Hopf bifurcation in the system regardless of the presence of nonlinearity in FI. However, for the given value of system parameters, the region corresponding to supercritical bifurcation can be increased by properly selecting the nonlinear parameters.

The paper is arranged as follows. We outline the extended model of MBMS with nonlinear FI in Sec. 2. Later on, the linear stability analysis of fixed-points is investigated in Sec. 3. This is followed by a detailed analysis of the current system in Sec. 4. Results and detailed discussions from the previous sections are presented in Sec. 5. Finally, some conclusions are drawn in Sec. 6.

2 Mathematical Formulation

To establish the effect of nonlinearity in the friction isolator and the contribution of various system parameters to linear and nonlinear instability of PD-controlled MBMS, a nonlinear two DOF model, illustrated in Fig. 1, is obtained by extending the linear model discussed in Ref. [22]. In this model, the precision motion stage is modeled as a rigid mass mt, whereas the combined mass of the friction isolator and bearing is modeled as mb. The nonlinear interactionsbetween mb and mt are represented by a nonlinear spring with a stiffness function K(.),andalineardamperwithdampingcoefficientcfi. Also, in this model, u1 and r(t) represent the output feedback control force from the PD controller and input reference/set-point signal to the PD controller, respectively. Therefore, if X1(t) and X2(t) represent the position of MBMS and FI, respectively, then the equations of motion for the system are given by

Fig. 1
Schematic of the MBMS with nonlinear friction isolator
Fig. 1
Schematic of the MBMS with nonlinear friction isolator
Close modal
mtX¨1+K(X1X2)+cfi(X˙1X˙2)=u1
(1a)
mbX¨2K(X1X2)+cfi(X˙2X˙1)=Ff
(1b)
In the above governing equations, the feedback force u1 is defined as
u1=kp*βkd*β˙1
(2)
where kp* is the proportional gain, kd* is the differential gain, β is the tracking error defined as the difference between the position and reference input, i.e., β=X1r and, Ff represents the friction force between the bearing and the support platform. For the current analysis, we use the LuGre friction model to represent the friction force between the surfaces as it incorporates hysteresis effects, premotion friction, and viscous friction altogether. Furthermore, the frictional force in the LuGre friction model relates to the relative velocity between the surfaces and on an internal state variable z, which can be interpreted as the average transverse deflection of microscopic asperities/bristles of the contact surfaces [2326]. The friction force in the LuGre friction model is described by the following equation [23]:
Ff=σ0*z+σ1*z˙+σ2*Vr
(3)
whereas z, the average asperity/bristle deflection, is determined by the following equation [23,27]:
z˙=Vrσ0*|Vr|G(Vr)z=Vr(1σ0*sgn(Vr)G(Vr)z)
(4)
In the above equations, Vr is the relative velocity between the contact surfaces, σ0* is the asperity/bristle's stiffness, σ1* is the damping coefficient of the asperity/bristle, whereas σ2* is the friction due to viscosity between the contact surfaces, and G(Vr)>0 is a positive valued function to define the Stribeck effect. In the current work, this positive valued function of relative velocity is chosen in the form of the following equation [27,28]:
G(Vr)=fC*+(fS*fC*)ea˜|Vr|
(5)
where a˜ is the slope parameter. Now with the definition of friction force, Eq. (1) can be rewritten in terms of e as
mtβ¨+kp*β+kd*β˙+g(ββb)+cfi(β˙β˙b)=mtv¨rv
(6a)
mbβ¨bg(ββb)+cfi(β˙bβ˙)=(σ0*z+σ1*z˙+σ2*Vr)mbv¨rv
(6b)
In the above-modified equations of motion, βb is defined as βb=X2r. Since earlier studies suggest that the nonlinearities involved with the LuGre friction model are primarily the combination of quadratic and cubic terms, we assume that the stiffness function involved with the friction isolator has a similar nonlinear restoring force characteristic as our primary system, i.e., the combination of quadratic and cubic terms [29,30]. Therefore, provided that kfl*,kfq*, and kfc* represent the linear, quadratic, and cubic stiffness of FI, respectively, Eq. (6) can be written as
mtβ¨+kp*β+kd*β˙+kfl*(ββb)+kfq*(ββb)2+kfc*(ββb)3+cfi(β˙β˙b)=mtβ¨
(7a)
mbβ¨b+kfl*(βbβ)kfq*(βbβ)2+kfc*(βbβ)3+cfi(β˙bβ˙)=(σ0*z+σ1*z˙+σ2*Vr)mbβ¨b
(7b)
Next, we define the following scales and nondimensional parameters in the system.
x=βX0,xb=βbX0,z˜=zX0,X0=gωp2,ωp=kp*mt,τ=ωpt,ζ=kd*2mtωp,vr=VrX0ωp,σ0=σ0*mtωp2,σ1=σ1*mtωp,σ2=σ2*mtωp,fc=fc*mtX0ωp2,fs=fs*mtX0ωp2,a=a˜ωpX0,κ=cfi2mtωp,kr=kfikp,mr=mtmb,krq=kfq*X0kp*,krc=kfc*X02kp*
(8)
In the above nondimensional scales, g is the acceleration due to gravity. Using the above-mentioned scales and nondimensional parameters, the governing equations of motion (Eqs. (7), (4), and (5)) can be nondimensionalized for the case of constant reference velocity, i.e., v¨rv=0, to get
x¨+2ζx˙+x+kr(xxb)+krq(xxb)2+krc(xxb)3+2κ(x˙x˙b)=0
(9a)
x¨b+krmr(xbx)krqmr(xbx)2+krcmr(xbx)3+2κmr(x˙bx˙)
(9b)
(σ0z˜+σ1vr(1σ0sgn(vr)G(vr)z˜)+σ2vr)z˜˙=vr(1σ0sgn(vr)G(vr)z˜)
(9c)
For ease of analysis of the current system, Eq. (9) can be rewritten compactly in the state space form as
x˙1=x2
(10a)
x˙2=2ζx2x1kr(x1x3)krq(x1x3)2krc(x1x3)32κ(x2x4)
(10b)
x˙3=x4
(10c)
x˙4=2κmr(x4x2)krmr(x3x1)+krqmr(x3x1)2krcmr(x3x1)3mr(σ0x5+σ1vr(1σ0x5G(vr)sgn(vr))+σ2vr)
(10d)
x˙5=vr(1σ0x5G(vr)sgn(vr))
(10e)
with [x1,x2,x3,x4,x5]=[x(τ),x˙(τ),xb(τ),x˙b(τ),z˜(τ)]. Next, if vrv represents the nondimensional constant reference velocity, vr can be rewritten as vr=x˙b+vrv=x4+vrv. For the analytical treatment of the current system, we follow the procedure outlined in Refs. [22,27] and, accordingly, modify Eq. (10) for the case of pure slipping motion sgn(vr)=1 as
x˙1=x2
(11a)
x˙2=krcx13x12krq+3x12krcx3x1krx1+2x1krqx33x1krcx322x2κ2x2ζkrqx32+2κx4+krx3+krcx33
(11b)
x˙3=x4
(11c)
x˙4=krcmrx13+x12krqmr3x12krcmrx3+x1krmr2x1krqmrx3+3x1krcmrx32+mrσ1σ0x5g3x44+x43mrσ1vrvσ0x5g3+x43mrσ1σ0x5g2+x42mrσ1vrvσ0x5g2+x42mrσ1σ0x5g12x4κmr+x4mrσ1σ0x5g0x4mrσ2+x4mrσ1vrvσ0x5g1+mrσ1vrvσ0x5g0mrσ2vrv+krqmrx32+2κmrx2mrσ1vrvkrmrx3krcmrx33mrσ0x5x4mrσ1
(11d)
x˙5=σ0x5g3x44x43σ0x5g2x43vrvσ0x5g3x42vrvσ0x5g2x42σ0x5g1+x4x4vrvσ0x5g1x4σ0x5g0+vrvvrvσ0x5g0.s
(11e)
where gi for i=0,1,2,3 are defined in Appendix  A. Before proceeding further, we need to determine the fixed-points of Eq. (11), and are given by
x1s=g0σ2vrv+1g0,x2s=0,x4s=0,x5s=1σ0g0
(12)
whereas x3s is governed by the roots of the following cubic equation
Ax3s3+Bx3s2+Cx3s+D=0
(13)
where A,B,C,D are the function of system parameters and vrv,0, and are defined as
A=krcg03,B=3krcσ2vrvg03+3krcg02krqg03,C=6krcg02σ2vrv2krqg02+3krcg03σ22vrv2+krg03+3krcg02krqg03σ2vrv,D=(σ2vrvg0+1)(krcg02σ22vrv2krqg02σ2vrv+krg02+g02+2krcσ2vrvg0krqg0+krc)
(14)
Therefore, depending on the coefficients of the cubic equation, we may have either only one set of fixed-points (one real with two complex conjugates roots of Eq. (13)) or three sets of fixed-points (three distinct real roots of Eq. (13)). For simplicity, we consider the scenario where Eq. (13) has only one real and two complex roots, or three repeated roots so that the analytical analysis will be carried out based on a set of unique real fixed-points. To ensure that Eq. (13) will have only one real root, the discriminant of polynomials should be less than zero. The discriminant for the cubic equation (Eq. (13)) is defined as
Ds=27A2D2+18ABCD4AC34B3D+B2C2
(15)

This solvability condition further puts the restriction on the ranges of krq,krc to get unique real fixed-points for the given value of kr. Therefore, we have ensured that the numerical values of krq and krc are chosen such that we get unique real fixed-points of the system.

Next, a small parameter ϵ(ϵ1) is introduced in Eq. (11) through following transformation
xi(τ)=xis+ϵyi(τ),(fori=1,2,,5)
(16)
where yi(τ)s are O(1) perturbations. Thus, Eq. (11) can be written as
y˙1=y2
(17a)
y˙2=y1kr(y1y3)2ζy22κ(y2y4)+h1h2(y1y3)+3h1krq(y1y3)ϵ(y32h22y3y1h2+y12h2)ϵ2krc(y1y3)3
(17b)
y˙3=y4
(17c)
y˙4=αy4h4+g0y5ασ1vrvσ0+2κα(y2y4)+y1krαy3krαy5ασ0+(3y1αh123y3αh12)krc+(2y1αh1+2y3αh1)krq+ϵy42(ασ1σ2vrvh7+h6h7)3krc(y3y1)2αh1+krq(y3y1)2α+y4y5ασ1σ0h5+ϵ2(y43ασ1h0h8y42y5ασ1σ0h7+krcα(y3y1)3)
(17d)
y˙5=vrvg0σ0y5g1vrvy4h0ϵ(g0y5y4σ0+y5y4vrvg1σ0+y42vrvg2h0+y42g1h0)ϵ2(y5y42σ0vrvg2+y43vrvg3h0+y5y42σ0g1+y43g2h0)
(17e)

where h0=1g0,h1=h0+xs3+σ2vrv,h2=3krch1+krq,h3=h1xs3σ2vrv,h4=σ2+σ1g1h3vrv,h5=g0+vrvg1,h6=ασ1(h1xs3),h7=vrvg2+g1, and h8=g2+g3vrv. Since nonlinearities in these equations appear as coefficients of higher orders of ϵ(>0), the unperturbed system can be obtained by setting ϵ=0 in Eq. (17) for the linear stability analysis.

3 Linear Stability of Fixed-Points and Existence of Critical Points

The fixed-point's linear stability, along with the existence of Hopf bifurcation in the system, is investigated in this section. As mentioned earlier, The linearized system of equations can be obtained by setting ϵ=0 in Eq. (17) to get
y˙1=y2
(18a)
y˙2=y1kr(y1y3)2ζy22κ(y2y4)+h1h2(y1y3)+3h1krq(y1y3)
(18b)
y˙3=y4
(18c)
y˙4=mry4h4+g0y5mrσ1vrvσ0+2κmr(y2y4)+y1krmry3krmry5mrσ0+3y1mrh123y3mrh12krc+(2y1mrh1+2y3mrh1)krq
(18d)
x˙5=vrvg0σ0y5g1vrvy4h0
(18e)
Accordingly, the characteristic equation for the system is obtained by substituting the synchronous solution for yi(for i=1.5) in the form of
(y1(τ)y2(τ)y3(τ)y4(τ)y5(τ))=(y10y20y30y40y50)eλτ
(19)
into Eq. (18) and focusing on the nontrivial solutions for the system to get
λ5+f1λ4+f2λ3+f3λ2+f4λ+f5=0
(20)
where fi(i=1,2,,5) are defined in Appendix-B and are functions of control gain parameters, the reference signal, and the primary system's characteristics. Note that the fixed-point's stability depends on the roots of the above characteristic equation for the given set of operational and system parameters. If all five roots of the above characteristic equation have a negative real part, the system is considered stable in a linear regime otherwise unstable. Since we can only vary the gain parameters and the reference signal for a given set of system parameters, the system's stability depends on these operational parameters. Therefore, there is a threshold/critical value for these operational parameters, which corresponds to a change in the system's stability. These threshold/critical points correspond to the Hopf bifurcation in the system and can be identified by setting λ=iω with ω>0 in the- characteristic equation (Eq. (20)). This substitution leads to a set of two algebraic equations by equating real and imaginary parts to zero,
f1ω4f3ω2+f5=0
(21)
ω5f2ω3+f4ω=0
(22)
For the current study, we select ζ (differential gain) and vrv (velocity signal) as our operational parameters. Accordingly, we solve the above simultaneous equations for ζ and vrv. Due to the complexity involved in fis(i=1,2,3,4,5), we solve these two equations numerically to get the threshold/critical values of ζ and vrv, i.e., ζi,cr and vrv,cr. Furthermore, the solution of the linearized equations (Eq. (18)) at the critical point is given by
y(τ)=A1r1eiωτ+A2r2eiωτ
(23)
In Eq. (23),y(τ)=[y1(τ),y2(τ),y3(τ),y4(τ),y5(τ)]T and r1 and r2 are the complex conjugate right eigenvectors of the characteristic matrix corresponding to eigenvalues λ=iω and λ=iω, respectively. Also, it should be noted that for y(τ) to be a real-valued vector, A1 and A2 have to be complex conjugate constants. The right eigenvector r1 corresponding to λ=iω is defined as
r1=[1iωRe1+im1Re2+im2Re3+im3]
(24)

Since the expressions for Rem and Imm(m=1,2,3) are lengthy, and hence, not reported in the work for the sake of brevity. In the next step, we analyze the system using the perturbation method, more specifically, the method of multiple scales (MMS).

4 The Method of Multiple Scales

The linear stability analysis of our system in Sec. 3 determines the fixed-points' local stability for a given set of operational and system parameters. For a given set of parameters, if a small perturbation to a fixed-point dies out with time, then it is locally stable; however, if it increases with time, then it is globally unstable. The sensitivity of fixed-points toward initial perturbation in a locally stable point and its time evolution depends on the nature of the existing system's nonlinearity. If all perturbations die out with time, irrespective of their magnitude, then a locally stable point is considered as a globally stable point for the fixed-point. However, the small perturbation dies out, and the large perturbation settles down to a limit cycle close to the critical point in a locally stable point. A globally stable point is different from a locally stable point, which further leads to the existence of bistable regions. Since such a phenomenon relies on the system's nonlinearity, linear stability analysis is insufficient to analyze the bistable regime. Therefore, we perform a thorough nonlinear analysis of the system at Hopf points to establish the globally stable region of the fixed-points.

For the current analysis, we use a perturbation method, more specifically, the method of multiple scales, to determine the amplitudes and locations of limit cycles along with the nature of the Hopf bifurcation at the Hopf point. To begin, we first define different time scales as
T0=τ,T1=ϵτ,T2=ϵ2τ,
(25)
Accordingly, the time-derivative operators also get modified to
ddτ=D0+ϵD1+ϵ2D2+O(ϵ3)
(26)
d2dτ2=D0,0+2ϵD0,1+ϵ2(2D0,2+D1,1)+O(ϵ3)
(27)
where Dn=Tn and Dm,n=2TmTn. In the next step, we assume the solution for Eq. (17) as a series of ϵtillO(ϵ2) as
y(τ)=y0(T0,T1,T2)+ϵy1(T0,T1,T2)+ϵ2y2(T0,T1,T2)=y0+ϵy1+ϵ2y2
(28)
with
y(τ)=y1(τ)y2(τ)y3(τ)y4(τ)y5(τ),yj(T0,T1,T2)=[y1,m(T0,T1,T2)y2,m(T0,T1,T2)y3,m(T0,T1,T2)y4,m(T0,T1,T2)y5,m(T0,T1,T2)]j=0,1,2
(29)
Next, we select the nondimensional reference velocity vrv as our bifurcation parameter, and perturb vrv as
vrv=vrv,cr+ϵ2k1
(30)
where vrv,cr is the value of vrv at the Hpof/critical point with ζ=ζcr.k1 is O(1) quantity and is chosen such that perturbed vrv remains in the unstable region. Also, due to the perturbation in vrv, the terms depending on vrv will get perturbed as well in the form of hi(vrv)=hi(vrv,cr)+ϵ2dhidvrvvrv=vrv,crk1+O(ϵ4)hi,cr+ϵ2hi,pk1. Next, we substitute Eqs. (26)(30) in Eq. (17), expand in Taylor series for smaller values of ϵ, and equate the coefficients of different orders of ϵ to zero to get O(ϵ0)
D0x1,1x2,1=0
(31a)
D0x2,1+2(x2,1x4,1)κ(x3,1x1,1)kr+h1,cr(x3,1x1,1)krq+x1,1+(x3,1x1,1)h1,crh2,cr+2ζx2,1=0
(31b)
D0x3,1x4,1=0
(31c)
D0x4,12mr(x2,1x4,1)κ+mr(x3,1x1,1)kr2mrh1,cr(x3,1x1,1)krq+mrσ0x5,1mrσ1σ0x5,1vrv,cr3g2,cr+mrσ1σ0x5,1vrv,crh5,cr3mrkrcx1,1h1,cr2+3mrkrcx3,1h1,cr2+mrx4,1h4,cr+mrσ1σ0x5,1vrv,cr2h7,cr=0
(31d)
D0x5,1+x4,1vrv,crg1,crh0,cr+σ0x5,1vrv,crg0,cr=0
(31e)
O(ϵ1)
D0x1,2x2,2=D1x1,1
(32a)
D0x2,2+2(x2,2x4,2)κ(x3,2x1,2)kr+h1,cr(x3,2x1,2)krq+x1,2+(x3,2x1,2)h1,crh2,cr+2ζx2,2=D1x2,1+h2,cr(x1,1x3,1)2
(32b)
D0x3,2x4,2=D1x3,1
(32c)
D0x4,22mr(x2,2x4,2)κ+mr(x3,2x1,2)kr2mrh1,cr(x3,2x1,2)krq+mrσ0x5,2+mrσ1σ0x5,2vrv,cr2h7,crmrσ1σ0x5,2vrv,cr3g2,cr+mrσ1σ0x5,2vrv,crh5,cr3mrkrcx1,2h1,cr2+3mrkrcx3,2h1,cr2+mrx4,2h4,cr=mr(x1,1x3,1)2krqD1x4,1mrσ1σ0x4,1x5,1h5,crx4,12mrσ1σ2vrv,crh7,cr+x4,12h7,crh6,cr3mrkrch1,cr(x3,1x1,1)2
(32d)
D0x5,2+x4,2vrv,crg1,crh0,cr+σ0x5,2vrv,crg0,cr=x4,12vrv,crg2,crh0,crσ0x4,1x5,1g0,crD1x5,1σ0x4,1x5,1vrv,crg1,crx4,12g1,crh0,cr
(32e)
O(ϵ2)
D0x1,3x2,3=D2x1,1D1x1,2
(33a)
D0x2,3+2(x2,3x4,3)κ(x3,3x1,3)kr+h1,cr(x3,3x1,3)krq+x1,3+(x3,3x1,3)h1,crh2,cr+2ζx2,3=(x1,1x3,1)k1h1,pkrq3krc(x1,1x3,1)3+(x1,1x3,1)h1,crk1h2,p+2x1,2x3,1h2,crD1x2,2D2x2,1(x1,1x3,1)k1h1,ph2,cr2x3,1x3,2h2,cr2x1,1x1,2h2,cr+2x1,1x3,2h2,cr
(33b)
D0x3,3x4,3=D2x3,1D1x3,2
(33c)
D0x4,32mr(x2,3x4,3)κ+mr(x3,3x1,3)kr2mrh1,cr(x3,3x1,3)krq+mrσ0x5,3+mrx4,3h4,cr+mrσ1σ0x5,3vrv,cr2h7,crmrσ1σ0x5,3vrv,cr3g2,cr+mrσ1σ0x5,3vrv,crh5,cr3mrkrcx1,3h1,cr2+3mrkrcx3,3h1,cr2=2x4,1x4,2h7,crh6,cr+mrσ1σ0x4,12x5,1h7,crD1x4,2D2x4,1x4,13h8,cr+krcmr(x1,1x3,1)3+2mr(x3,1x3,2x1,1x3,2x1,2x3,1+x1,1x1,2x1,1k1h1,p+x3,1k1h1,p)krqmrσ1σ0x4,1x5,2h5,crmrσ1σ0x4,2x5,1h5,cr+3mrσ1σ0x5,1vrv,cr2k1g2,cr2mrσ1σ0x5,1vrv,crk1h7,crmrσ1σ0x5,1vrv,cr2k1h7,p+mrσ1σ0x5,1vrv,cr3k1g2,pmrσ1σ0x5,1vrv,crk1h5,p6mrkrcx1,1x1,2h1,cr+6mrkrcx1,1x3,2h1,cr+6mrkrcx1,2x3,1h1,cr6mrkrcx3,1x3,2h1,crmrσ1σ0x5,1k1h5,cr6mrkrcx3,1h1,crk1h1,p+6mrkrcx1,1h1,crk1h1,p2x4,1x4,2mrσ11σ2vrv,crh7,crmrx4,1k1h4,cr
(33d)
D0x5,3+x4,3vrv,crg1,crh0,cr+σ0x5,3vrv,crg0,cr=σ0x4,2x5,1g0,crD2x5,1σ0x4,1x5,2g0,crD1x5,2x4,1vrv,crk1g1,ph0,crσ0x5,1k1g0,crσ0x4,12x5,1g1,cr2x4,1x4,2g1,crh0,crσ0x4,12x5,1vrv,crg2,crx4,13g2,crh0,cr2x4,1x4,2vrv,crg2,crh0,crσ0x5,1vrv,crk1g0,pσ0x4,1x5,2vrv,crg1,crx4,13vrv,crg3,crh0,crσ0x4,2x5,1vrv,crg1,crx4,1k1g1,crh0,crx4,1vrv,crg1,crk1h0,p
(33e)
Since the equations corresponding to O(ϵ0) (Eq. (31)) are similar to the linearized equations, i.e., Eq. (18) with critical/Hopf control parameters, the solution for Eq. (31) can be written in the form of the generalized right-eigen vector as
y0(T0,T1,T2)=A1(T1,T2)r1eiωT0+A2(T1,T2)r2eiωT0
(34)
In the above-assumed solution form, A1 and A2 are complex and conjugate functions. To get the closed form solution for A1(T1,T2) and A2(T1,T2), we proceed to the higher order of ϵ. Substitution of Eq. (34) in O(ϵ1) equations, i.e., Eq. (32) leads to the appearance of e2iωT0,e2iωT0,eiωT0, and eiωT0 in the resultant equations. However, the appearance of secular terms eiωT0 and eiωT0 give rise to linear growth in the solution for y1, and have to be eliminated from the corresponding equations to get a bounded solution for y1. In particular, this step requires that the dot product of the vector consisting of the coefficients of eiωT0(eiωT0) with the generalized left eigenvectors corresponding to λ=iω(λ=iω) should be zero [31]. The generalized left eigenvectors l1 for the characteristic matrix corresponding to eiωT0 is
l1=[1Lre1+iLim1Lre2+iLim2Lre3+iLim3Lre4+iLim4.]
(35)
Since the expressions for the components of left eigenvectors are lengthy and involved, these are not reported here for the sake of brevity. Furthermore, the left-eigen vector l2 for the eigenvalue λ=iω is the complex conjugate of l1, and hence, not reported in paper. The coefficient vectors u1 and u2 corresponding to eiωT0 and eiωT0, respectively are
u1=A1(T1,T2)T1[1iωRe1+iI1Re2+im2Re3+im3],u2=A2(T1,T2)T1[1iωRe1im1Re2im2Re3im3]
(36)
It can be easily observed that u1 and u2 are complex conjugates of each other and have the same coefficients as the generalized right eigenvector. This observation can be further justified by the appearance of only quadratic nonlinear terms in O(ϵ1) equations. These quadratic nonlinear terms further give rise to either eiωT0(eiωT0) or constant terms, and hence, do not contribute to secular terms. As discussed above, to remove secular terms corresponding to eiωT0 at O(ϵ1) we use the solvability condition l1u1=0, which further implies
A1(T1,T2)T1[1ωLim1+Lre2Re1Lim2Im1+Lre3Re2Lim3Im2+Lre4Re3Lim4Im3]+i(ωLre1+Lre2Im1+Lim2Re1+Lre3Im2+Lim3Re2+Lre4Im3+Lim4Re3)]=0
(37)
A1(T1,T2)T1=0.
(38)
Equation (38) follows from Eq. (37) as the term inside the bracket in Eq. (37) will not be zero for the given values of parameters at the Hopf point. Furthermore, the removal of secular terms corresponding to eiωT0 using the solubility condition also leads to
A2(T1,T2)T1=0
(39)
However, these results (Eqs. (38) and (39)) suggest that A1 and A2 do not depend on T1 and are only the functions of T2. Furthermore, to obtain the nonzero solutions of A1 and A2, we proceed to the next order of ϵ, i.e., equations corresponding to O(ϵ2), for which the solution at the O(ϵ1), i.e., y1 is needed (as evident from Eq. (33)). To get the solution for y1, we substitute A1T1=0 and A2T1=0, along with the assumed form of solution for y0 (Eq. (34)) in Eq. (32), and use the Harmonic balance method. For this, we assume the following form of the solution for y1
y1(T0,T1,T2)=A12(T2)B11e2iωT0+A22(T2)B22e2iωT0+A1(T2)A2(T2)B12,
(40)
where coefficient vector B11,B22 and B12 are defined as
B11=[b11b12b13b14b15],B22=[b21b22b23b24b25],andB12=[b31b32b33b34b35]
(41)
We are substituting Eq. (40) in Eq. (32) and collecting e2iωT0,e2iωT0 and constant terms, we get 15 simultaneous algebraic equations in terms of bmn (for m=1,2,3 and n=1.5) . Since the closed-form expressions for bmn are lengthy, we do not report these expressions for brevity. We substitute y0 and y1 in terms of A1(T2) and A2(T2) in O(ϵ2) equations (Eq. (33)) and collecting the coefficients of secular terms, i.e., eiωT0 and eiωT0. If V1 and V2 are the secular coefficient vectors corresponding to eiωT0 and eiωT0, respectively, then secular terms from Eq. (33) can be removed by l1V1=0 and l2V2=0. These solvability conditions lead to two first-order complex conjugate differential equations governing the slow time evolution of A1 and A2. Therefore, we switch to polar coordinates by substituting
A1(T2)=R(T2)eiϕ(T1)2,andA2(T2)=R(T2)eiϕ(T1)2
(42)
into the complex differential equation resulting from l1V1=0. We separate the resultant equation into real and imaginary parts and solve for R(T2)/T2 and ϕ(T2)/T2 as
R(T2)T2=p11k1R+p12R3
(43)
ϕ(T2)T2=p21k1+p22R2
(44)
where p11,p12,p21, and p22 are functions of system parameters, ζcr,vrv,cr, and ω. As p11,p12,p21, and p22 are very long expressions of the above-mentioned parameters, we do not report these in the current work for brevity. The equations governing the R and ϕ in the original time scale (τ) can be expressed using Eq. (26) as
dRdτ=ϵRT1+ϵ2RT2=ϵ2(p11k1R+p12R3)
(45a)
ϕτ=ϵϕT1+ϵ2ϕT2=ϵ2(p21k1+p22R2)
(45b)

Further, xi(τ) can be obtained by utilizing Eqs. (16), (28), (34), (42), and (45).

5 Results and Discussions

In this section, we examine the analytical results presented in Secs. 3 and 4 through numerical simulations. For our numerical simulation, we use the parameter values listed in Table 1 . We first analyze the linear stability of the system in the parametric space of ζ and vrv followed by the validation of our analytical formulation. Later on, by utilizing our analytical findings, we present the different regions of sub and supercritical Hopf bifurcation on linear stability boundaries. Finally, a detailed bifurcation analysis is presented.

Table 1

Numerical values of system's parameters [32]

mt(kg)1.5kp2e4
mb(Ns/m)0.75X0(m)0.0007353
σ0*(N/m)2.2e6σ1*,σ2*(Ns/m)237,14.25
fc*(N)5.1fs*(N)6.5
ω0(rad/s)115.5κ0.001
σ0110σ11.37
σ20.0823fs0.44
fc0.35a2.5
mt(kg)1.5kp2e4
mb(Ns/m)0.75X0(m)0.0007353
σ0*(N/m)2.2e6σ1*,σ2*(Ns/m)237,14.25
fc*(N)5.1fs*(N)6.5
ω0(rad/s)115.5κ0.001
σ0110σ11.37
σ20.0823fs0.44
fc0.35a2.5

5.1 Linear Stability Curves.

In this section, we present the effect of the nonlinear components of friction isolator on the linear stability of the system. For this step, we plot the stability curves for different combinations of krq and krc on the operational parameter region of ζvrv and are shown in Figs. 24. For ease of understanding, the unstable and stable regions are denoted by “U” and “S,” respectively.

Fig. 2
Comparison of stability region (a) for different values of cubic nonlinearity and krq=0.2, and (b) for different values of quadratic nonlinearity and krc=0.2 in nonlinear FI with linear FI. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5,kr=0.5,krq=0.2 and mr=2.
Fig. 2
Comparison of stability region (a) for different values of cubic nonlinearity and krq=0.2, and (b) for different values of quadratic nonlinearity and krc=0.2 in nonlinear FI with linear FI. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5,kr=0.5,krq=0.2 and mr=2.
Close modal
Fig. 3
Comparison of stability region (a) in the absence of quadratic stiffness (krq=0) for different values of cubic nonlinearity (b) in the absence of cubic stiffness (krc=0) for different values of cubic nonlinearity in nonlinear FI with linear FI. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35, κ=0.001,a=2.5,kr=0.5, and mr=2.
Fig. 3
Comparison of stability region (a) in the absence of quadratic stiffness (krq=0) for different values of cubic nonlinearity (b) in the absence of cubic stiffness (krc=0) for different values of cubic nonlinearity in nonlinear FI with linear FI. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35, κ=0.001,a=2.5,kr=0.5, and mr=2.
Close modal
Fig. 4
Comparison of stability region (a) in the absence of cubic and linear stiffness (krc=kr=0) for different values of quadratic nonlinearity and (b) in the absence of quadratic and linear stiffness krq=kr= 0) for different values of cubic nonlinearity in nonlinear FI with linear FI. Other parameters are σ0=110, σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5, and mr=2.
Fig. 4
Comparison of stability region (a) in the absence of cubic and linear stiffness (krc=kr=0) for different values of quadratic nonlinearity and (b) in the absence of quadratic and linear stiffness krq=kr= 0) for different values of cubic nonlinearity in nonlinear FI with linear FI. Other parameters are σ0=110, σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5, and mr=2.
Close modal

As mentioned earlier, it is mathematically challenging to get analytical expressions for vrv,cr and ζcr, hence, we obtain stability boundaries numerically by solving Eqs. (21) and (22) along with Eq. (13) for the varying values of frequency ω in a range ω∈(ω1, ω2). Since ω1 and ω2 are functions of system parameters, their numerical values vary from one case to another. On solving Eqs. (21) and (22) along with Eq. (13), for a given range of frequency, we get negative values of ζ and vrv. However, as negative values of the control gain and the reference signal are not feasible, we plot the stability curves for the positive values of parameters.

From Figs. 24, we can easily observe that, compared to the case of linear FI, the inclusion of quadratic and cubic nonlinearities in the FI increases the fixed-point's stability significantly. This observation further implies that the nonlinearities in FI support a wider range of stable operating conditions. However, the relative effects of quadratic/cubic nonlinearity on the stability region for a given value of cubic/quadratic nonlinearity in various scenarios are different. For example, Fig. 2 shows the stability boundaries for different values of cubic/quadratic nonlinearity for a given nonzero value of quadratic/cubic and linear stiffness whereas Fig. 3 shows the stability boundaries with different values of cubic/quadratic nonlinearity in the absence of quadratic/cubic nonlinearities and nonzero linear stiffness, and Fig. 4 shows stability boundaries for different values of cubic/quadratic nonlinearity in the absence of quadratic/cubic and linear stiffness. From Figs. 2 and 3, it can be observed that the relative effects of cubic and quadratic nonlinearities on the fixed-point's stability are almost identical with or without the other component of nonlinearity. This further implies that the rate of increase in stability with the increase in cubic nonlinearity is approximately the same as with the increase in quadratic nonlinearity. However, we emphasize that the overall stability at the given value of cubic nonlinearity (Fig. 3(a)) is higher than that at quadratic nonlinearity (Fig. 3(b)). Furthermore, from Fig. 4, it can be observed that although the rate of increase in stability with quadratic nonlinearity (Fig. 4(b)) is much higher than cubic nonlinearity, the overall stability boundary for a given value of cubic nonlinearity is significantly larger than the case of quadratic nonlinearity. These observations further suggest that increasing the cubic stiffness of FI is more beneficial than increasing the quadratic stiffness. Having established the effect of nonlinear stiffness on the fixed-point's stability, we analyze the Hopf bifurcation on the stability curves using analytical results obtained by MMS. However, before this step, we must validate our analytical results, which can be done by comparing them with numerical simulations and presented next.

5.2 Validation of Method of Multiple Scales.

To evaluate the accuracy of the MMS, we compare the solution of the system obtained from the slow-flow equations to the one obtained from Eq. (11) using the matlabode solver “ode45.” We first present the time response of the motion stage with nonlinear FI for two sets of operational parameters close to the Hopf point. In particular, we respectively choose two nondimensional reference velocities with a smaller and a larger value, i.e., vrv=0.0495<vrv,cr=0.05,ζcr=0.11873, and vrv=0.09<vrv,cr=0.1,ζcr=0.09502. Since both sets of parameters are in the unstable regime, we obtain a gradually increasing periodic response of different amplitudes as shown in Fig. 5. We emphasize that these time responses are shifted to the origin set at the fixed-point (Eq. (12)). From Fig. 5, it can be easily observed that the analytical solution of the system from MMS exhibits an excellent match with the numerical solution of the system. We repeat the same steps for the motion stage with linear FI and two sets of operational parameters, viz., vrv=0.0495<vrv,cr=0.05,ζcr=0.24562 and vrv=0.09<vrv,cr=0.1,ζcr=0.19082117. The results are shown in Fig. 6 and we observe a good match between the two approaches for the motion stage with linear FI. Hence, our analytical solutions (Eq. (45)) are valid.

Fig. 5
Comparison of analytical and numerical simulation for (a) ζcr=0.11873,vrv=0.0495<vrv,cr= 0.05, (b)ζcr=0.09502,vrv=0.09<vrv,cr=0.1, PD controlled motion stage with nonlinear friction isolator. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5, mr=2,krq=0.2,krc=0.2, and kr=0.5.
Fig. 5
Comparison of analytical and numerical simulation for (a) ζcr=0.11873,vrv=0.0495<vrv,cr= 0.05, (b)ζcr=0.09502,vrv=0.09<vrv,cr=0.1, PD controlled motion stage with nonlinear friction isolator. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5, mr=2,krq=0.2,krc=0.2, and kr=0.5.
Close modal
Fig. 6
Comparison of analytical and numerical simulation for (a) ζcr=0.24562,vrv=0.0495<vrv,cr=0.05, (b)ζcr=0.19082117,vrv=0.09<vrv,cr=0.1, PD controlled motion stage with linear friction isolator. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5,mr=2, and kr=0.5.
Fig. 6
Comparison of analytical and numerical simulation for (a) ζcr=0.24562,vrv=0.0495<vrv,cr=0.05, (b)ζcr=0.19082117,vrv=0.09<vrv,cr=0.1, PD controlled motion stage with linear friction isolator. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5,mr=2, and kr=0.5.
Close modal

5.3 Hopf Bifurcation.

In this section, we analyze the different regions of super and subcritical Hopf bifurcation on the stability lobes. When the system changes its stability through the Hopf bifurcation, the fixed-points settle down to a limit cycle close to Hopf point. Furthermore, the location of a limit cycle with respect to the Hopf point decides the nature of Hopf bifurcation. More specifically, in the case of supercritical Hopf bifurcation, these limit cycles exist in the unstable region only, whereas the existence of limit cycles close to the Hopf point in the stable regime signifies subcritical Hopf bifurcation. We emphasize that the presence of supercritical Hopf bifurcation leads to the fixed-point's global stability of the stable region, whereas subcritical bifurcation leads to a bistable region in the system. Therefore, it is an essential step toward the understanding of the criticality of Hopf bifurcation on the stable curves.

From Figs. 24, we can observe that the location of an unstable region with respect to a Hopf point is not uniform in case of perturbation in vrv,cr, i.e., the unstable vrv is higher for vrv,cr of low values and lower for vrv,cr of high values. This further implies that the sign of k1 in Eq. (30) will vary, i.e., positive for low values of vrv,cr and negative for high values of vrv,cr. However, in the case of perturbation in ζcr, the unstable region always lies on the left of the Hopf point, which further leads to a consistent sign (-ve) for k1. Therefore, for the sake of simplicity in determining the transition points from super to subcritical Hopf bifurcation or vice versa, we perturb ζ as ζ=ζcrϵ2k1, and follow the procedure mentioned in Sec. 4 to get
dRdτ=ϵ2(q11k1R+q12R3)
(46a)
ϕτ=ϵ2(q21k1+q22R2)
(46b)
where q11,q12,q21, and q22 depend on the system characteristics, critical operational parameters, and ω. Furthermore, the amplitude of the periodic solutions near the Hopf points can be estimated by setting R˙=0 in Eq. (46a), i.e., the nontrivial solution of Eq. (46a) and is given by
R=q11k1q12
(47)

Equation (47) plays an essential role in determining the criticality of Hopf bifurcation. In Eq. (47), if q11q12 is negative then for R to be a real quantity k1 should be positive. This further implies that limit cycles will exist in a linearly unstable region only, and the Hopf bifurcation will be supercritical. However, if for another set of critical parameter values q11q12 get positive then for a real value of R,k1 should be negative. Hence, limit cycles exist in the linear stable region, and the Hopf-bifurcation will be subcritical. After determining the criteria for subcritical and supercritical Hopf bifurcations, we evaluate q11q12 at every Hopf point on the stability curve and decide the characteristic of Hopf bifurcation. Figures 7(a) and 7(b) show the characteristic of Hopf bifurcation on the stability boundary for the MBMS with nonlinear and linear FI, respectively. From both figures, we can observe that supercritical Hopf bifurcations occur at low values of vrv. At the same time, fixed-points lose stability through subcritical Hopf bifurcations for high values of vrv.

Fig. 7
Criticality of Hopf bifurcation in the motion stage with (a) with nonlinear FI (krc=0.2,krq=0.2) and (b) linear FI (krc=krq=0). Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44, fc=0.35,κ=0.001,a=2.5,mr=2, and kr=0.5.
Fig. 7
Criticality of Hopf bifurcation in the motion stage with (a) with nonlinear FI (krc=0.2,krq=0.2) and (b) linear FI (krc=krq=0). Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44, fc=0.35,κ=0.001,a=2.5,mr=2, and kr=0.5.
Close modal

At first glance on Figs. 7(a) and 7(b), it appears that there is no effect of nonlinearity of FI on the subcritical and supercritical Hopf bifurcation regions on the stability curve, and transition point remains the same with nonlinear and linear FI. Therefore, to demonstrate the effect of nonlinear parameters of FI (krq,krc) on the criticality of Hopf bifurcation, we plot q11q12 with different values of ζcr and vrv,cr for different sets of krq and krc, and the results are shown in Figs. 8(a) and 8(b), respectively.

Fig. 8
Variation of q11/q12 with (a) ζcr and (b) vrv,cr for different values of krq and krc. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5,mr=2, and kr=0.5.
Fig. 8
Variation of q11/q12 with (a) ζcr and (b) vrv,cr for different values of krq and krc. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5,mr=2, and kr=0.5.
Close modal

From Fig. 8(a), it can be noted that the inclusion of nonlinearity in FI reduces the range of ζcr corresponding to supercritical Hopf bifurcation. However, this can be further justified by the fact that the nonlinearity shrinks the unstable region by decreasing ζcr values which lead to a decrease in the effective range of ζcr for supercritical Hopf bifurcation. Instead, Figs. 8(a) and 8(b) provide more information about the effect of krq and krc on the criticality of Hopf bifurcation. From Figs. 8(a) and 8(b), we can easily observe that the inclusion of nonlinearity in FI can decrease or increase the region of supercritical Hopf bifurcation depending on the numerical values of krq and krc. The optimization of these values for a larger region of supercritical Hopf bifurcation and globally stable region is left for future work.

We emphasize that these analytical findings only provide the amplitude of limit cycles close to the Hopf point and characteristics of Hopf bifurcation on the stability boundaries. Therefore, to observe the global behavior of the system in the unstable region, we employ the numerical bifurcation analysis and present in the Sec. 5.4. Note that this step not only provides information about the large amplitude response of the system but also further verifies our analytical findings.

5.4 Bifurcation Analysis.

To perform the numerical bifurcation analysis for the motion stage with nonlinear and linear FI, we solve the system of ODEs given by Eq. (10) using matlabode solver “ode45.” The bifurcation plots, showing the extreme points of x1, i.e., the error amplitude of the motion stage (corresponding to x2=0), for the motion stage with nonlinear and linear FI have been shown in Figs. 9 and 10, respectively. The numerical bifurcation analysis can be performed by making either of the operational parameters a constant and varying the another. However, to show the existence of super- and subcritical Hopf bifurcations for lower and higher values of vrv, respectively, as observed in Sec. 5.3, we have fixed ζ and varied vrv. To plot these numerical bifurcation diagrams, we vary vrv in upward and downward directions so that the system loses and gains stability through Hopf bifurcation. For completeness, we have also plotted these numerical bifurcation diagrams for two different values of ζ. To get a better picture of the dynamics of the motion stage with nonlinear and linear FI, the bifurcation diagrams close to the Hopf points are shown in the inset of Figs. 9 and 10. From these numerical bifurcation diagrams, we can easily observe the existence of stable limit cycles with fixed-point solutions at higher values of vrv for a given value of ζ, which implies Hopf bifurcation is subcritical by nature. However, for lower values of vrv, stable limit cycles exist in the unstable region only, which indicates supercritical bifurcation. Both of these observations are consistent with our analytical findings in Sec. 5.3. Furthermore, in the case of the motion stage with nonlinear

Fig. 9
Numerical bifurcation diagram of motion stage with nonlinear FI with vrv as bifurcation parameter for (a) ζ=0.06 and (b) ζ=0.1. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44, fc=0.35,κ=0.001,a=2.5,mr=2,krq=0.2,krc=0.2, and kr=0.5.
Fig. 9
Numerical bifurcation diagram of motion stage with nonlinear FI with vrv as bifurcation parameter for (a) ζ=0.06 and (b) ζ=0.1. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44, fc=0.35,κ=0.001,a=2.5,mr=2,krq=0.2,krc=0.2, and kr=0.5.
Close modal
Fig. 10
Numerical bifurcation diagram of motion stage with linear FI with vrv as bifurcation parameter for (a) ζ=0.06 and (b) ζ=0.1. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44, fc=0.35,κ=0.001,a=2.5,mr=2, and kr=0.5.
Fig. 10
Numerical bifurcation diagram of motion stage with linear FI with vrv as bifurcation parameter for (a) ζ=0.06 and (b) ζ=0.1. Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44, fc=0.35,κ=0.001,a=2.5,mr=2, and kr=0.5.
Close modal

FI, the response amplitude for higher values of vrv is relatively smaller than the ones corresponding to the motion stage with linear FI. For a better understanding of the dynamics of a PD-controlled motion stage with nonlinear and linear FI, the zoomed views of Figs. 9 and 10 have been shown in Figs. 11 and 12, respectively. For the sake of brevity, we only present these zoomed views for differential gains of higher values, i.e., ζ=0.1. The corresponding representative phase portraits for different values of vrv have been shown inside these zoomed figures. From Figs. 11 and 12, we can easily observe that in both cases, close to Hopf points, stable period-1 solutions lose stability through period-doubling bifurcation. This further leads to the appearance of period-2 solutions, which can also be observed from the phase portraits (Figs. 11(a) and 11(d) and 12(a) and 12(d)). Furthermore, in the case of nonlinear FI, the system exhibits only period-4 solutions away from the Hopf points, and there is no exchange in the stability of limit cycles away from the Hopf points (Figs. 11(b) and 11(c)). However, in the case of linear FI, apart from the coexistence of period-1 and period-2 solutions (as can be seen by phase portraits for vrv=0.02 in Fig. 12(a)), there is a continuous exchange of stability between period-1 and period-2 solutions as shown in Figs. 12(a),12(d). Also, when comparing Figs. 11 and 12 in terms of subplots (i) and (ii), we observe that the branch of stable period1 solutions close to Hopf point is significantly smaller in case of linear FI when compared to the case of nonlinear FI. This observation further signifies the importance of nonlinear FI over linear FI.

Fig. 11
Zoomed view of Fig. 9(b) for different range of vrv. Other parameters are σ0=110,σ1=1.37, σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5,mr=2,krq=0.2,krc=0.2, and kr=0.5.
Fig. 11
Zoomed view of Fig. 9(b) for different range of vrv. Other parameters are σ0=110,σ1=1.37, σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5,mr=2,krq=0.2,krc=0.2, and kr=0.5.
Close modal
Fig. 12
Zoomed view of Fig. 10(b) for different range of vrv. Other parameters are σ0=110,σ1=1.37, σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5,mr=2, and kr=0.5.
Fig. 12
Zoomed view of Fig. 10(b) for different range of vrv. Other parameters are σ0=110,σ1=1.37, σ2=0.0823,fs=0.44,fc=0.35,κ=0.001,a=2.5,mr=2, and kr=0.5.
Close modal

We perform the quantitative match between MMS results and numerical simulations for completeness. For this step, we use the fixed-arc-length continuation scheme [33] to get the branch of limit cycles close to the Hopf point and later compare it with the branch of limit cycles obtained using the slow-flow equation emerging from the Hopf point. These results are shown in Figs. 13(a) and 13(b) for MBMS with nonlinear and linear FI, respectively. Since the analysis has been performed close to the Hopf points, we can observe an excellent match between MMS and numerical simulations close to the Hopf point, which further verifies our analytical approach.

Fig. 13
Comparison of bifurcation diagram from numerical simulation and MMS with vrv as bifurcation parameter for the motion stage (a) with nonlinear FI (krq=0.2,krc=0.2,ζ=0.1), and (b) with linear FI (krq=0,krc=0,ζ=0.19). Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35, a=2.5,κ=0.001,kr=0.5,mr=2.
Fig. 13
Comparison of bifurcation diagram from numerical simulation and MMS with vrv as bifurcation parameter for the motion stage (a) with nonlinear FI (krq=0.2,krc=0.2,ζ=0.1), and (b) with linear FI (krq=0,krc=0,ζ=0.19). Other parameters are σ0=110,σ1=1.37,σ2=0.0823,fs=0.44,fc=0.35, a=2.5,κ=0.001,kr=0.5,mr=2.
Close modal

6 Conclusion

In this study, we analyzed the linear and nonlinear dynamics of a PD-controlled MBMS with a FI using analytical and numerical methods. Contrary to earlier studies where the nonlinearity in the friction isolator was ignored, the effects of nonlinearity from the friction isolator on the dynamics of the motion stage have been explored in this work. The dynamical effect of friction was captured through the LuGre friction model. A parametric study on the linear stability analysis revealed that compared to the linear friction isolator, the nonlinearity in the friction isolator increases the fixed-points' local stability in the operating parameter space of reference signal and differential gain. This further implied that the system's stability is underestimated when using the linear FI model. The nonlinearity in the FI should be considered in the modeling for a better prediction of steady operating conditions. The nonlinear analysis of MBMS with FI was performed analytically using MMS and numerical simulations. The analytical findings were verified by comparing them against numerical solutions, and a good match was observed. Both analytical and numerical simulations revealed the existence of supercritical and subcritical Hopf bifurcation. Furthermore, the parametric analysis on the criticality of Hopf bifurcation revealed the sensitivities of subcritical and supercritical Hopf bifurcation regions in terms of the nonlinearities of the friction isolator. On exploring the dynamics of MBMS with nonlinear FI in the unstable regime, we observed period-doubling bifurcations, period-4 solutions, and quasi-periodic solutions.

Finally, these findings suggest that the consideration of nonlinearity in the FI model is an essential step to get an accurate picture of global dynamics of the motion stage with a FI.

Acknowledgment

This work was funded by National Science Foundation (NSF) Award CMMI #1855390: Towards a Fundamental Understanding of a Simple, Effective, and Robust Approach for Mitigating Friction in Nanopositioning Stages, and CMMI #2000984: Nonlinear Dynamics of Pneumatic Isolators in Ultra-Precision Manufacturing Machines.

Funding Data

  • National Science Foundation (NSF) (Award No. CMMI #1855390; doi: 10.13039/100000001).

Conflict of Interest

The authors declare that they have no conflict of interest.

Data Availability Statement

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Mathematical Model

A.1 Expression Used in Eq. (11)

g0=1g(vrv)=1g,g1=1g2gvrv,g2=1g3[(gvrv)2g22gvrv2],g3=1g4[(gvrv)3ggvrv2gvrv2+g263gvrv3]

Linear Stability Analysis

B.1 Expressions Used in Eq. (20)

f1=(vrvσ0g0+2ζ+2κmr+2κ+mrh4)f2=vrv4g1h0mrσ1σ0g22mrh1krq+2ζvrvσ0g0+3mrkrch12+krmrvrvg1h0mrσ0+2κmrvrvσ0g0+2κmrh4vrv3g1h0mrσ1σ0h7+2κvrvσ0g0vrv2g1h0mrσ1σ0h5+krh1h2+4ζκmr+mrh4vrvσ0g0h1krq+2ζmrh4+1f3=2κvrv4g1h0mrσ1σ0g22κvrv3g1h0mrσ1σ0h72κvrv2g1h0mrσ1σ0h5+2ζvrv4g1h0mrσ1σ0g22ζvrv3g1h0mrσ1σ0h72ζvrv2g1h0mrσ1σ0h5+krvrvσ0g0+6ζmrkrch124ζmrh1krq+2ζkrmrh1h2mrh4h1krqmrh4+krmrh4+mrh4+2κmr+3mrkrch12vrvσ0g02mrh1krqvrvσ0g0+2κmrh4vrvσ0g02κvrvg1h0mrσ0+2ζmrh4vrvσ0g0+4ζκmrvrvσ0g02ζvrvg1h0mrσ0+vrvσ0g0+krmrvrvσ0g0h1krqvrvσ0g0h1h2vrvσ0g0f4=krvrv4g1h0mrσ1σ0g2krvrv3g1h0mrσ1σ0h7krvrv2g1h0mrσ1σ0h5+2ζkrmrvrvσ0g0+krmrh4vrvσ0g0krvrvg1h0mrσ0+3mrkrch122mrh1krqh1h2vrv4g1h0mrσ1σ0g2+h1h2vrv3g1h0mrσ1σ0h7+h1h2vrv2g1h0mrσ1σ0h5h1krqvrv4g1h0mrσ1σ0g2+krmr+h1krqvrv3g1h0mrσ1σ0h7+h1krqvrv2g1h0mrσ1σ0h5+mrh4vrvσ0g0+2κmrvrvσ0g0vrvg1h0mrσ0+vrv4g1h0mrσ1σ0g2vrv3g1h0mrσ1σ0h7vrv2g1h0mrσ1σ0h5+6ζmrkrch12vrvσ0g04ζmrh1krqvrvσ0g0+h1h2vrvg1h0mrσ0h1krqmrh4vrvσ0g0+h1krqvrvg1h0mrσ0h1h2mrh4vrvσ0g0f5=3mrkrch12vrvσ0g0+krmrvrvσ0g02mrh1krqvrvσ0g0

References

1.
Kim
,
W.-J.
, and
Trumper
,
D. L.
,
1998
, “
High-Precision Magnetic Levitation Stage for Photolithography
,”
Precis. Eng.
,
22
(
2
), pp.
66
77
.10.1016/S0141-6359(98)00009-9
2.
Kim
,
W. J.
,
Verma
,
S.
, and
Shakir
,
H.
,
2007
, “
Design and Precision Construction of Novel Magneticlevitation-Based Multi-Axis Nanoscale Positioning Systems
,”
Precis. Eng.
,
31
(
4
), pp.
337
350
.10.1016/j.precisioneng.2007.02.001
3.
Salapaka
,
S. M.
, and
Salapaka
,
M. V.
,
2008
, “
Scanning Probe Microscopy
,”
IEEE Control Syst. Mag.
,
28
(
2
), pp.
65
83
.10.1109/MSC.2007.914688
4.
Sebastian
,
A.
,
Pantazi
,
A.
,
Moheimani
,
S. R.
,
Pozidis
,
H.
, and
Eleftheriou
,
E.
,
2008
, “
Achieving Subnanometer Precision in a Mems-Based Storage Device During Self-Servo Write Process
,”
IEEE Trans. Nanotechnol.
,
7
(
5
), pp.
586
595
.10.1109/TNANO.2008.926441
5.
Altintas
,
Y.
,
Verl
,
A.
,
Brecher
,
C.
,
Uriarte
,
L.
, and
Pritschow
,
G.
,
2011
, “
Machine Tool Feed Drives
,”
CIRP Annals
,
60
(
2
), pp.
779
796
.10.1016/j.cirp.2011.05.010
6.
Futami
,
S.
,
Furutani
,
A.
, and
Yoshida
,
S.
,
1990
, “
Nanometer Positioning and Its Microdynamics
,”
Nanotechnology
,
1
(
1
), pp.
31
37
.10.1088/0957-4484/1/1/006
7.
Kim
,
M.-S.
, and
Kim
,
J.-H.
,
2011
, “
Design of a Gain Scheduled Pid Controller for the Precision Stage in Lithography
,”
Int. J. Precis. Eng. Manuf.
,
12
(
6
), pp.
993
1000
.10.1007/s12541-011-0132-6
8.
Liu
,
J. Y.
,
Yin
,
W. S.
, and
Zhu
,
Y.
,
2011
, “
Application of Adaptive Fuzzy Pid Controller to Nano-Scale Precision Motion Stage System
,”
Control Eng. China
,
18
(
2
), pp.
254
257
.
9.
Al-Bender
,
F.
, and
Swevers
,
J.
,
2008
, “
Characterization of Friction Force Dynamics
,”
IEEE Control Syst. Mag.
,
28
(
6
), pp.
64
81
.
10.
Marques
,
F.
,
Flores
,
P.
,
Pimenta Claro
,
J. C.
, and
Lankarani
,
H. M.
,
2016
, “
A Survey and Comparison of Several Friction Force Models for Dynamic Analysis of Multibody Mechanical Systems
,”
Nonlinear Dyn.
,
86
(
3
), pp.
1407
1443
.
11.
Armstrong-Hélouvry
,
B.
,
Dupont
,
P.
, and
De Wit
,
C. C.
,
1994
, “
A Survey of Models, Analysis Tools and Compensation Methods for the Control of Machines With Friction
,”
Automatica
,
30
(
7
), pp.
1083
1138
.10.1016/0005-1098(94)90209-7
12.
Hensen
,
R. H.
,
Van de Molengraft
,
M. J. G.
, and
Steinbuch
,
M.
,
2003
, “
Friction Induced Hunting Limit Cycles: A Comparison Between the Lugre and Switch Friction Model
,”
Automatica
,
39
(
12
), pp.
2131
2137
.10.1016/S0005-1098(03)00234-6
13.
Maeda
,
Y.
, and
Iwasaki
,
M.
,
2013
, “
Rolling Friction Model-Based Analyses and Compensation for Slow Settling Response in Precise Positioning
,”
IEEE Trans. Ind. Electron.
,
60
(
12
), pp.
5841
5853
.10.1109/TIE.2012.2229676
14.
Dong
,
X.
,
Yoon
,
D.
, and
Okwudire
,
C. E.
,
2017
, “
A Novel Approach for Mitigating the Effects of Pre-Rolling/Pre-Sliding Friction on the Settling Time of Rolling Bearing Nanopositioning Stages Using High Frequency Vibration
,”
Precis. Eng.
,
47
, pp.
375
388
.10.1016/j.precisioneng.2016.09.011
15.
Dejima
,
S.
,
Gao
,
W.
,
Katakura
,
K.
,
Kiyono
,
S.
, and
Tomita
,
Y.
,
2005
, “
Dynamic Modeling, Controller Design and Experimental Validation of a Planar Motion Stage for Precision Positioning
,”
Precis. Eng.
,
29
(
3
), pp.
263
271
.10.1016/j.precisioneng.2004.11.005
16.
Yao
,
J.
,
Deng
,
W.
, and
Jiao
,
Z.
,
2015
, “
Adaptive Control of Hydraulic Actuators With Lugre Model-Based Friction Compensation
,”
IEEE Trans. Ind. Electron.
,
62
(
10
), pp.
6469
6477
.10.1109/TIE.2015.2423660
17.
Wang
,
X.
, and
Wang
,
S.
,
2012
, “
High Performance Adaptive Control of Mechanical Servo System With Lugre Friction Model: Identification and Compensation
,”
J. Dyn. Syst., Meas., Control
,
134
(
1
), p.
011021
.10.1115/1.4004785
18.
Barreto S
,
J. C. L.
,
Conceicao
,
A. G. S.
,
Dorea
,
C. E. T.
,
Martinez
,
L.
, and
de Pieri
,
E. R.
,
2014
, “
Design and Implementation of Model-Predictive Control With Friction Compensation on an Omnidirectional Mobile Robot
,”
IEEE/ASME Trans. Mechatron.
,
19
(
2
), pp.
467
476
.10.1109/TMECH.2013.2243161
19.
Dong
,
X.
,
Liu
,
X.
,
Yoon
,
D.
, and
Okwudire
,
C. E.
,
2017
, “
Simple and Robust Feedforward Compensation of Quadrant Glitches Using a Compliant Joint
,”
CIRP Ann.
,
66
(
1
), pp.
353
356
.10.1016/j.cirp.2017.04.048
20.
Dong
,
X.
, and
Okwudire
,
C. E.
,
2018
, “
An Experimental Investigation of the Effects of the Compliant Joint Method on Feedback Compensation of Pre-Sliding/Pre-Rolling Friction
,”
Precis. Eng.
,
54
, pp.
81
90
.10.1016/j.precisioneng.2018.05.004
21.
Dong
,
X.
,
Okwudire
,
C.
,
Wang
,
J.
, and
Barry
,
O.
,
2019
, “
On the Friction Isolator for Precision Motion Control and Its Dynamics
,”
ASME
Paper No. DETC2019-98354.
22.
Gupta
,
S. K.
,
Wang
,
J.
, and
Barry
,
O. R.
,
2020
, “
Nonlinear Vibration Analysis of a Servo Controlled Precision Motion Stage With Friction Isolator
,”
Int. J. Non-Linear Mech.
,
126
, p.
103554
.10.1016/j.ijnonlinmec.2020.103554
23.
Canudas De Wit
,
C.
,
Olsson
,
H.
,
Astrom
,
K. J.
, and
Lischinsky
,
P.
,
1995
, “
A New Model for Control of Systems With Friction
,”
IEEE Trans. Autom. Control
,
40
(
3
), pp.
419
425
.10.1109/9.376053
24.
Johanastrom
,
K.
, and
Canudas-De-Wit
,
C.
,
2008
, “
Revisiting the Lugre Friction Model
,”
IEEE Control Syst. Mag.
,
28
(
6
), pp.
101
114
.10.1109/MCS.2008.929425
25.
Pennestrì
,
E.
,
Rossi
,
V.
,
Salvini
,
P.
, and
Valentini
,
P. P.
,
2016
, “
Review and Comparison of Dry Friction Force Models
,”
Nonlinear Dyn.
,
83
(
4
), pp.
1785
1801
.10.1007/s11071-015-2485-3
26.
Marques
,
F.
,
Flores
,
P.
,
Pimenta Claro
,
J. C.
, and
Lankarani
,
H. M.
, 2019, “
Modeling and Analysis of Friction Including Rolling Effects in Multibody Dynamics: A Review
,”
Multibody Syst. Dyn.
,
45
(
2
), pp.
223
244
.10.1007/s11044-018-09640-6
27.
Saha
,
A.
, and
Wahi
,
P.
,
2014
, “
An Analytical Study of Time-Delayed Control of Friction-Induced Vibrations in a System With a Dynamic Friction Model
,”
Int. J. Non-Linear Mech.
,
63
, pp.
60
70
.10.1016/j.ijnonlinmec.2014.03.012
28.
Saha
,
A.
,
2013
, “
Analysis and Control of Friction-Induced Vibrations by Time-Delayed Feedback
,” Ph.D. thesis,
Indian Institute of Technology
,
Kanpur, India
.
29.
Habib
,
G.
,
Detroux
,
T.
,
Viguie
,
R.
, and
Kerschen
,
G.
,
2015
, “
Nonlinear Generalization of Den Hartog's Equal-Peak Method
,”
Mech. Syst. Signal Process.
,
52–53
, pp.
17
28
.
30.
Habib
,
G.
, and
Kerschen
,
G.
,
2016
, “
A Principle of Similarity for Nonlinear Vibration Absorbers
,”
Phys. D: Nonlinear Phenom.
,
332
, pp.
1
8
.10.1016/j.physd.2016.06.001
31.
Balachandran
,
B.
, and
Nayfeh
,
A. H.
,
1992
, “
Cyclic Motions Near a Hopf Bifurcation of a Four-Dimensional System
,”
Nonlinear Dyn.
,
3
(
1
), pp.
19
39
.10.1007/BF00045469
32.
Wang
,
J.
,
Dong
,
X.
,
Barry
,
O. R.
, and
Okwudire
,
C.
,
2019
, “Friction-Induced Instability and Vibration in a Precision Motion Stage With a Friction Isolator,”
J. Vib. Control
, 28(15–16), pp. 1879–1893.10.1177/1077546321999510
33.
Wahi
,
P.
, and
Chatterjee
,
A.
,
2008
, “
Self-Interrupted Regenerative Metal Cutting in Turning
,”
Int. J. Non-Linear Mech.
,
43
(
2
), pp.
111
123
.10.1016/j.ijnonlinmec.2007.10.010