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

Authors: Ru-Gui Wang; Gan-Wei Cai; Xiao-Rong  Zhou
Journal: Mathematical Problems in Engineering
 (2011)

---

## 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 oft 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.


---

## Body

## 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 Figure1. 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 u1=[U1U2U3U4]T, and XYZ 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 Figure3, where point O is the inner circle geometric center of the motor stator, point O1 is the outer circle geometric center of the rotor journal, point O2 is the outer circle center of the journal under the deformation of shaft or bearing, the coordinate of point O3(x,y) is the outer circle geometric center of the rotor, δ is the length of air gap, and e1 is the air-gap eccentricity. e01=e112+e122 is the static eccentricity, which is caused by the rotor gravity and mismachining tolerance of motor. ε01 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, e1=x2+y2,(2.1)x=U1+e11+ε01cos⁡β1,y=U2+e12+ε01sinβ1,
where u1 and u2 are the components of vibration eccentricity in the x and y direction, and e0=u12+u22. β1=(1-s1)ω01t is the rotational angle of rotor with respect to the stator of controllable motor. ω01 is the synchronous speed of rotation of controllable motor. s1is the slide ratio.Diagram of air-gap eccentric vibration of controllable motor.
(a)(b)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:(2.2)T1=12∫0lm1(x)[Ẇ1(x,t)]2dx+12∫0lJ01(x)[V̇1(x,t)]2dx=12u̇1TM̅1u̇1,(2.3)V1=12u1T(K̅11+K̅12)u1+e̅1T(K̅11+K̅12)u1+12e̅1T(K̅11+K̅12)e̅1+(e̅1+u1)k̅01+pR01L01Λ012∫02π[F+scos⁡(ω01t-pα)+F-scos⁡(ω01t+pα)+F+rcos⁡⁡(ω01t-pα-φ10)+F-rcos⁡(ω01t+pα-φ20)]2dα,
where W1(x,t)and V1(x,t) the transverse displacements and the angle of elastic torsion of any points in the controllable motor element, respectively (see the appendix), and l=l1+l2+l3 is the length of controllable motor shaft (see the appendix);m1(x), including the rotor mass m0 which is at x=l1, is the mass distribution function of the controllable motor shaft;J01(x), including the moment of inertia of rotorJ0which is atx=l1, is the moment of inertia distribution function of the controllable motor shaft. R01 is the inner radius of the controllable motor stator, L01 is the effective length of the rotor, Λ01=μ0/σ is the even air gap permeance of the controllable motor, μ0 is the magnetic permeability coefficient of air, σ=kμδ0, kμ is saturation, δ0 is the uniform air-gap size, kμ=1+δFe/k1δ0, k1 is the calculation air-gap coefficient of the even air gap, δFe is the equivalent air-gap of ferromagnetic materials, φ10 is the phase angle of the positive-sequence current of rotor lagging behind the positive-sequence current of stator, and φ20 is the phase angle of the negative-sequence current of rotor lagging behind the negative-sequence current of stator. F+s,F-s,F+r, and F-r are the positive-sequence and negative-sequence components of the magnetomotive amplitude of stator and rotor respectively (see the appendix). M̅1 is the mass 4×4 matrix of controllable motor element, K̅11 is the stiffness 4×4 matrix of controllable motor element in connection with the structural parameters of the rotor, K̅12 is the stiffness 4×4 matrix of controllable motor element in connection with the electromagnetic parameters of the rotor, and k̅01 is the 4-order vector in connection with the eccentric motor (see the appendix). e̅=[e11+ε01cosβ1e12+ε01sinβ100]T 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 Figure4. 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:
(2.4)T3=12∫0LρA(x̅)([V̇3(x̅,t)]2+[Ẇ3(x̅,t)]2)dx̅=12u̇3TM̅3u̇3,(2.5)V3=12∫0LEJ(x̅)[W3′′(x̅,t)]2dx̅+12∫0LEA(x̅)[V3′(x̅,t)]2dx̅+12∫0L{EA(x̅)[V3′(x̅,t)+12[W3′(x̅,t)]2}[W3′(x̅,t)]2dx̅=12u̇3TK̅3u̇3+12∫0L{EA(x̅)[V3′(x̅,t)+12[W3′(x̅,t)]2}[W3′(x̅,t)]2dx̅,
where ρ is the average mass density of beam element, A(x̅) is the cross-section area function, and V3(x̅,t) and W3(x̅,t) are, respectively, the longitudinal displacement and the transversal displacement of any point in the beam element (see the appendix). E is the modulus of elasticity, J(x̅) is the moment of inertia distribution function of the element. M̅3 and K̅3 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 Figure1.    ③ ,    ④ ,    ⑤ ,    ⑥ ,    ⑦    and   ⑧    are the serial number of the beam element as shown in Figure 5. U11, U12,U15,U16, U19, U20, U23, U24, U27, and u28are elastic displacement, U9,U13, U14, U17, U21, U22, U25, U29, U30, and u31 are elastic rotational angle, u10, u18, u26, and u32 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 asu. Assume that Ri  is the transformation matrix of the ith element between the element coordinates and the global coordinates, and Bi  is the coordinate matrix of the ith 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)ddt(∂T∂u̇i)-∂T∂ui+∂V∂ui=P+Q.The electromechanical coupling nonlinear dynamic equation of the system can be obtained as(2.7)Mü+Cu̇+(K+K0)u=P-Mür-(K1e+K2e)e-k0+ε,
where u is the generalized coordinates array of the system, M and C are, respectively, the n×n mass matrix and damping matrix of the system, K and K1e are the n×n stiffness matrixes in connection with the structural parameters of system, K0 and K2e are the n×n stiffness matrixes in connection with the electromagnetic parameters of system, k0 is the n orders array in connection with the electromagnetic parameters of system, P is the external force array of the system, ür is the rigid acceleration array of system in the global coordinates, e is an array in connection with the eccentricity of rotor, ε is the nonlinear term and also is a small parameter, and (2.8)ε=∑n=18uTGnKnu+12∑n=18uTKnuGn+12∑g∑luTGgluKglu+12∑g∑lGgluuTKglu,(g,l=2,3,4,6,7,8),Gn,Kn,Gkl and Kkl are the n×n matrixes in connection with the structural parameters of the linkage mechanism,(2.9)M=∑i=1nBiTRiTM̅iRiBi,K=K11e+K21e+∑i=3nBiTRiTK̅3RiBi,K0=K12e+K22e,K1e=K11e+K12e,K2e=K21e+K22e,K11e=B1TR1TK̅11R1B1,K12e=B1TR1TK̅12R1B1,K21e=B2TR2TK̅21R2B2,K22e=B2TR2TK̅22R2B2,e=B1TR1Te̅1+B2TR2Te̅2,e̅1=[e11+ε01cos⁡β1e12+ε01sinβ100]T,e̅2=[e21+ε02cos⁡β2e22+ε02sinβ200]T,k0=B1TR1Tk̅01+B2TR2Tk̅02,Gn=∑i=3nBiTRiTg̅n，(g̅n)i=1,(i=1,2,…,8),Kn=∑i=3nBiTRiT[∫0LEAg̅nTg̅K̅adx̅]RiBi,g̅=[-1L0001L000]T,(K̅a)ij=γi′γj′,i,j=2,3,4,6,7,8,Ggl=∑i=3nBiTRiTG̅glRiBi,(G̅gl)gl=(G̅gl)lg=1,g,l=2,3,4,6,7,8,Kgl=∑i=3nBiTRiTK̅glRiBi,(K̅gl)ij=(K̅gl)ji=∫0LEAγg′γl′γi′γj′dx̅,i,j=2,3,4,6,7,8,n is the number of the elements, and here n=8. γi are the shape functions (see the appendix).In (2.7), -Mür-(K1e+K2e)e-k0 are periodic terms. They can be expanded to the form of Fourier series as follows [12]:(2.10)(Ff)i=∑k=1mFfkicos⁡(kvft+τfki),(f=1,2,…8),
where v1 is the working frequency of rotor of controllable motor 1, v2 is the working frequency of rotor of controllable motor 2, and v3 is the smallest common multiple of v1 and v2. v4=v1, v5=v2, v6 is the rotation frequency of magnetic field of stator of the controllable motor 1, v7 is the rotation frequency of magnetic field of stator of the controllable motor 2, v5 is the smallest common multiple of v6 and v7, and v8 is the smallest common multiple ofv1, v2,v6, and v7. Ffki express the amplitude values, τfki are the corresponding phase angles, and m 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)u=ϕη,
where ϕ is the modal transfer matrix, η is the modal coordinate vector corresponding to ϕ. Substituting (3.1) into (2.7), and premultiply the equation by ϕT, then (2.7) can be transferred as(3.2)η̈+ω2η=P0+∑f=18F0f+ε[12∑g∑lG0glηηTK0glη-K00η+∑n=18ηTG0nK0nη+12∑n=18ηTK0nηG0n+12∑g∑lηTG0glηK0glη+12∑g∑lG0glηηTK0glη-C0η̇],(g,l=2,3,4,6,7,8)
where
(3.3)ω2=[ω12⋯0⋮ω22⋱⋮0⋯ωn2],C0=[2ζ1ω1⋯0⋮2ζ2ω2⋱⋮0⋯2ζnωn],P0=ϕTP, ∑f=18F0f=∑f=18ϕTFf, K0n=ϕTKnϕ, K00=ϕTK0ϕ, G0n=ϕTGn, M0=ϕTM, G0gl=ϕTGglϕ, K0gl=ϕTKglϕ. ζn and ωn are the n-order damping ration mean and the n-order instantaneous natural frequency mean of the n-order canonical mode of the system in a period of motion, and ζn can be obtained through experiment.Equation (3.2) also can be expressed as(3.4)η̈r+ωr2ηr=P0r+∑f=18F0fr+ε(-2ζrωrη̇r-∑sαsηs+∑s,tδηsηt+∑s,t,uγstuηsηtηu)(g,l=2,3,4,6,7,8;r,s,t,u=1,2,…,n),
where(3.5)∑sαsηs=ξK0η,∑s,tδstηsηt=∑n=18ξηTGnKnη+12∑n=18ξηTKnηGn,∑s,t,uνstuηsηtηu=12∑g∑lξηTGglηKglη+12∑g∑lξGglηηTKglη,αs,δst, and γstu are the coefficients of ηs,ηsηt, and ηsηtηu, respectively. ξ is the nth-order vector, and the rth 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 thatvp≈ωr, the primary resonance will take place in the system.(2)
On condition thatωr≈kvp(k≠1), the ultra-harmonic resonance will take place in the system.(3)
On condition thatωr≈1/3vp or ωr≈1/2vp, the subharmonic resonance will take place in the system.(4)
On condition thatωs≈ωt+ωu, the inner resonance will take place in the system.(5)
On condition thatωr≈|±jvp±kvq| or ωr≈|±2jvp±kvq|, 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:(3.6)u=∑r=12ηrϕ(r),(3.7)ηr=2arcos⁡(ωrt+θr)+∑f=18∑k=1m(Λrfkcos⁡(ikvft+τfki)),(r=1,2),
where ar and θr can be resolved by the method of Newton-Raphson, ωr is the r-order instantaneous natural frequency mean of the system, and Λrfk=Γrfkexp(iτfki) and Γrfk=Frfk/(2(ωr2-k2vf2)).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 t of motion of the system is divided into S time steps Δt under given trajectory, that is to say, (4.1)T=S⋅Δt.
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 t 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(4.2)u̇t+Δt=u̇t+((1-ς)üt+ςüt+Δt)Δt,(4.3)ut+Δt=ut+u̇tΔt+((12-ϑ)üt+ϑut+Δt)Δt2,
where ς and ϑ are the Newmark parameters, namely, the constant of integration and ς≥0.5 and ϑ≥0.25(0.5+ς)2.Because the multi-DOF nonlinear coupling dynamic (2.7) is time variant, the equation at the moment of t+Δt can be expressed as(4.4)Mt+Δtüt+Δt+Ct+Δtu̇t+Δt+(Kt+Δt+(K0)t+Δt)ut+Δt=Pt+Δt-Mt+Δtür-(K1e+K2e)t+Δte-(k0)t+Δt+(ε)t+Δt,(k,l=2,3,4,6,7,8).Substitutingüt+Δt calculated according to (4.3) into (4.4), the ut+Δt expressed by ut, u̇t and üt can be obtained. The calculation accuracy is related to time stepΔt.


## 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, crankL1=200mm and L4=150mm, coupler L2=L3=400mm, and frame L5=400mm as shown in Figure 1. The material of link is aluminium. The density of aluminium ρ=2700kg/m3 and the Young’s modulus E=70GPa. The lumped mass of the intersection between the crank and the coupler is m01=0.142kg. The lumped mass of the intersection between the two couplers is m02=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, PN=0.75kW, the rated voltage, Ukn=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, e01=0.75mm and the rotational eccentricity, ε01=0.5mm. The magnetic permeability coefficient of air is μ0=4π×10-7H/m, the length of the even air is gap δ0=0.25mm, and the saturation is kμ=1.2. The number of excitation windings of the motor is W=924, and the coefficient is Kw=0.92. The peak value of field current, I=3.58A. The number of magnetic pole-pairs of the motor is p=1.m1=2 and m2=2 are the number of phases of the stator and rotor, respectively. The reactance of field windings is xm=594.35Ω. The reduction value of resistance and equivalent self-induction reactance of rotor, respectively, are r′=27.24Ω and x′=0.0196Ω. The slide ratio is s=0.0713. The control voltage is Uk=26V. The mass of the motor rotor is m0=2.2kg. The moment of inertia of the motor rotor is J0=0.018kg·m2. The length of the motor shaft is l=363mm(l1=130mm,l1=53mm, andl1=180mm), the effective length of the rotor L01=140mm, and the inner radius of the motor stator is R01=23mm.②   The Parameters of the Control Motor 2.The rated power of the YS8024 motor is PN=0.75Kw, the rated current is IN=3.48/2.01A, the rated voltage is Ukn=220V, and the rated rotational speed is nN=1440r/min. The static-geometric eccentricity is e02=0.73mm and the rotational eccentricity is ε02=0.5mm. The magnetic permeability coefficient of air is μ0=4π×10-7H/m the length of the even air gap is δ0=0.25mm, and the saturation is kμ=1.2. The number of excitation windings of the motor is W=824 and the coefficient is Kw=1. The peak value of field current is I=3.58A. The number of the magnetic pole-pair of compounded magnetic field is p=2. m1=3 and m2=0.5 are the number of phases of the stator and rotor respectively. The reactance of field windings is xm=600Ω. The reduction value of resistance and equivalent self-induction reactance of rotor, respectively, are r′=30Ω and x′=0.02Ω. The slide ratio is s=0.15. The control voltage is Uk=95V. The mass of the motor rotor is m0=2.93kg. The moment of inertia of the motor rotor J0=0.021kg·m2. The length of the motor shaft l=208mm(l1=100mm,l2=33mm, and l3=175mm), the effective length of the rotor L01=100mmand the inner radius of the motor stator is R01=38mm.As shown in Figure1, 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 linksL2, and L3 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.Simulation and experiment curve of dynamic response for midpoint of linkL2 and L3.
(a)(b)


## 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.



---

*Source: 254547-2011-06-23.xml*