An Evolutionary Analytic Method of Multi-DOF Nonlinear Coupling Dynamic Model for Controllable Close-Chain Linkage Mechanism System

Abstract

The 2-DOF controllable close-chain linkage mechanism is investigated in this paper. Based on the characteristics of the multi-DOF nonlinear coupling dynamic equation of the system established by the finite element method, an analytic method of multiple-scales Newmark is presented after thinking about the method of perturbation and the method of numerical analysis. Firstly, the first-order approximate solution of the dynamic responses of the system at the time of 𝑡 is calculated by the multiple scales method. Then, taken the first-order approximate solution as the initialization of the generalized coordinate of the system, the stable dynamic response of the system is obtained by the implicit Newmark method. The simulation and experimental results are given in the end. The studies indicate that the method of multiple-scales Newmark is correct and practicable to study the dynamic characteristics of such kind of multi-DOF nonlinear coupling system.

1. Introduction

Multi-DOF controllable linkage mechanism, which can accurately actualize the given trajectory, velocity, and acceleration, will has a wide outlook of application in robots, automatic production lines, and so on [1–5]. The dynamic equation of such kind of linkage mechanism is nonlinear coupling time-variant from the dynamics modeling process of the linkage system. The numerical Newmark method can calculate the dynamic responses, but that method could not analyze internal relation between the dynamic characteristics and scale and electromagnetism parameters. Though the multiple-scales method can analyze the vibratory mechanism, it usually adopts first approximations and quadratic or higher order approximation is very complex [6–10]. An evolutionary analytic method of multiple-scales-Newmark method is firstly presented to study the dynamic characteristics of such kind of multi-DOF coupling system synthetically using the property that the implicit Newmark method can calculate the unconditioned stable solution when the Newmark parameters are specific constants the features of the multiple scales method.

2. Nonlinear Dynamic Equation of System

The 2-DOF controllable close-chain linkage mechanism is investigated in this paper. The analysis diagram of system is shown in Figure 1. Based on the air-gap field of nonuniform airspace of controllable motors of the linkage mechanism caused by the eccentricity of rotor, the controllable motor element [11] as shown in Figure 2, which defined the transverse vibration and torsional vibration of the controllable motors as its nodal displacement, was established. In the diagram, numbers 1, 2, 3, and 4 denote four nodes of the element. So the transverse vibration and the torsional vibration can be expressed by the generalized coordinate vector 𝐮₁=[𝑈₁𝑈₂𝑈₃𝑈₄]^𝑇, and 𝑋𝑌𝑍 is the coordinate system of controllable motor element.

Figure 1 

Diagram of the 2-DOF close-chain linkage mechanism.

Figure 2 

Diagram of the controllable motor element.

According to the mechatronics analysis dynamics, the air-gap eccentric vibration is shown in Figure 3, where point O is the inner circle geometric center of the motor stator, point O₁ is the outer circle geometric center of the rotor journal, point O₂ is the outer circle center of the journal under the deformation of shaft or bearing, the coordinate of point O₃(𝑥,𝑦) is the outer circle geometric center of the rotor, 𝛿 is the length of air gap, and 𝑒₁ is the air-gap eccentricity. 𝑒₀₁=𝑒²¹¹+𝑒²¹² is the static eccentricity, which is caused by the rotor gravity and mismachining tolerance of motor. 𝜀₀₁ is the rotational eccentricity, which is caused by centering error between the outer circle center of the journal and the outer circle geometric center of the rotor. The static eccentricity and the rotational eccentricity are considered at the same time in this paper. The static and rotational eccentricities also can be ignored, and only the eccentricity caused by rotor vibration can be considered in engineering application latterly. Then, 𝑒₁=√𝑥²+𝑦²,𝑥𝑦==𝑈𝑈₁₂+𝑒₁₁+𝜀₀₁cos𝛽₁,+𝑒₁₂+𝜀₀₁sin𝛽₁,(2.1) where 𝑢₁ and 𝑢₂ are the components of vibration eccentricity in the 𝑥 and 𝑦 direction, and 𝑒₀=𝑢²¹+𝑢²². 𝛽₁=(1−𝑠₁)𝜔₀₁𝑡 is the rotational angle of rotor with respect to the stator of controllable motor. 𝜔₀₁ is the synchronous speed of rotation of controllable motor. 𝑠₁is the slide ratio.

(a)

(b)

Figure 3 

Diagram of air-gap eccentric vibration of controllable motor.

According to the theory of electromechanics and electromechanical analysis dynamics, as far as controllable motor is concerned, the voltage and current between the two windings are asymmetric, that is to say, the controllable motor works in elliptic magnetic field. According to this real running state of controllable motor, the kinetic energy and potential energy of controllable motor can be obtained as follows:𝑇₁=12^𝑙0𝑚₁(𝑥)𝑊̇₁(𝑥,𝑡)²𝑑𝑥+12^𝑙0𝐽₀₁(𝑥)𝑉₁̇(𝑥,𝑡)²𝑑𝑥=12𝐮̇^𝑇1𝑉𝐌₁𝐮̇₁,(2.2)₁=12𝐮^𝑇1𝐊₁₁+𝐊₁₂𝐮¹+𝐞^𝑇1𝐊₁₁+𝐊₁₂𝐮¹+12𝐞^𝑇1𝐊₁₁+𝐊₁₂𝐞₁+𝐞₁+𝐮₁𝐤₀₁+𝑝𝑅⁰¹𝐿₀₁Λ₀₁2^02𝜋𝐹_+𝑠+𝐹_−𝑠cos𝜔₀₁𝑡−𝑝𝛼+𝐹_−𝑟cos𝜔₀₁𝑡+𝑝𝛼+𝐹_+𝑟cos𝜔₀₁𝑡−𝑝𝛼−𝜑₁₀cos𝜔₀₁𝑡+𝑝𝛼−𝜑₂₀²𝑑𝛼,(2.3) where 𝑊₁(𝑥,𝑡)and 𝑉₁(𝑥,𝑡) the transverse displacements and the angle of elastic torsion of any points in the controllable motor element, respectively (see the appendix), and 𝑙=𝑙₁+𝑙₂+𝑙₃ is the length of controllable motor shaft (see the appendix);𝑚₁(𝑥), including the rotor mass 𝑚₀ which is at 𝑥=𝑙₁, is the mass distribution function of the controllable motor shaft;𝐽₀₁(𝑥), including the moment of inertia of rotor𝐽₀which is at𝑥=𝑙₁, is the moment of inertia distribution function of the controllable motor shaft. 𝑅₀₁ is the inner radius of the controllable motor stator, 𝐿₀₁ is the effective length of the rotor, Λ₀₁=𝜇₀/𝜎 is the even air gap permeance of the controllable motor, 𝜇₀ is the magnetic permeability coefficient of air, 𝜎=𝑘_𝜇𝛿₀, 𝑘_𝜇 is saturation, 𝛿₀ is the uniform air-gap size, 𝑘_𝜇=1+𝛿_𝐹𝑒/𝑘₁𝛿₀, 𝑘₁ is the calculation air-gap coefficient of the even air gap, 𝛿_𝐹𝑒 is the equivalent air-gap of ferromagnetic materials, 𝜑₁₀ is the phase angle of the positive-sequence current of rotor lagging behind the positive-sequence current of stator, and 𝜑₂₀ is the phase angle of the negative-sequence current of rotor lagging behind the negative-sequence current of stator. 𝐹_+𝑠,𝐹_−𝑠,𝐹_+𝑟, and 𝐹_−𝑟 are the positive-sequence and negative-sequence components of the magnetomotive amplitude of stator and rotor respectively (see the appendix). 𝐌₁ is the mass 4×4 matrix of controllable motor element, 𝐊₁₁ is the stiffness 4×4 matrix of controllable motor element in connection with the structural parameters of the rotor, 𝐊₁₂ is the stiffness 4×4 matrix of controllable motor element in connection with the electromagnetic parameters of the rotor, and 𝐤₀₁ is the 4-order vector in connection with the eccentric motor (see the appendix). 𝐞=[𝑒₁₁+𝜀₀₁cos𝛽₁𝑒₁₂+𝜀₀₁sin𝛽₁00]^𝑇 is a matrix in connection with the static and rotary eccentricity of the controllable motor.

In general, the links of the 2-DOF controllable linkage mechanism are slim bars, so they are adapted to be simulated using beam element as shown in Figure 4. In dynamic analysis of the beam element, the coupling terms of the elastic motion and the rigid body motion in the Coriolis acceleration and transport acceleration are neglected in studying the absolute acceleration of any point in the beam element. In calculation of strain energy, the shearing deformation energy and yield deformation energy are also omitted. The material of the link is adopted as metal. Therefore, the kinetic energy and potential energy, respectively, are as follows: 𝑇₃=12^𝐿0𝜌𝐴𝑥𝑉₃̇𝑥,𝑡²+𝑊̇₃𝑥,𝑡²𝑑𝑥=12𝐮̇^𝑇3𝑉𝐌₃𝐮̇₃,(2.4)₃=12^𝐿0𝐸𝐽𝑥𝑊^3𝑥,𝑡²𝑑𝑥+12^𝐿0+1𝐸𝐴𝑥𝑉^3𝑥,𝑡²𝑑𝑥2^𝐿0𝐸𝐴𝑥𝑉^3𝑥,𝑡+12𝑊^3=1𝑥,𝑡²𝑊^3𝑥,𝑡²𝑑𝑥2𝐮̇^𝑇3𝐊₃𝐮̇₃+12^𝐿0𝐸𝐴𝑥𝑉^3𝑥,𝑡+12𝑊^3𝑥,𝑡²𝑊^3𝑥,𝑡²𝑑𝑥,(2.5) where 𝜌 is the average mass density of beam element, 𝐴(𝑥) is the cross-section area function, and 𝑉₃(𝑥,𝑡) and 𝑊₃(𝑥,𝑡) are, respectively, the longitudinal displacement and the transversal displacement of any point in the beam element (see the appendix). 𝐸 is the modulus of elasticity, 𝐽(𝑥) is the moment of inertia distribution function of the element. 𝐌₃ and 𝐊₃ are the mass matrix and stiffness matrix of the beam element, respectively, and they are 8×8 matrix (see the appendix).

Figure 4 

Beam element model of linkage mechanism.

We regard the links BC and CD as two beam elements, respectively, and regard the crank AB and DE as one element respectively as shown in Figure 1.    ③ ,    ④ ,    ⑤ ,    ⑥ ,    ⑦    and   ⑧    are the serial number of the beam element as shown in Figure 5. 𝑈₁₁, 𝑈₁₂,𝑈₁₅,𝑈₁₆, 𝑈₁₉, 𝑈₂₀, 𝑈₂₃, 𝑈₂₄, 𝑈₂₇, and 𝑢₂₈are elastic displacement, 𝑈₉,𝑈₁₃, 𝑈₁₄, 𝑈₁₇, 𝑈₂₁, 𝑈₂₂, 𝑈₂₅, 𝑈₂₉, 𝑈₃₀, and 𝑢₃₁ are elastic rotational angle, 𝑢₁₀, 𝑢₁₈, 𝑢₂₆, and 𝑢₃₂ are the curvature. The serial numbers of the controllable motor elements are   ①   and   ② .

Figure 5 

Diagram of finite element analysis of mechanism.

The global coordinate system is established, and the system elastic displacement vector in the global coordinate system is defined as 𝐮. Assume that 𝐑_𝑖  is the transformation matrix of the 𝑖th element between the element coordinates and the global coordinates, and 𝐁_𝑖  is the coordinate matrix of the 𝑖th element between the local numbering and the system numbering (see the appendix).

The linear viscous damper damping is adopted in this paper, according to the second Lagrange equation𝑑𝑑𝑡𝜕𝜕𝑢𝑇̇_𝑖−𝜕𝜕𝑢𝑇_𝑖+𝜕𝑉𝜕𝑢_𝑖=𝐏+𝐐.(2.6)

The electromechanical coupling nonlinear dynamic equation of the system can be obtained as 𝐌𝐮̈+𝐂𝐮̇+𝐊+𝐊₀𝐮=𝐏−𝐌𝐮̈_𝑟−𝐊^𝑒1+𝐊^𝑒2𝐞−𝐤₀+𝜀,(2.7) where 𝐮 is the generalized coordinates array of the system, 𝐌 and 𝐂 are, respectively, the 𝑛×𝑛 mass matrix and damping matrix of the system, 𝐊 and 𝐊^𝑒1 are the 𝑛×𝑛 stiffness matrixes in connection with the structural parameters of system, 𝐊₀ and 𝐊^𝑒2 are the 𝑛×𝑛 stiffness matrixes in connection with the electromagnetic parameters of system, 𝑘₀ is the 𝑛 orders array in connection with the electromagnetic parameters of system, 𝐏 is the external force array of the system, 𝐮̈_𝑟 is the rigid acceleration array of system in the global coordinates, 𝐞 is an array in connection with the eccentricity of rotor, 𝜀 is the nonlinear term and also is a small parameter, and 𝜀=⁸_𝑛=1𝐮^𝑇𝐆_𝑛𝐊_𝑛𝐮+12⁸_𝑛=1𝐮^𝑇𝐊_𝑛𝐮𝐆_𝑛+12_𝑔_𝑙𝐮^𝑇𝐆_𝑔𝑙𝐮𝐊_𝑔𝑙𝐮+12_𝑔_𝑙𝐆_𝑔𝑙𝐮𝐮^𝑇𝐊_𝑔𝑙𝐮,(𝑔,𝑙=2,3,4,6,7,8),(2.8)𝐆_𝑛,𝐊_𝑛,𝐆_𝑘𝑙 and 𝐊_𝑘𝑙 are the 𝑛×𝑛 matrixes in connection with the structural parameters of the linkage mechanism,𝐌=^𝑛_𝑖=1𝐁^𝑇𝑖𝐑^𝑇𝑖𝐌_𝑖𝐑_𝑖𝐁_𝑖,𝐊=𝐊^𝑒11+𝐊^𝑒21+^𝑛_𝑖=3𝐊_𝑒1𝐁^𝑇𝑖𝐑^𝑇𝑖𝐊₃𝐑_𝑖𝐁_𝑖,𝐊₀=𝐊^𝑒12+𝐊^𝑒22,𝐊_𝑒21=𝐊^𝑒11+𝐊^𝑒12,𝐊^𝑒2=𝐊^𝑒21+𝐊^𝑒22,𝐊^𝑒11=𝐁^𝑇1𝐑^𝑇1𝐊₁₁𝐑₁𝐁₁,𝐊^𝑒12=𝐁^𝑇1𝐑^𝑇1𝐊₁₂𝐑₁𝐁₁,=𝐁^𝑇2𝐑^𝑇2𝐊₂₁𝐑₂𝐁₂,𝐊^𝑒22=𝐁^𝑇2𝐑^𝑇2𝐊₂₂𝐑₂𝐁₂,𝐞=𝐁^𝑇1𝐑^𝑇1𝐞₁+𝐁^𝑇2𝐑^𝑇2𝐞₂,𝐞₁=𝑒₁₁+𝜀₀₁cos𝛽₁𝑒₁₂+𝜀₀₁sin𝛽₁00^𝑇,𝐆𝐞_2𝑛=𝑒₂₁+𝜀₀₂cos𝛽₂𝑒₂₂+𝜀₀₂sin𝛽₂00^𝑇,𝐤₀=𝐁^𝑇1𝐑^𝑇1𝐤₀₁+𝐁^𝑇2𝐑^𝑇2𝐤₀₂,=^𝑛_𝑖=3𝐁^𝑇𝑖𝐑^𝑇𝑖𝐠_𝑛0𝑥0𝑓𝑓0𝑐𝐠_𝑛_𝑖=1,(𝑖=1,2,…,8),𝐊_𝑛=^𝑛_𝑖=3𝐁^𝑇𝑖𝐑^𝑇𝑖^𝐿0𝐸𝐴𝐠^𝑇𝑛𝐠𝐊_𝑎𝑑𝑥𝐑_𝑖𝐁_𝑖,𝐠=−1𝐿000𝐿1000^𝑇𝐆,𝐊_𝑎_𝑖𝑗=𝛾^𝑖𝛾^𝑗,𝑖,𝑗=2,3,4,6,7,8,_𝑔𝑙=^𝑛_𝑖=3𝐊𝐁^𝑇𝑖𝐑^𝑇𝑖𝐆_𝑔𝑙𝐑_𝑖𝐁_𝑖,𝐆_𝑔𝑙_𝑔𝑙=𝐆_𝑔𝑙_𝑙𝑔=1,𝑔,𝑙=2,3,4,6,7,8,_𝑔𝑙=^𝑛_𝑖=3𝐁^𝑇𝑖𝐑^𝑇𝑖𝐊_𝑔𝑙𝐑_𝑖𝐁_𝑖,𝐊_𝑔𝑙_𝑖𝑗=𝐊_𝑔𝑙_𝑗𝑖=^𝐿0𝐸𝐴𝛾^𝑔𝛾^𝑙𝛾^𝑖𝛾^𝑗𝑑𝑥,𝑖,𝑗=2,3,4,6,7,8,(2.9)

 𝑛 is the number of the elements, and here 𝑛=8. 𝛾_𝑖 are the shape functions (see the appendix).

In (2.7), −𝐌𝐮̈_𝑟−(𝐊^𝑒1+𝐊^𝑒2)𝐞−𝐤₀ are periodic terms. They can be expanded to the form of Fourier series as follows [12]:𝐅_𝑓_𝑖=^𝑚_𝑘=1𝐹_𝑓𝑘𝑖cos𝑘𝑣_𝑓𝑡+𝜏_𝑓𝑘𝑖,(𝑓=1,2,…8),(2.10) where 𝑣₁ is the working frequency of rotor of controllable motor 1, 𝑣₂ is the working frequency of rotor of controllable motor 2, and 𝑣₃ is the smallest common multiple of 𝑣₁ and 𝑣₂. 𝑣₄=𝑣₁, 𝑣₅=𝑣₂, 𝑣₆ is the rotation frequency of magnetic field of stator of the controllable motor 1, 𝑣₇ is the rotation frequency of magnetic field of stator of the controllable motor 2, 𝑣₅ is the smallest common multiple of 𝑣₆ and 𝑣₇, and 𝑣₈ is the smallest common multiple of𝑣₁, 𝑣₂,𝑣₆, and 𝑣₇. 𝐹_𝑓𝑘𝑖 express the amplitude values, 𝜏_𝑓𝑘𝑖 are the corresponding phase angles, and 𝑚 is the number of terms of the Fourier expansion formula.

3. First Approximate Solutions

Substituting (2.10) into (2.7) and assuming that the linear transfer functions are𝐮=𝝓𝜼,(3.1) where 𝝓 is the modal transfer matrix, 𝜼 is the modal coordinate vector corresponding to 𝝓. Substituting (3.1) into (2.7), and premultiply the equation by 𝝓^𝑇, then (2.7) can be transferred as𝜼̈+𝝎²𝜼=𝐏₀+⁸_𝑓=1𝐅_0𝑓+𝜀−𝐊₀₀𝜼+⁸_𝑛=1𝜼^𝑇𝐆_0𝑛𝐊_0𝑛𝜼+12⁸_𝑛=1𝜼^𝑇𝐊_0𝑛𝜼𝐆_0𝑛+12_𝑔_𝑙+1𝜼^𝑇𝐆_0𝑔𝑙𝜼𝐊_0𝑔𝑙𝜼2_𝑔_𝑙𝐆_0𝑔𝑙𝜼𝜼^𝑇𝐊_0𝑔𝑙𝜼−𝐂₀𝜼̇,(𝑔,𝑙=2,3,4,6,7,8)(3.2) where 𝝎²=⎡⎢⎢⎢⎢⎣𝜔²¹⋮𝜔⋯0₂₂0⋯⋱𝜔⋮_2𝑛⎤⎥⎥⎥⎥⎦,𝐂₀=⎡⎢⎢⎢⎢⎣2𝜁¹⋮𝜔₁2𝜁₂⋯0𝜔₂0⋯⋱2𝜁⋮_𝑛𝜔_𝑛⎤⎥⎥⎥⎥⎦,(3.3)

𝐏₀=𝝓^𝑇𝐏, ∑^8𝑓=1𝐅_0𝑓=∑^8𝑓=1𝝓^𝑇𝐅_𝑓, 𝐊_0𝑛=𝝓^𝑇𝐊_𝑛𝝓, 𝐊₀₀=𝝓^𝑇𝐊₀𝝓, 𝐆_0𝑛=𝝓^𝑇𝐆_𝑛, 𝐌₀=𝝓^𝑇𝐌, 𝐆_0𝑔𝑙=𝝓^𝑇𝐆_𝑔𝑙𝝓, 𝐊_0𝑔𝑙=𝝓^𝑇𝐊_𝑔𝑙𝝓. 𝜁_𝑛 and 𝜔_𝑛 are the 𝑛-order damping ration mean and the 𝑛-order instantaneous natural frequency mean of the 𝑛-order canonical mode of the system in a period of motion, and 𝜁_𝑛 can be obtained through experiment.

Equation (3.2) also can be expressed as𝜂̈_𝑟+𝜔^2𝑟𝜂_𝑟=𝑃_0𝑟+⁸_𝑓=1𝐅_0𝑓𝑟+𝜀−2𝜁_𝑟𝜔_𝑟𝜂̇_𝑟−_𝑠𝛼_𝑠𝜂_𝑠+_𝑠,𝑡𝛿𝜂_𝑠𝜂_𝑡+_{𝑠,𝑡,𝑢}𝛾_𝑠𝑡𝑢𝜂_𝑠𝜂_𝑡𝜂_𝑢(𝑔,𝑙=2,3,4,6,7,8;𝑟,𝑠,𝑡,𝑢=1,2,…,𝑛),(3.4) where_𝑠𝛼_𝑠𝜂_𝑠=𝜉𝐊₀𝜼,_𝑠,𝑡𝛿_𝑠𝑡𝜂_𝑠𝜂_𝑡=⁸_𝑛=1𝜉𝜂^𝑇𝐆_𝑛𝐊_𝑛𝜼+12⁸_𝑛=1_{𝑠,𝑡,𝑢}𝜉𝜼^𝑇𝐊_𝑛𝜼𝐆_𝑛,𝜈_𝑠𝑡𝑢𝜂_𝑠𝜂_𝑡𝜂_𝑢=12_𝑔_𝑙𝜉𝜂^𝑇𝐆_𝑔𝑙𝜼𝐊_𝑔𝑙𝜼+12_𝑔_𝑙𝜉𝐆_𝑔𝑙𝜼𝜼^𝑇𝐊_𝑔𝑙𝜼,(3.5)𝛼_𝑠,𝛿_𝑠𝑡, and 𝛾_𝑠𝑡𝑢 are the coefficients of 𝜂_𝑠,𝜂_𝑠𝜂_𝑡, and 𝜂_𝑠𝜂_𝑡𝜂_𝑢, respectively. 𝜉 is the 𝑛th-order vector, and the 𝑟th element of 𝜉 is 1 and the other elements are zero. Therefore, the system is affected by multifrequency excitations.

The method of multiple time scales is employed to study the nonlinear (3.4), and the frequent factors of system can be obtained. and some resonance phenomena will take place in the system under certain conditions.(1)On condition that 𝑣_𝑝≈𝜔_𝑟, the primary resonance will take place in the system. (2)On condition that 𝜔_𝑟≈𝑘𝑣_𝑝(𝑘≠1), the ultra-harmonic resonance will take place in the system.(3) On condition that 𝜔_𝑟≈1/3𝑣_𝑝 or 𝜔_𝑟≈1/2𝑣_𝑝, the subharmonic resonance will take place in the system.(4) On condition that 𝜔_𝑠≈𝜔_𝑡+𝜔_𝑢, the inner resonance will take place in the system.(5)On condition that 𝜔_𝑟≈|±𝑗𝑣_𝑝±𝑘𝑣_𝑞| or 𝜔_𝑟≈|±2𝑗𝑣_𝑝±𝑘𝑣_𝑞|, the combination resonance will take place in the system.(6) On condition that two types of the resonance take place in the same time, the multiple resonance of the system will take place in the system.

Considering the first two orders modal of vibration, the first approximation solutions of the multi-DOF nonlinear coupling system in the generalized coordinate are also obtained by the method of multiple scales as follows:𝑢=²_𝑟=1𝜂_𝑟𝜙^(𝑟),(3.6)𝜂_𝑟=2𝑎_𝑟cos𝜔_𝑟𝑡+𝜃_𝑟+⁸_𝑚𝑓=1_𝑘=1Λ_𝑟𝑓𝑘cos𝑖𝑘𝑣_𝑓𝑡+𝜏_𝑓𝑘𝑖,(𝑟=1,2),(3.7) where 𝑎_𝑟 and 𝜃_𝑟 can be resolved by the method of Newton-Raphson, 𝜔_𝑟 is the 𝑟-order instantaneous natural frequency mean of the system, and Λ_𝑟𝑓𝑘=Γ_𝑟𝑓𝑘exp(𝑖𝜏_𝑓𝑘𝑖) and Γ_𝑟𝑓𝑘=𝐹_𝑟𝑓𝑘/(2(𝜔^2𝑟−𝑘²𝑣^2𝑓)).

Therefore, the resonance characteristics of the multi-DOF nonlinear coupling system can be analyzed using the method of multiple scales, and the first approximate solutions of the dynamic model of the system also can be obtained by that method.

4. Newmark Stable Numerical Solution

Because the implicit Newmark integration method can calculate the unconditioned stable solution when the Newmark parameters 𝜍 and 𝜗 are specific constants [13, 14], the method is adopted to obtain the stable solution of the multi-DOF nonlinear coupling system. Therefore, the stable solutions can be obtained using the implicit Newmark integration method based on the first approximate solutions calculated by the method of multiple scales.

In order to calculate the dynamic response of the system as shown in (2.7), the period 𝑡 of motion of the system is divided into 𝑆 time steps Δ𝑡 under given trajectory, that is to say, 𝑇=𝑆⋅Δ𝑡.(4.1) And the first approximate solution of dynamic response of the system is solved by the multiple-scales method firstly as shown in (3.7). The first approximate solution when 𝑡 is at a moment is looked as the initial value,and then the dynamic response of system is solved by immediate integration implicit Newmark method.

When we calculate the response of the system by implicit Newmark method, the hypothesis about the response of some moment must be taken as𝐮̇_𝑡+Δ𝑡=𝐮̇_𝑡+(1−𝜍)𝐮̈_𝑡+𝜍𝐮̈_𝑡+Δ𝑡Δ𝑡,(4.2)𝐮_𝑡+Δ𝑡=𝐮_𝑡+𝐮̇_𝑡Δ𝑡+12−𝜗𝐮̈^𝑡+𝜗𝐮_𝑡+Δ𝑡Δ𝑡²,(4.3) where 𝜍 and 𝜗 are the Newmark parameters, namely, the constant of integration and 𝜍≥0.5 and 𝜗≥0.25(0.5+𝜍)².

Because the multi-DOF nonlinear coupling dynamic (2.7) is time variant, the equation at the moment of 𝑡+Δ𝑡 can be expressed as𝐌_𝑡+Δ𝑡𝐮̈_𝑡+Δ𝑡+𝐂_𝑡+Δ𝑡𝐮̇_𝑡+Δ𝑡+𝐊_𝑡+Δ𝑡+𝐊₀_𝑡+Δ𝑡𝐮_𝑡+Δ𝑡=𝐏_𝑡+Δ𝑡−𝐌_𝑡+Δ𝑡𝐮̈_𝑟−𝐊^𝑒1+𝐊^𝑒2_𝑡+Δ𝑡𝐞−𝐤₀_𝑡+Δ𝑡+(𝜀)_𝑡+Δ𝑡,(𝑘,𝑙=2,3,4,6,7,8).(4.4)

Substituting 𝐮̈_𝑡+Δ𝑡 calculated according to (4.3) into (4.4), the 𝐮_𝑡+Δ𝑡 expressed by 𝐮_𝑡, 𝐮̇_𝑡 and 𝐮̈_𝑡 can be obtained. The calculation accuracy is related to time stepΔ𝑡.

5. Simulation and Experimental Results

Every link of the linkage mechanism is homogeneous. The width and thickness are 30 mm and 2 mm, respectively. The lengths of each links are, crank 𝐿₁=200mm and 𝐿₄=150mm, coupler 𝐿₂=𝐿₃=400mm, and frame 𝐿₅=400mm as shown in Figure 1. The material of link is aluminium. The density of aluminium 𝜌=2700kg/m³ and the Young’s modulus 𝐸=70GPa. The lumped mass of the intersection between the crank and the coupler is 𝑚₀₁=0.142kg. The lumped mass of the intersection between the two couplers is 𝑚₀₂=0.092kg. The controllable motor 1 is 90ZYT motor and motor 2 is YS8024 motor. The motors are custom made by motor manufacturer. The parameters of motors are offered by the manufacturer.

The parameters of the motors are as follows:

①   The Parameters of the Control Motor 1. The rated power of the 90ZYT motor, 𝑃_𝑁=0.75kW, the rated voltage, 𝑈_𝑘𝑛=220V, the rated rotational speed of the motor is 1500 r/min, and the stall torque is 2.0N⋅m. The static-geometric eccentricity of the motor, 𝑒₀₁=0.75mm and the rotational eccentricity, 𝜀₀₁=0.5mm. The magnetic permeability coefficient of air is 𝜇₀=4𝜋×10⁻⁷H/m, the length of the even air is gap 𝛿₀=0.25mm, and the saturation is 𝑘_𝜇=1.2. The number of excitation windings of the motor is 𝑊=924, and the coefficient is 𝐾_𝑤=0.92. The peak value of field current, 𝐼=3.58A. The number of magnetic pole-pairs of the motor is 𝑝=1.𝑚₁=2 and 𝑚₂=2 are the number of phases of the stator and rotor, respectively. The reactance of field windings is 𝑥_𝑚=594.35Ω. The reduction value of resistance and equivalent self-induction reactance of rotor, respectively, are 𝑟^=27.24Ω and 𝑥^=0.0196Ω. The slide ratio is 𝑠=0.0713. The control voltage is 𝑈_𝑘=26V. The mass of the motor rotor is 𝑚₀=2.2kg. The moment of inertia of the motor rotor is 𝐽₀=0.018kg⋅m². The length of the motor shaft is 𝑙=363mm(𝑙₁=130mm,𝑙₁=53mm, and𝑙₁=180mm), the effective length of the rotor 𝐿₀₁=140mm, and the inner radius of the motor stator is 𝑅₀₁=23mm.

②   The Parameters of the Control Motor 2. The rated power of the YS8024 motor is 𝑃_𝑁=0.75Kw, the rated current is 𝐼_𝑁=3.48/2.01A, the rated voltage is 𝑈_𝑘𝑛=220V, and the rated rotational speed is 𝑛_N=1440r/min. The static-geometric eccentricity is 𝑒₀₂=0.73mm and the rotational eccentricity is 𝜀₀₂=0.5mm. The magnetic permeability coefficient of air is 𝜇₀=4𝜋×10⁻⁷H/m the length of the even air gap is 𝛿₀=0.25mm, and the saturation is 𝑘_𝜇=1.2. The number of excitation windings of the motor is 𝑊=824 and the coefficient is 𝐾_𝑤=1. The peak value of field current is 𝐼=3.58A. The number of the magnetic pole-pair of compounded magnetic field is 𝑝=2. 𝑚₁=3 and 𝑚₂=0.5 are the number of phases of the stator and rotor respectively. The reactance of field windings is 𝑥_𝑚=600Ω. The reduction value of resistance and equivalent self-induction reactance of rotor, respectively, are 𝑟^=30Ω and 𝑥^=0.02Ω. The slide ratio is 𝑠=0.15. The control voltage is 𝑈_𝑘=95V. The mass of the motor rotor is 𝑚₀=2.93kg. The moment of inertia of the motor rotor 𝐽₀=0.021kg⋅m². The length of the motor shaft 𝑙=208mm(𝑙₁=100mm,𝑙₂=33mm, and 𝑙₃=175mm), the effective length of the rotor 𝐿₀₁=100mmand the inner radius of the motor stator is 𝑅₀₁=38mm.

As shown in Figure 1, The initial angles of the two cranks are 0°. The calculation initial value of the response is the one when the two cranks move to the initial position of them after the system comes to the stabilized state. Given 𝜍=0.25 and 𝜗=0.5, the period of motion of the linkage mechanism is divided into S=100 time steps.

The responses of transversal and longitudinal displacement of the midpoints of the links 𝐿₂, and 𝐿₃ can be calculated by the method of multiple-scales Newmark mentioned above, then the dynamic responses on the direction perpendicular to the links’ axis can be simulated by the Multiple-Scales-Newmark method mentioned above as shown in Figure 6, and the corresponding experimental curve of dynamic response of the links’ midpoint on the same parameter conditions can be obtained through experiment as shown in Figure 6. The dynamic responses of midpoints of the links are measured by the dynamic test system in experiment. Comparing the simulation with experimental figures, one can find that the multiple-scales Newmark method studied in this paper is correct and practicable.

(a)

(b)

Figure 6 

Simulation and experiment curve of dynamic response for midpoint of link 𝐿2 and 𝐿3.

6. Conclusions

The numerical method can calculate the dynamic responses, but that method could not analyze internal relationship between the dynamic characteristics and the scale and electromagnetism parameters. Though the multiple-scales method can analyze the dynamics mechanism, it usually adopt first approximations and second, or higher-order approximation is very complex. An evolutionary analytic method of multiple-scales Newmark is firstly presented in this paper to study the dynamic characteristics of such kind of multi-DOF nonlinear coupling system synthetically using the property that the implicit Newmark method can calculate the unconditioned stable solution when the parameters of Newmark are specific constants with the features of the multiple scales method. The numerical simulation and experimental results indicate that the multiple-scales Newmark method mentioned in this paper is correct and practicable. The study provides the basis of further study on the dynamic control of such kind of mechanism.

Appendix

 

According to the finite element method, 𝑊₁(𝑥,𝑡) and 𝑉₁(𝑥,𝑡) can be expressed as𝑊𝑉₁(𝑥,𝑡)=_𝑖𝜓_𝑖(𝑥)𝑈_𝑖(𝑡),𝑖=1,2,₁(𝑥,𝑡)=_𝑖𝜓_𝑖(𝑥)𝑈_𝑖(𝑡),𝑖=3,4,(A.1)

where 𝜓_𝑖(𝑥) are the shape functions, 𝜓_1,2(𝑥)=⎧⎪⎪⎨⎪⎪⎩1−10𝑒³¹1−10𝑒₃₂+15𝑒⁴¹−6𝑒⁵¹,𝑥≤𝑙₁,0,+15𝑒⁴²−6𝑒⁵²,𝑙𝑙₁₁<𝑥≤𝑙₁+𝑙₂,𝜓₃+𝑙₂<𝑥≤𝑙₁+𝑙₂+𝑙₃,(𝑥)=⎧⎪⎪⎨⎪⎪⎩0,𝑥≤𝑙¹1−𝑒₂,𝑒₃,𝑙₁<𝑥≤𝑙₁+𝑙₂,𝜓₄,𝑙₁+𝑙₂<𝑥≤𝑙₁+𝑙₂+𝑙₃,(𝑥)=⎧⎪⎪⎨⎪⎪⎩0,𝑥≤𝑙¹0,𝑙₁,1−𝑒₃<𝑥≤𝑙₁+𝑙₂,,𝑙₁+𝑙₂<𝑥≤𝑙₁+𝑙₂+𝑙₃,(A.2)where 𝑒₁=𝑥/𝑙₁, 𝑒₂=(𝑙₁+𝑙₂−𝑥)/𝑙₂, 𝑒₃=(𝑙₁+𝑙₂+𝑙₃−𝑥)/𝑙₃, 𝑙₁, 𝑙₂, and 𝑙₃ are, respectively, the length between point 1 and point 2, point 2 and point 3, point 3 and point 4. 𝑈_𝑖(𝑡) is the nodal displacement as shown in Figure 2.

The mass matrix 𝐌₁ of the electromotor element is 4×4 matrix, 𝐌₁₁₁=^𝑙0𝑚₁(𝑥)𝜓₁(𝑥)𝜓₁(𝑥)𝑑𝑥,𝐌₁₂₂=^𝑙0𝑚₁(𝑥)𝜓₂(𝑥)𝜓₂(𝑥)𝑑𝑥,𝐌₁_𝑘𝑝=^𝑙0𝑚₁(𝑥)𝜓_𝑘(𝑥)𝜓_𝑝(𝑥)𝑑𝑥,(𝑘,𝑝=3,4),(A.3)

and the other components are equal to zeros.

The stiffness matrixes of the electromotor element are 𝐊_𝟏𝟏₁₁=^𝑙0𝐸₁𝐼₁(𝑥)𝜕²𝜒₁(𝑥)𝜕𝑥²²𝑑𝑥,𝐊_𝟏𝟏₁₂=𝐊₁₁₂₁=^𝑙0𝐸₁𝐼₁(𝑥)𝜕𝜒¹(𝑥)𝜕𝑥𝜕𝜒²(𝑥)𝜕𝑥𝑑𝑥,𝐊_𝟏𝟏₂₂=^𝑙0𝐸₁𝐼₁(𝑥)𝜕²𝜒₂(𝑥)𝜕𝑥²²𝑑𝑥,𝐊_𝟏𝟏_𝑘𝑝=^𝑙0𝐺₁𝐽₁(𝑥)𝜕𝜒^𝑘(𝑥)𝜕𝑥⋅𝜕𝜒^𝑝(𝑥)𝜕𝑥𝑑𝑥,𝑘,𝑝=3,4,𝐊_𝟏𝟐₁₁=𝑝𝑅⁰¹𝑙₀₁Λ₀₁4𝜎²²^02𝜋(1+cos2𝛼)𝐹^+𝑠+𝐹_+𝑟cos𝜔₀₁𝑡−𝑝𝛼+𝐹_−𝑠cos𝜔₀₁𝑡+𝑝𝛼cos𝜔₀₁𝑡−𝑝𝛼−𝜑₁₀+𝐹_−𝑟cos𝜔₀₁𝑡+𝑝𝛼−𝜑₂₀²𝑑𝛼,𝐊_𝟏𝟐₁₂=𝐊_𝟏𝟐₂₁=𝑝𝑅⁰¹𝑙₀₁Λ₀₁4𝜎²²^02𝜋(sin2𝛼)𝐹_+𝑠+𝐹_+𝑟cos𝜔₀₁𝑡−𝑝𝛼+𝐹_−𝑠cos𝜔₀₁𝑡+𝑝𝛼cos𝜔₀₁𝑡−𝑝𝛼−𝜑₁₀+𝐹_−𝑟cos𝜔₀₁𝑡+𝑝𝛼−𝜑₂₀²𝑑𝛼,𝐊_𝟏𝟐₂₂=𝑝𝑅⁰¹𝑙₀₁Λ₀₁4𝜎²²^02𝜋(1−cos2𝛼)𝐹_+𝑠+𝐹_+𝑟cos𝜔₀₁𝑡−𝑝𝛼+𝐹_−𝑠cos𝜔₀₁𝑡+𝑝𝛼cos𝜔₀₁𝑡−𝑝𝛼−𝜑₁₀+𝐹_−𝑟cos𝜔₀₁𝑡+𝑝𝛼−𝜑₂₀²𝑑𝛼,𝐊_𝟏𝟐_𝑘𝑝=0,𝑘,𝑝=3,4,𝐊_𝟏𝟏₁₃=𝐊_𝟏𝟏₁₄=𝐊_𝟏𝟏₂₃=𝐊_𝟏𝟏₂₄=𝐊_𝟏𝟏₃₁=𝐊_𝟏𝟏₃₂=𝐊_𝟏𝟏₄₁=𝐊_𝟏𝟏₄₂=𝐊_𝟏𝟐₁₃=𝐊_𝟏𝟐₁₄=𝐊_𝟏𝟐₂₃=𝐊_𝟏𝟐₂₄=𝐊_𝟏𝟐₃₁=𝐊_𝟏𝟐₃₂=𝐊_𝟏𝟐₄₁=𝐊_𝟏𝟐₄₂=0,𝐤₀₁₁=𝑝𝑅⁰¹𝑙₀₁Λ₀₁2𝜎^02𝜋(cos𝛼)𝐹^+𝑠+𝐹_+𝑟cos𝜔₀₁𝑡−𝑝𝛼+𝐹_−𝑠cos𝜔₀₁𝑡+𝑝𝛼cos𝜔₀₁𝑡−𝑝𝛼−𝜑₁₀+𝐹_−𝑟cos𝜔₀₁𝑡+𝑝𝛼−𝜑₂₀²𝑑𝛼,𝐤₀₁₂=𝑝𝑅⁰¹𝑙₀₁Λ₀₁2𝜎^02𝜋(sin𝛼)𝐹^+𝑠+𝐹_+𝑟cos𝜔₀₁𝑡−𝑝𝛼+𝐹_−𝑠cos𝜔₀₁𝑡+𝑝𝛼cos𝜔₀₁𝑡−𝑝𝛼−𝜑₁₀+𝐹_−𝑟cos𝜔₀₁𝑡+𝑝𝛼−𝜑₂₀²𝑑𝛼,𝐤₀₁₃=𝐤₀₁₄=0.(A.4)

According to the finite element method, 𝑉₃(𝑥,𝑡) and 𝑊₃(𝑥,𝑡) can be expressed as𝑊𝑉₃𝑥,𝑡=_𝑖𝛾_𝑖𝑥𝑢_𝑖(𝑡),(𝑖=1,5),₃𝑥,𝑡=_𝑖𝛾_𝑖𝑥𝑢_𝑖(𝑡),(𝑖=2,3,4,6,7,8),(A.5)

where 𝑥 is the coordinate of beam element in local coordinate system, and the shape functions are as follows: 𝛾𝛾₁𝑥=1−𝑒,₂𝛾₃𝑥=1−10𝑒³+15𝑒⁴−6𝑒⁵,𝛾₄𝑥=𝐿×𝑒−6𝑒³+8𝑒⁴−3𝑒⁵,𝑥=𝐿²×𝑒²−3𝑒³+3𝑒⁴−5𝑒⁵𝛾₅2,𝛾𝑥=𝑒,₆𝛾₇𝑥=10𝑒³−15𝑒⁴+6𝑒⁵,𝛾𝑥₈=𝐿×−4𝑒³+7𝑒⁴−3𝑒⁵,𝑥=𝐿²×𝑒³−2𝑒⁴+𝑒⁵2,(A.6) where 𝑒=𝑥/𝐿, 𝐿 is the length of the beam element.

The mass matrix and stiffness matrix of the beam element are as follows:𝐌_𝟑_𝑖𝑗=𝜌𝐴𝑥^𝐿0𝛾_𝑖𝑥𝛾_𝑗𝑥𝑑𝑥+𝑝𝐼^𝐿0𝛾^𝑖𝑥𝛾^𝑗𝑥𝑑𝑥,(𝑖,𝑗=2,3,4,6,7,8),𝐌_𝟑_𝑖𝑗=𝜌𝐴𝑥^𝐿0𝛾_𝑖𝑥𝛾_𝑗𝑥𝑑𝑥,(𝑖,𝑗=1,5),(A.7)𝐊₃_𝑖𝑗=𝐸𝐽𝑥^𝐿0𝛾^𝑖𝑥𝛾^𝑗𝑥𝑑𝑥,(𝑖,𝑗=2,3,4,6,7,8),𝐊₃_𝑘𝑝=𝐸𝐴𝑥^𝐿0𝛾^𝑘𝑥𝛾^𝑝𝑥𝑑𝑥,(𝑘,𝑝=1,5),(A.8)

where 𝐼 is the moment of inertia of cross-section.

The positive-sequence and negative-sequence components of the magnetomotive amplitude of stator and rotor, respectively, are𝐹_+𝑠=121+𝑎^e𝐹_𝑚𝐹_−𝑠,=121−𝑎^e𝐹_𝑚𝐹,(A.9)_+𝑟=12𝑚𝑚²₁𝑥_𝑚√(𝑟^/𝑠)²+(𝑥^)²𝐹_−𝑟1+𝑎₂𝐹_𝑚,=12𝑚𝑚²₁𝑥_𝑚√[𝑟^/(2−𝑠)]²+(𝑥^)²1+𝑎₂𝐹_𝑚,(A.10)

where 𝐹_𝑚=0.9(𝑊𝐾_𝑤)/𝑝𝐼, 𝑊 is the number of turns of field winding, 𝐾_𝑤 is the coefficient of field winding, 𝑝 is the number of pole pairs, 𝐼 is the peak value of field current. 𝑎_e=𝑈_𝑘/𝑈_𝑘𝑛 is the effective signal coefficient, 𝑈_𝑘𝑛 is the rated voltage, and 𝑈_𝑘 is the actual control voltage. 𝑚₁ and 𝑚₂ are the number of phases of the stator and rotor, respectively, 𝑥_𝑚 is the reactance of field windings, 𝑟′ and 𝑥′ are respectively the reduction value of resistance and equivalent self-induction reactance of rotor windings, and the slip ratio 𝑠 is the difference between the rotor speed and forward revolving field. 2−𝑠 is the slip ratio between the rotor speed and backward revolving field.

The transformation matrix of the 𝑖th element 𝐑_𝑖 is as follows

  𝐑𝐑₁=𝐑₂=⎡⎢⎢⎢⎢⎢⎣101101₃⎤⎥⎥⎥⎥⎥⎦,=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣cos𝛽¹−sin𝛽₁−sin𝛽₁cos1𝛽₁0(4×4)1cos𝛽₁−sin𝛽₁−sin𝛽₁cos𝛽₁0(4×4)1𝐑₄1⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,=𝐑₅=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣cos𝛽²−sin𝛽₂−sin𝛽₂cos1𝛽₂0(4×4)1cos𝛽₂−sin𝛽₂−sin𝛽₂cos𝛽₂0(4×4)1𝐑₆1⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,=𝐑₇=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣cos𝛽³−sin𝛽₃−sin𝛽₃cos1𝛽₃0(4×4)1cos𝛽₃−sin𝛽₃−sin𝛽₃cos𝛽₃0(4×4)1𝐑₈1⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣cos𝛽⁴−sin𝛽₄−sin𝛽₄cos1𝛽₄0(4×4)1cos𝛽₄−sin𝛽₄−sin𝛽₄cos𝛽₄0(4×4)11⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,(A.11)

where 𝛽_𝑖 is the position angle of the link 𝐿_𝑖 with the level position.

The coordinate matrix of the 𝑖th element 𝐁_𝟏 and 𝐁_𝟐 is 4×32 matrixes, and 𝐁₃, 𝐁₄, 𝐁₅, 𝐁₆, 𝐁₇, and 𝐁₈ are 8×24 matrixes, and𝐁𝐁₁₁₁=𝐁₁₂₂=𝐁₁₃₃=𝐁₁₄₄=1,₂𝐁₁₅=𝐁₂₂₆=𝐁₂₃₇=𝐁₂₄₈=1,₃𝐁_3,9=𝐁₃_4,10=𝐁₃_5,11=𝐁₃_6,12=𝐁₃_7,13=1,₄𝐁_1,11=𝐁₄_2,12=𝐁₄_3,14=𝐁₄_5,15=𝐁₄_6,16=𝐁₄_7,17=𝐁₄_8,18=1,₅𝐁_1,15=𝐁₅_2,16=𝐁₅_3,17=𝐁₅_4,18=𝐁₅_5,19=𝐁₅_6,20=𝐁₅_7,21=1,₆𝐁_1,19=𝐁₆_2,20=𝐁₆_3,21=𝐁₆_3,22=𝐁₆_5,23=𝐁₆_6,24=𝐁₆_7,25=𝐁₆_8,26=1,₇_1,23=𝐁𝐁₇_2,24=𝐁₇_3,25=𝐁₇_4,26=𝐁₇_5,27=𝐁₇_6,28=𝐁₇_7,29=1,₈_1,27=𝐁₈_2,28=𝐁₈_3,30=𝐁₈_7,31=𝐁₈_8,32=1,(A.12)

and the other components of 𝐁_𝑖 are 0.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (Grant no. 51075077), the Scientific Research Foundation of Guangxi University (Grant no. X081018), and the Teams for Innovation in the Construction of Talent Highlands in Guangxi Institutions of Higher Learning. The authors gratefully acknowledge the supporting agencies.