Next Article in Journal
A Scientometric Analysis and Visualization of Forest Soil Contamination Research from Global Perspectives
Previous Article in Journal
Comparison of Planted Pine versus Natural Mix Forests in Nepal
Previous Article in Special Issue
A Comparative Analysis of Tannin and Commercial Fire Retardants in Wood Fire Protection
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Stability Analysis of Planetary Rotor with Variable Speed Self Rotation and Uniform Eccentric Revolution in the Rubber Tapping Machinery

1
Rubber Research Institute, Chinese Academy of Tropical Agricultural Science, Haikou 570101, China
2
Mechanical Equipment Subcenter of National Center of Important Tropical Crops Engineering and Technology Research, Haikou 571101, China
*
Author to whom correspondence should be addressed.
Forests 2024, 15(6), 1071; https://doi.org/10.3390/f15061071
Submission received: 12 May 2024 / Revised: 14 June 2024 / Accepted: 17 June 2024 / Published: 20 June 2024
(This article belongs to the Special Issue Advances in the Study of Wood Mechanical and Physical Properties)

Abstract

:
Natural rubber is a critical material that is essential to industry and transportation. In order to reduce the cost of rubber tapping and improve the efficiency and profitability of rubber production, the 4GXJ-2 portable electric rubber cutter and automatic rubber tapping robot have been developed. In their vibration tool holder, the planetary rotor with variable speed self rotation and uniform eccentric revolution is the most important transmission component, and its instability will cause irregular vibration of the tapping tool, thereby reducing the accuracy of vibration cutting and increasing noise. Base on the ANCF (Absolute Nodal Coordinate Formulation) 3D-beam element and 3D REF (3D Ring on Elastic Foundation), a novel eccentric 3D REF model of a planetary rotor is proposed. By introducing multiple coordinate systems, the coupled motion of uniform eccentric revolution, variable speed self rotation and flexible deformation is decomposed and the influences of these motions on the centrifugal force and Coriolis force are more clearly derived. The model is degraded and validated by comparing with other examples of a rotating circular ring model and uniformly eccentrically revolving annular plate. According to the Floquet theory and Runge−Kutta method, the unstable region of revolution speed of a planetary rotor in rubber tapping machinery is predicted as [817 rad/s, 909 rad/s], [1017 rad/s, 1095 rad/s] and [1263 rad/s,1312 rad/s]. Compared with the rubber-tapping experiment of rubber tapping machinery, the validity of the proposed model is further verified. This model provides important design references for the speed settings of those rubber tapping machines.

1. Introduction

Compared to chemical synthetic rubber, natural rubber possesses better comprehensive properties and has become an irreplaceable and important material in many fields, such as transportation, medical treatment, and mechanical manufacturing [1,2,3]. As shown in Figure 1a, more than 98% of natural rubber is obtained by rubber tapping [4]. The bark of a rubber tree has biological characteristics of lamination and uneven thickness [5], and it is necessary to cut off as much of the outer secondary phloem layer as possible without damaging the inner secondary phloem layer and cambium layer during the rubber tapping process. Without a labor-saving structure and auxiliary locator, the traditional tools (Figure 1b) limit work efficiency and increase the labor intensity and skill requirements of rubber workers in the process of manual rubber tapping [6]. Because of the influence of temperature and air humidity, rubber tapping usually needs to be performed in the small hours [7], which further increases the difficulty of manual rubber tapping. With the increasing cost of labor and the long-term low prices of natural rubber, the research and development of mechanized rubber tapping tools have become key issues related to the natural rubber industry.
In recent years, many scholars have attempted to use intelligent rubber tapping machinery, including rubber tapping robots [8,9,10,11,12,13] and fixed rubber tapping machines [14,15,16,17] instead of manual rubber tapping. They want to realize the precise tapping depth control of intelligent rubber tapping machinery by introducing industrial ranging technology involving machine vision [18,19,20], laser ranging [21,22] and ultrasonic ranging [9]. But these technologies can only guide machines to accurately determine the position of the outer surface of the bark and cannot further quickly detect bark thickness. Therefore, the tapping depth of intelligent rubber tapping machinery will be relatively conservative, resulting in its rubber production being less than 60% of manual rubber tapping. Given this, the portable electric tools for manual rubber tapping currently are also a good option in regard to enhancing the efficiency of rubber harvesting operations. All of those rubber tapping machineries can reduce the physical burden on workers and improve the overall productivity of the rubber harvesting process.
Vibration cutting is a process of multiple micro cutting operations. In the cutting cycle, the tools obtain significant instantaneous velocity and acceleration and generate high energy, reducing the plastic deformation and friction coefficient of the cutting area [23]. During vibration cutting, the shear angle shows a significant increasing trend, which reduces cutting force [24]. Vibration cutting only performs effective cutting in a short period of time and has a relatively long heat dissipation time, so the temperature in the cutting area is relatively low. Because of these advantages, vibration cutting is widely used in fields such as mechanical processing [25,26], agricultural production [27] and surgical operation [28].
Based on the principle of vibration cutting, the tapping parts of the 4GXJ-2 portable electric rubber cutter (Figure 2a) and automatic rubber tapping robot (Figure 2b) are developed by our team at the Chinese Academy of Tropical Agricultural Science Rubber Research Institute with Self Moving Technology Co., Ltd (Beijin, China). This vibration tool holder (Figure 2c) uses the blushless motor to drive the planetary rotor to eccentrically revolve at a constant speed and circumferentially impact the driving fork, and in this way the circular motion of the blushless motor is converted into the high-frequency reciprocating motion of the tool (Figure 2d). To reduce friction with the drive fork, the planetary rotor is designed to be self rotatable during eccentric revolution. When the planetary rotor becomes unstable, the tool will also experience unstable vibrations. This unstable vibration will form vibration ripples on the machining surface, leading to an increase in surface roughness [29]. It can also accelerate tool wear and may even cause the cutting edge to crack [30]. Meanwhile, the noise can be caused and may have harm on the operator.
To simplify the structure, a small rolling bearing is used as the planetary rotor in the vibration tool holder. More research focuses on the dynamic stability of rolling bearing in the fixed axis rotation state. Chinta et al. [31] established a rolling bearing dynamics model using a set of dimensionless differential equations and described the unbalanced response of the system. Tiwari et al. [32] established a nonlinear dynamic model of rolling bearing and analyzed the influence of the radial clearance size on the dynamic stability. Rafsanjani et al. [33] proposed an analytical model to study the nonlinear dynamic behavior of rolling bearing systems and used the Floquet theory to study the chaotic behavior of systems. Tu [34] considered the dynamic contact relationship between the internal components of rolling bearings and established its nonlinear dynamic model to discuss the changes in dynamic characteristics under variable speed rotation. Because the two highly complex rotational movements make system inertia, stiffness, and damping become more complex, the dynamic stability of rolling bearing with variable speed self rotation and eccentric revolution are not further researched.
The REF model is widely used in the research into dynamic characteristics of rotating flexible cylindrical parts of flexible gear [35], tires [36] and rolling bearing [37,38]. Based on the ANCF, the 3D-beam element and 3D REF [39], a novel eccentric 3D REF model is proposed to analyze the dynamic stability of the planetary rotor. By introducing multiple coordinate systems, the coupled motion of uniform eccentric revolution, variable speed self rotation and flexible deformation is decomposed and the influence of these motions on the centrifugal force and Coriolis force is more clearly derived.
The paper is structured as follows. Section 1 is the introduction. Section 2 is the materials and methods. In Section 3, the model is degraded and validated by comparing with other examples, and the stability of the planetary rotor in the vibration tool holder is analyzed and its unstable eccentric revolving speed is predicted. Section 4 is then the summary.

2. Methods and Materials

2.1. Overview of the Novel Eccentric 3D-REF Model

According to the structure of rolling bearing (Figure 3a), a novel eccentric 3D-REF model is established (Figure 3b). In the eccentric 3D-REF model, the inner ring, rolling ball and retainer of rolling bearing are simplified as spring/damping groups; the outer ring is simplified as a spatial ring with a rectangular cross-section. The model is divided into n ANCF 3D-beam elements along the circumferential directions of the basal axis with n + 1 node N i ( i = 1 , 2 , n + 1 ), and the radius of the basal axis is R (Figure 3c).
The Cartesian coordinate system O ¯ X ¯ Y ¯ Z ¯ is a fixed coordinate system, with its origin O ¯ being the center of eccentric revolution. The Cartesian coordinate system O X ˜ Y ˜ Z ˜ and O X Y Z are the co-moving coordinate system, and their origin O is the center of self rotation. The O ¯ X ¯ Y ¯ Z ¯ eccentrically revolves around the O ¯ at angular speed Ω ¯ , while the O X ˜ Y ˜ Z ˜ rotates around O at angular speed Ω ˜ . The dynamic process of an element can be decomposed into three movements:
(1)
The element undergoes flexible deformation in O X Y Z ;
(2)
The element rigidly self rotates in O X ˜ Y ˜ Z ˜ ;
(3)
The element rigidly reaches eccentric revolution in O ¯ X ¯ Y ¯ Z ¯ .
For the i - t h element (Figure 3d), its lateral thickness and radial thickness are a and b , respectively. The orthogonal curve coordinate system N i ξ ζ η is the elemental coordinate system, in which the ξ , ζ and η are the dimensionless coordinates of along the radial, circumferential, and lateral directions, respectively.
The base axis of the element is fixed on O by equivalent uniformly distributed radial, circumferential, and lateral spring-damping groups with stiffness k x and damping c x ; for the arbitrary node N i , the R i is the position coordinates and the R x i is the unit tangent vector along the x -aix in N i ξ ζ η ( x = ξ ,   ζ   and   η ).

2.2. The ANCF 3D-Beam Element

In the O X Y Z , the position field function of the ANCF 3D-beam element can be expressed by ξ , ζ and η as
X = A 0 + A 1 ζ + A 2 ζ 2 + A 3 ζ 3 + A 4 ξ + A 5 η + A 6 ξ ζ + A 7 ζ η Y = B 0 + B 1 ζ + B 2 ζ 2 + B 3 ζ 3 + B 4 ξ + B 5 η + B 6 ξ ζ + B 7 ζ η Z = C 0 + C 1 ζ + C 2 ζ 2 + C 3 ζ 3 + C 4 ξ + C 5 η + C 6 ξ ζ + C 7 ζ η
According to Equation (1), the elemental position field function can be expressed as
R = S ξ , ζ , η q i t
where the q i t are the elemental generalized coordinates and the S ξ , ζ , η is the elemental shape function.
The q i t is a 24 × 1 column matrix, and it can be written as
q i = R i T   R ξ i T   R ζ i T   R η i T   R i + 1 T   R ξ i + 1 T   R ζ i + 1 T   R η i + 1 T T
The initial elemental generalized coordinates q i 0 can be expressed as
q i 0 = R i 0 T   R ξ i 0 T   R ζ i 0 T   R η i 0 T   R i + 1 0 T   R ξ i + 1 0 T   R ζ i + 1 0 T   R η i + 1 0 T T
where
R i 0 = R cos 2 π i 1 n R sin 2 π i 1 n 0 T
R ξ i 0 = cos 2 π i 1 n sin 2 π i 1 n 0 T
R ξ i 0 = sin 2 π i 1 n cos 2 π i 1 n 0 T
R η i 0 = 0 0 1 T
The elemental shape function S ξ , η is a 3 × 24 matrix, and it can be expressed as
S = S 1 S 2 S 3 = S 1 I b S 2 I 2 π R n S 3 I a S 4 I S 5 I b S 6 I 2 π R n S 7 I a S 8 I
S 1 = 1 3 ζ 2 + 2 ζ 3 S 2 = ξ ξ ζ S 3 = ζ 2 ζ 2 + ζ 3 S 4 = η ζ η S 5 = 3 ζ 2 2 ζ 3 S 6 = ξ ζ S 7 = ζ 3 ζ 2 S 8 = ζ η
where the I is a third-order identity matrix.
Assuming that the initial and deformed global generalized coordinates of the planetary rotor are q 0 and q t , respectively, by introducing the Boolean matrix B i , the initial and deformed elemental generalized coordinates of q i 0 and q i can be represented as
q i 0 = B i q 0   and   q i = B i q t
In the O X ˜ Y ˜ Z ˜ , after the planetary rotor rotates ϕ ˜ t counterclockwise around O , the elemental position field function can be expressed as
R ˜ = T ˜ ( t ) R = T ˜ ( t ) S ξ , ζ , η B i q
In the O ¯ X ¯ Y ¯ Z ¯ , after the center point O of planetary rotor rotates ϕ ¯ t counterclockwise around O ¯ , the elemental position field function can be expressed as
R ¯ = T ¯ ( t ) R ˜ + T ¯ ( t ) r
where the r is the initial vector of O ¯ O and r = λ 0 0 T ; the T ˜ and T ¯ are the transformation matrices of rotation and eccentric revolution, respectively, which can be expressed as
T ˜ t = cos ϕ ˜ t sin ϕ ˜ t 0 sin ϕ ˜ t cos ϕ ˜ t 0 0 0 1   and   T ¯ t = cos ϕ ¯ t sin ϕ ¯ t 0 sin ϕ ¯ t cos ϕ ¯ t 0 0 0 1
For the convenience of later description, it is defined as
A ξ = A ξ ,   A ζ = A ζ ,   A η = A η ,   A ˙ = d A d t   and   A ¨ = d 2 A d t 2
where the A is an arbitrary vector or matrix.

2.3. The Nonlinear Dynamic Differential Equation of Novel Eccentric 3D REF Model

2.3.1. Kinetic Energy of Planetary Rotor

The elemental kinetic energy can be expressed as
T e = 1 2 ρ V i 0 1 0 1 1 1 R ¯ ˙ T R ¯ ˙ d ξ d ζ d η = 1 2 M i V e R ¯ ˙ T R ¯ ˙ d V e
where the ρ is the density; the V i is the elemental volume, V i = a 2 + 2 a R π / n ; the M i is elemental mass.
According to Equation (15), it can be observed that
d d t T e q ˙ T T e q T = M i V e d d t R ¯ ˙ T q ˙ T R ¯ ˙ d V e M i V e R ¯ ˙ T q T R ¯ ˙ d V e
According to Equations (12) and (13), it can be observed that
R ¯ ˙ = T ¯ ˙ T ˜ R + T ¯ T ˜ ˙ R + T ¯ T ˜ R ˙ + T ¯ ˙ r = T ¯ ˙ T ˜ S ξ , ζ , η B i q + T ¯ T ˜ ˙ S ξ , ζ , η B i q + T ¯ T ˜ S ξ , ζ , η B i q ˙ + T ¯ ˙ r
Substituting Equation (17) into Equation (16), it can be expressed as
d d t T e q ˙ T T e q T = M e 1 q ¨ + 2 ϕ ¯ ˙ + ϕ ˜ ˙ M e 2 q ˙ ϕ ¯ ˙ + ϕ ˜ ˙ 2 M e 1 q ϕ ¯ ¨ + ϕ ˜ ¨ M e 2 q ϕ ¯ ˙ 2 B i T S T T ˜ T r + ϕ ¯ ¨ B i T S T T ˜ ^ r
where
M e 1 = M i V e B i T S T S B i d V e , M e 2 = M i V e B i T S T I ^ S B i d V e
I ^ = 0 1 0 1 0 0 0 0 1 ,   T ˜ ^ = sin ϕ ˜ cos ϕ ˜ 0 cos ϕ ˜ sin ϕ ˜ 0 0 0 1
The derivation process of Equations (18)–(20) is listed in Appendix A.
In the Refs. [39,40], the rotation motion is represented by a rotation matrix of the global generalized coordinates, and the rotation matrix can not be easily eliminated. It will cause a large number of matrices with trigonometric functions to participate in the numerical calculation of dynamic equations, and the truncation error of these trigonometric functions will reduce the calculation accuracy. However, in this paper, the rotation matrices of the position field are used to represent the rotation and eccentric revolution. In the calculation process, many rotation matrices will be eliminated by T ˜ T T ˜ = I and T ¯ T T ¯ = I so that fewer matrices with trigonometric functions can participate in the calculation process.
The global kinetic energy of the planetary rotor can be expressed as
T = i = 1 n T e
When the revolution is at a constant speed and the rotation is at a variable speed, it can be obtained as ϕ ¯ ¨ = Ω ¯ ˙ = 0 , ϕ ¯ ˙ = Ω ¯ , ϕ ˜ ¨ = Ω ˜ ˙ and ϕ ˜ ˙ = Ω ˜ . With M 1 = i = 1 n M e 1 and M 2 = i = 1 n M e 2 . According to Equation (21), it can be obtained as
d d t T q ˙ T T q T = i = 1 n d d t T e q ˙ T T e q T = M q ¨ + C T q ˙ + K T q + Q T
M = M 1 C T = 2 Ω ¯ + Ω ˜ M 2 K T = Ω ¯ + Ω ˜ 2 M 1 Ω ˜ ˙ M 2 Q T = i = 1 n 0 1 M i Ω ¯ 2 V e B i T S T T ˜ T r d V e

2.3.2. Strain Energy of Planetary Rotor

In the orthogonal curvilinear coordinate system N i ξ ζ η , the Lagrange strain tensor can be expressed as
ε = 1 2 J 0 T J T J J 0 1 I
where the J is the deformation gradient tensor and the J 0 is initial deformation gradient tensor, which can be seen as the residual strain in the bending configuration.
Because flexible deformation only occurs in coordinate system O X Y Z , the J and J 0 can be expressed as
J = R ξ R ζ R η = S 1 ξ B i q S 1 ζ B i q S 1 η B i q S 2 ξ B i q S 2 ζ B i q S 2 η B i q S 3 ξ B i q S 3 ζ B i q S 3 η B i q
J 0 = R ξ 0 R ζ 0 R η 0 = S 1 ξ B i q 0 S 1 ζ B i q 0 S 1 η B i q 0 S 2 ξ B i q 0 S 2 ζ B i q 0 S 2 η B i q 0 S 3 ξ B i q 0 S 3 ζ B i q 0 S 3 η B i q 0
By calculations, it can be seen that the J 0 is an orthogonal constant matrix, so it can be set to
J 0 1 = J 0 T = C 11 C 12 0 C 21 C 22 0 0 0 1
According to Equations (25) and (27), Equation (24) can be rewritten as
ε = ε 11 ε 12 ε 13 ε 21 ε 22 ε 23 ε 31 ε 32 ε 33
where
ε 11 = C 11 2 Λ 1 + 2 C 11 C 21 Λ 4 + C 21 2 Λ 2 1 ε 22 = C 12 2 Λ 1 + 2 C 12 C 22 Λ 4 + C 22 2 Λ 2 1 ε 33 = Λ 3 1 ε 13 = ε 31 = C 11 Λ 6 + C 21 Λ 5 ε 23 = ε 32 = C 12 Λ 6 + C 22 Λ 5 ε 12 = ε 21 = C 11 C 12 Λ 1 + C 12 C 21 + C 11 C 22 Λ 4 + C 21 C 22 Λ 2
And
Λ = J T J = Λ 1 Λ 4 Λ 6 Λ 4 Λ 2 Λ 5 Λ 6 Λ 5 Λ 3
i.e.,
Λ 1 = q T B i T S 1 ξ T S 1 ξ B i q + q T B i T S 2 ξ T S 2 ξ B i q + q T B i T S 3 ξ T S 3 ξ B i q Λ 2 = q T B i T S 1 ζ T S 1 ζ B i q + q T B i T S 2 ζ T S 2 ζ B i q + q T B i T S 3 ζ T S 3 ζ B i q Λ 3 = q T B i T S 1 η T S 1 η B i q + q T B i T S 2 η T S 2 η B i q + q T B i T S 3 η T S 3 η B i q Λ 4 = q T B i T S 1 ξ T S 1 ζ B i q + q T B i T S 2 ξ T S 2 ζ B i q + q T B i T S 3 ξ T S 3 ζ B i q Λ 5 = q T B i T S 1 ζ T S 1 η B i q + q T B i T S 2 ζ T S 2 η B i q + q T B i T S 3 ζ T S 3 η B i q Λ 6 = q T B i T S 1 ξ T S 1 η B i q + q T B i T S 2 ξ T S 2 η B i q + q T B i T S 3 ξ T S 3 η B i q
The column matrix form of the Lagrange strain tensor ε ^ can be expressed as
ε ^ = ε 11 ε 22 ε 33 ε 23 ε 13 ε 12 T
The strain energy of the planetary rotor can be expressed as
U E q = i = 1 n 1 2 V i V e ε ^ T E ε ^ d V e
where
E = E 1 ν 1 + ν 1 2 ν E ν 1 ν 1 + ν 1 2 ν E ν 1 ν 1 + ν 1 2 ν 0 0 0 E ν 1 ν 1 + ν 1 2 ν E 1 ν 1 + ν 1 2 ν E ν 1 ν 1 + ν 1 2 ν 0 0 0 E ν 1 ν 1 + ν 1 2 ν E ν 1 ν 1 + ν 1 2 ν E 1 ν 1 + ν 1 2 ν 0 0 0 0 0 0 E 1 + ν 0 0 0 0 0 0 E 1 + ν 0 0 0 0 0 0 E 1 + ν
where the E is the elastic modulus, the ν is the Poisson’s ratio.
The generalized elastic force and stiffness of the planetary rotor can be expressed as
Q E = U E q q T
K E = 2 U E q q T q
The derivation process of Equations (35) and (36) are listed in Appendix B.

2.3.3. Boundary Condition

The spring/damping groups are introduced on the base axis of the planetary rotor to simulate boundary conditions. According to Hooke’s law, the density function of strain energy of the spring/damping groups can be written as
U S e = 1 2 k ξ u 2 + 1 2 k ζ v 2 + 1 2 k η w 2 + u c ξ u ˙ + v c ζ v ˙ + w c η w ˙
where the u , v and w are the radial, circumferential and lateral displacement of arbitrary point on the base axis, respectively.
Based on geometric relationships, it can be inferred that
R ˜ R ˜ 0 ζ = 0 η = 0 = u R ˜ ξ 0 + v R ˜ ζ 0 + w R ˜ η 0 ζ = 0 η = 0
So
R ˜ ˙ ζ = 0 η = 0 = u ˙ R ˜ ξ 0 + v ˙ R ˜ ζ 0 + w ˙ R ˜ η 0 ζ = 0 η = 0
Deformation potential energy of all spring/damping groups can be expressed as
W S q ,   q ˙ = 1 2 i = 1 n 2 π R n 0 1 Δ R ˜ T T ˜ k κ T ˜ k T Δ R ˜ + Δ R ˜ T T ˜ k c T ˜ k T R ˜ ˙ ζ = 0 η = 0 d ξ
where κ = diag k ξ k ζ k η , c = diag c ξ c ζ c η T k = R ˜ ξ 0 R ˜ ζ 0 R ˜ η 0 and Δ R ˜ = R ˜ R ˜ 0 .
The generalized elastic force and stiffness of all spring/damping groups can be expressed as
Q S q ,   q ˙ = W S q T d d t W S q ˙ T = i = 1 n 0 1 2 π R n B i T S T T k κ T k T S   B i q q 0 B i T S T T k c T k T S   B i q ˙ ζ = 0 η = 0 d ξ
K S = 2 U S q T q = i = 1 n 0 1 2 π R n B i T S T T k κ T k T S   B i ζ = 0 η = 0 d ξ C S = 2 U S q T q ˙ = i = 1 n 0 1 2 π R n B i T S T T k c T k T S   B i ζ = 0 η = 0 d ξ

2.3.4. External Load

The generalized cutting load can be expressed as
Q C = B i T S T F C ζ = 1 ,   η = 0 i = α ,   ξ = β
where the F C is the vector of the cutting load and the parameters α and β can determined based on the following methods.
According to the reciprocating cutting method, there are two corresponding contact situations between the planetary rotor and the drive fork, as shown in Figure 4.
In the contact situation 1, the point O is on the right side of the axis Y ¯ . Because the position vector of point O is d cos Ω ¯ t d sin Ω ¯ t 0 T in O ¯ X ¯ Y ¯ Z ¯ , it indicates the occurrence of contact situation 1 with cos Ω ¯ t 0 . The X ¯ value of node N α and N α + 1 must be the two smallest values in all nodes, as shown in Figure 4a. The α can be determined by searching for the X ¯ value of all nodes. At this time, F C = F C 0 0 T . If the nodal generalized coordinates of N α and N α + 1 are set to
N α = e 1 e 2 e 3 e 4 e 5 e 6 e 7 e 8 e 9 e 10 e 11 e 12 T N α + 1 = e 1 e 2 e 3 e 4 e 5 e 6 e 7 e 8 e 9 e 10 e 11 e 12 T
The parameter β can be determined by position field function of N α and N α + 1 , as follows
1 3 β 2 + 2 β 3 e 2 + 2 π R n β 2 β 2 + β 3 e 5 + 3 β 2 2 β 3 e 2 + 2 π R n β 2 + β 3 e 5 = d sin Ω ¯ t
Similarly, for the planetary rotor in the contact situation 2, there is cos Ω ¯ t 0 and F C = F C 0 0 T . Moreover, the X ¯ value of node N α 1 and N α must be the two largest values in all nodes, as shown in Figure 4b. The parameter β can be calculated by Equation (44).

2.3.5. Nonlinear Dynamic Differential Equation

The Lagrange equations can be expressed as
d d t L q ˙ T L q T + Q C = 0
where the L is Lagrange function of the planetary rotor, and it can be expressed as
L q ˙ , q = T q ˙ , q U E q W S q
According to Equations (22), (40), (43), (45) and (46), it can be seen that
M q ¨ 1 + C T t q ˙ 2 + K T t q 3 + Q S q , q ˙ 4 + Q E q 5 + Q T t 6 + Q C t 7 = 0
For the convenience of later description, we have numbered each term of Equation (47).

2.4. Numerical Calculation Method for Characteristic Multipliers

In this section, the Runge−Kutta method with Gill coefficients is used to derive the characteristic multipliers Φ T of the nonlinear dynamic differential equations with time-varying coefficients; then, the stability of planetary rotors is explained according to the Floquet theory.
Because some coefficients in Equation (47) are time-varying, its analytical expression cannot be obtained. Equation (47) can be rewritten as
y ˙ + ψ t y = f t ,   y
where
y ˙ = q ˙ q ¨ , y = q q ˙ , ψ t = C T t + C S M 1 M 1 0 1 K T t + K S 0 0 M 1 f t ,   y = C T t + C S M 1 M 1 0 1 Q T t Q C t Q E I y Q s I 1 y ,   I 2 y 0
with
I 1 = I 12 n m + 1 × 12 n m + 1 0 12 n m + 1 × 12 n m + 1 0 12 n m + 1 × 12 n m + 1 0 12 n m + 1 × 12 n m + 1 ,   I 2 = 0 12 n m + 1 × 12 n m + 1 0 12 n m + 1 × 12 n m + 1 0 12 n m + 1 × 12 n m + 1 I 12 n m + 1 × 12 n m + 1
By using the fourth-order Runge−Kutta method with Gill coefficients, it can be seen that
y t i + 1 = y t i + Δ t 6 χ 1 + 2 2 χ 2 + 2 + 2 χ 3 + χ 4
where the Δ t = t i + 1 t i and
χ 1 = ψ t i y i χ 2 = ψ ^ 1 t i y i χ 3 = ψ ^ 2 t i y i χ 4 = ψ ^ 3 t i y i
with
ψ t i y i = f t i ,   y i ψ ^ 1 t i = ψ t i + h 2 + 1 2 Δ t ψ t i ψ ^ 2 t i = ψ t i + h 2 + 2 1 2 Δ t ψ t i + 2 2 2 ψ ^ 1 t i ψ ^ 3 t i = ψ t i + h 2 2 Δ t ψ ^ 1 t i + 2 + 1 2 Δ t ψ ^ 2 t i
According to Equations (49)–(51), it can be observed that
y t i + 1 = Ξ t i y t i
where
Ξ t i = I + Δ t 6 ψ t i + 2 2 ψ ^ 1 t i + 2 + 2 ψ ^ 2 t i + ψ ^ 3 t i
where the I is a 24n (m + 1) order identity matrix.
Combining Equation (53), it can be inferred that
y T = k = 1 p Ξ T k Δ t y 0 = Φ T y 0
where the T is a motion period of planetary rotor and the k is a natural number.
If the real part of eigenvalues of Φ T is λ i , Equation (54) can be transformed into y T = λ y 0 ( λ = diag λ 1 λ 2 λ 24 n m + 1 ), and it can be derived that y kT = λ k y 0 . If the planetary rotor is stable, it will be satisfied that y kT y 0 with k .
According to the Floquet theory, it can be inferred that
Stable state: λ = 0
Unstable state: λ > 0

2.5. Initial Parameters of Planetary Rotor

In simplifying the bearing structure (Deep Groove Ball Bearing 625Z), the parameters of the planetary rotor are shown in Table 1.
Based on the finite element model, it can be the calculated that the 1st lateral frequency is f η 1 = 47.2   Hz , the 1st circumferential frequency is f ζ 1 = 17.1   Hz and the 1st radial frequency is f ξ 1 = 410.4   Hz .
According to the numerical relationship [41] between the axial/circumferential stiffness k η / k ζ and the 1st lateral/circumferential frequency f z 1 / f ξ 1 , it can be seen that
k η = 1 2 ρ a h 2 π f η 1 2 2060   Pa
k ζ = 1 2 ρ a h 2 π f ζ 1 2 270   Pa
According to the numerical relationship [42] between the radial stiffness k ξ and f ξ 1 , it can be observed that
f ξ 1 = 1 2 π 1 2 ρ A E a h 3 12 R + 2 k ξ + 2 k ζ
Through calculation Equation (57), the radial stiffness can be estimated as k ξ 6.22 × 10 6   Pa .
Because the bearing-damping is typically in the order of 0.25 2.5 × 10 5 times the linearized stiffness [43,44], which is chosen c ζ c η 0 and c ξ 150   Pa s .
The reason for the self rotation of the planetary rotor is the friction force when it comes into contact with the drive fork. Ignoring the sliding friction and contact deformation, it can be observed as
Ω ¯ R + a 2 + d cos Ω ¯ t = R + a 2 Ω ˜
The ideal rotating speed Ω ˜ can be rewritten as
Ω ˜ = Ω ¯ + 2 ε d 2 + ε a cos Ω ¯ t
with ε a = a / R and ε d = d / R .
As shown in Figure 5, by testing the output current of the brushless motor of 4GXJ-2 portable electric rubber cutter under low-speed revolution conditions ( Ω ¯ = 10   rad / s ), it can be observed that the output current of the brushless motor approximates the cosine function in the initial state; then, when the cutter starts to cut into the bark, the output current of the motor rapidly increases, and when the cutter is in the smooth-cutting state, the output current of the brushless motor approximates the cosine function.
Therefore, the cutting load can be expressed as
F C = F ¯ C + f cos Ω ¯ t
According to the testing of mechanical sensors, the average cutting force F ¯ C 48   N and the maximum cutting force is 62 N, i.e., f 14   N .

3. Verification and Experiment

3.1. Theoretical Verification

In this chapter, the eccentric 3D-REF mode is degraded to analyze the vibration characteristics of the uniformly rotating ring on elastic foundation or uniformly eccentric revolution annular plates. By comparing the results from other references, the validity of the proposed model can be verified.

3.1.1. Vibration of Uniformly Rotating Circular Ring Model

If the revolving speed Ω ¯ = 0 , the parameters of damping groups c ξ = c ζ = c η = 0 and cutting load F C = 0 , while the rotating speed Ω ˜ is a constant value, the novel eccentric 3D REF model can be degenerated into the free uniformly rotating circular ring model. Its degenerate dynamic equation can be rewritten as
M 1 q ¨ 1 + 2 Ω ˜ M 2 q ˙ 2 Ω ˜ 2 M q 3 + K S q 4 + Q E q 5 = 0
In Equation (61), the complete Item.1 and Item.5 of Equation (47) are retained, and so are the relevant parts of rotational speed in Item.2 and Item.3. In this degradation example, the validity of these terms is verified.
There are no coefficients of variation for time in Equation (61). Therefore, the equilibrium position q S can be obtained by solving Ω ˜ 2 M 1 q S + K S q S + Q E q S = 0 , and the vibration equation of the uniformly rotating circular ring model can be written as
ω 2 M 1 q ¨ + 2 ω Ω ˜ M 2 Ω ˜ 2 M 1 + K S + K E q S A = 0
where the ω is the natural frequency and the A is the eigenvector.
Because of the influence of the uniform rotation, in the O X ˜ Y ˜ Z ˜ , the natural frequency will split into two traveling waves as the rotating speed increases. According to the Doppler effect [39], the relationship between the natural circular frequencies in the O X Y Z and circular frequencies of the traveling wave in the O X ˜ Y ˜ Z ˜ can be written as
ω ˜ f = ω 1 + N Ω ˜ ω ˜ b = ω 2 N Ω ˜
where the ω 1 and ω 2 are a set solution of ω of Equation (63), which is also a set of natural circular frequencies in the O X Y Z ; the ω ˜ f and ω ˜ b are the circular frequencies of the forward and backward traveling waves in the O X ˜ Y ˜ Z ˜ , respectively.
The parameters of the rotating circular ring model are shown in Table 2 [45]. The forward and backward traveling wave circular frequencies ( ω ˜ f and ω ˜ b ) of this uniformly rotating circular ring model [45] and proposed model are shown in Figure 6. The calculation results of the two models are similar, which indicates the validity of the proposed model in calculating the rotating motion.

3.1.2. Rigid-Body Vibration Mode of Uniformly Eccentric Revolving Annular Plate

With Ω ˜ = 0 and F C = 0 , the novel eccentric 3D REF model is used to analyze the rigid-body vibration mode of the uniformly eccentric revolving annular plate, and its dynamic equation can be rewritten as
M 1 q ¨ 1 + 2 Ω ¯ M 2 q ˙ 2 Ω ¯ 2 M 1 q 3 + Q E q 5 + Q T 6 = 0
In Equation (64), the complete Item.1, Item.5 and Item.6 of Equation (47) are retained, as are the relevant parts of eccentric revolving speed in Item.2 and Item.3. In this degradation example, the validity of these terms is verified.
There is no variable coefficient of time in Equation (64), and the equilibrium position q S can be obtained by solving Ω ¯ 2 M 1 q S + K S q S + Q E q S + Q T = 0 . The vibration equation of the uniformly eccentric revolving annular plate can be written as
ω 2 M 1 q ¨ + ω 2 Ω ¯ M 2 Ω ¯ 2 M 1 + K S + K E q S A = 0
According to the parameters in Table 3 [46], the degradation model of the uniformly eccentric revolving annular plate was established. As shown in Figure 7, the calculation results of the rigid-body vibration frequency ω in the O X Y Z of the degradation model for this plate have a good consistency with it in Ref. [46], which indicates that the proposed model can effectively analyze the influences of eccentric revolution on vibration characteristics.

3.2. Stability Analysis and Experimental Analysis of the Planetary Rotor in the Rubber Tapping Machinery

Under the initial parameter conditions, the value of λ under different revolving speeds Ω ¯ is shown in Figure 8. When the revolving speed Ω ¯ is within the range of 0~1500 rad/s, this planetary rotor has three unstable regions of (817, 909), (1017, 1095) and (1263, 1312). In the unstable regions with lower revolving speed ranges, the maximum value of λ is larger.
Under different revolving speeds, the ten-cycle trajectory of nodes on the outer edge of the planetary rotor in O ¯ X ¯ Y ¯ Z ¯ is shown in Figure 9. When the revolving speed Ω ¯ is within the stable range (Figure 9b,d,f), the nodal radial jitter amplitude is less than 0.2mm and the ten-cycle nodal trajectories are regular and similar. However, the ten-cycle nodal trajectory will exhibit irregular tremors, accompanied by the revolving speed Ω ¯ being within the three unstable speeds (Figure 9a,c,e). When the revolving speed Ω ¯ are 850 rad/s, 1050 rad/s and 1290 rad/s, the maximum radial amplitude of nodal jitter are 3.2 mm, 2.7mm and 1.8 mm, respectively. At lower revolving speeds, the node jitter will be more pronounced. This indicates that the planetary rotor becomes more unstable with the increasing of λ .
According to the unstable speed domain predicted by the model, the motor speed of the 4GXJ-2 portable electric rubber cutter (which is the revolving speed Ω ¯ of planetary rotor) is adjusted to 750 rad/s, 850 rad/s, 950 rad/s, 1050 rad/s and 1150 rad/s, respectively, and rubber tapping experiments are conducted on the bark of rubber trees, as shown in Figure 10. At speeds of 750 rad/s (Figure 10a), 950 rad/s (Figure 10c) and 1150 rad/s (Figure 10e), the cutting surface of the bark is flat with little burrs. However, at speeds of 850 rad/s (Figure 10b) and 1050 rad/s (Figure 10d), which is in unstable regions of revolving speed, the cutting surface of bark is rough and presents a serrated shape; the jagged shape in cutting surface of 850 rad/s is more pronounced than that of 1050 rad/s. The experimental results are basically consistent with the calculated results of the model. The validity of the novel eccentric 3D REF model can be verified by the rubber tapping experiments.

4. Conclusions

For the stability analysis of a planetary rotor with variable speed self rotation and uniform eccentric revolutions in rubber tapping machinery, the novel eccentric 3D REF model of a planetary rotor is established by using an ANCF 3D beam element. Based on the degradation form of this model, the vibration characteristics of the rotating circular rings and eccentric rotating ring plates are analyzed. Then, combined with the Floquet theory and Runge−Kutta method, the stability of the vibration tool holder in a rubber tapping machine is analyzed and its unstable speed is predicted.
In the modeling process, an ANCF 3D beam element is first used to establish a 3D REF model. Compared with the spatial Euler−Bernoulli beam element, this element does not need to calculate high-order curvature terms, greatly reducing the nonlinearity of the model and increasing computational efficiency. By introducing multiple coordinate systems, the coupled motion of uniform eccentric revolution, variable speed self rotation and flexible deformation is decomposed; moreover, the centrifugal force and Coriolis force of these motions are more clearly derived. Furthermore, the rotation matrix of the position field is introduced to replace the rotation matrix of the generalized coordinates to describe eccentric revolution and self rotation, resulting in the numerical calculations of the trigonometric function being greatly reduced and avoiding numerous truncation errors interferring with the calculation results.
For the degradation models of the uniformly rotating circular ring and uniformly eccentric revolving annular plate, their calculation results of natural frequencies are similar to those of existing references. This indicates that the novel eccentric 3D REF model can validly analyze the dynamic characteristics of annular components with self rotational motion and eccentric running.
Through a dynamic stability characteristics analysis of a planetary rotor with the numerical calculation method for characteristic multipliers, the planetary rotor in a 4GXJ-2 portable electric rubber cutter (Rubber Research Institute, Chinese Academy of Tropical Agricultural Science, Haikou, China) will experience instability under the eccentric revolving speeds of (1263, 1312) (817, 909), (1017, 1095) and (1263, 1312) (unit: rad/s). By analyzing the nonlinear forced transient response of nodes within 10 revolution cycles, the prediction of the unstable revolution speed is further verified, and it is also gathered that the planetary rotor becomes more unstable with the increasing of λ . In the experiment of rubber tapping a 4GXJ-2 portable electric rubber cutter, the above conclusion is further demonstrated.
Through the calculation of a novel eccentric 3D REF model, the speed control switch of a 4GXJ-2 portable electric rubber cutter is optimized to two gears of low speed (750 rad/s) and high speed (1150 rad/s); the revolving speed of a vibration tool holder of an automatic rubber tapping robot is designed as 1500 rad/s.

Author Contributions

Conceptualization, J.C. and B.F.; methodology, B.F.; software, B.F.; validation, B.F., S.X. and X.S.; formal analysis, B.F.; investigation, B.F.; resources, J.C.; data curation, J.C.; writing—original draft preparation, B.F.; writing—review and editing, B.F; visualization, X.S.; supervision, S.X.; project administration, S.X. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the support of the Central Public-Interest Scientific Institution Basal Research Fund (No. 1630022022005); Supported by the Hainan Provincial Natural Science Foundation of China (No. 521QN0937) and Supported by the Hainan Provincial Natural Science Foundation of China (No. 323MS074).

Data Availability Statement

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

Acknowledgments

The financial support mentioned in the Funding part is gratefully acknowledged.

Conflicts of Interest

The authors declared no potential conflicts of interest with respect to the research, authorship, and publication of this article.

Abbreviations

ANCF: Absolute Nodal Coordinate Formulation; REF: Ring on Elastic Foundation.

Appendix A

According to Equation (17), it can be got
d d t R ¯ ˙ T q ˙ T R ¯ ˙ = d d t R ˙ T q ˙ T T ˜ T T ¯ T T ¯ ˙ T ˜ R + R ˙ T q ˙ T T ˜ T T ¯ T T ¯ T ˜ ˙ R + R ˙ T q ˙ T T ˜ T T ¯ T T ¯ T ˜ R ˙ + R ˙ T q ˙ T T ˜ T T ¯ T T ¯ ˙ r = R ˙ T q ˙ T T ˜ T T ¯ T T ¯ T ˜ R ¨ + 3 R ˙ T q ˙ T T ˜ T T ¯ T T ¯ ˙ T ˜ R ˙ + 3 R ˙ T q ˙ T T ˜ T T ¯ T T ¯ T ˜ ˙ R ˙ + 4 R ˙ T q ˙ T T ˜ ˙ T T ¯ T T ¯ ˙ T ˜ R + R ˙ T q ˙ T T ˜ T T ¯ ˙ T T ¯ ˙ T ˜ R + R ˙ T q ˙ T T ˜ ˙ T T ¯ T T ¯ T ˜ ˙ R + R ˙ T q ˙ T T ˜ T T ¯ T T ¯ ¨ T ˜ R + R ˙ T q ˙ T T ˜ T T ¯ T T ¯ T ˜ ¨ R
R ¯ ˙ T q T R ¯ ˙ = R T q T T ˜ ˙ T T ¯ T T ¯ T ˜ R ˙ + R T q T T ˜ T T ¯ ˙ T T ¯ T ˜ R ˙ + 2 R T q T T ˜ ˙ T T ¯ T T ¯ ˙ T ˜ R + R T q T T ˜ ˙ T T ¯ T T ¯ T ˜ ˙ R + R T q T T ˜ T T ¯ ˙ T T ¯ ˙ T ˜ R + R T q T T ˜ T T ¯ ˙ T T ¯ ˙ r + R T q T T ˜ ˙ T T ¯ T T ¯ ˙ r
According to Equation (14), it can be observed as
T ˜ T T ¯ T T ¯ T ˜ = 1 0 0 0 1 0 0 0 1 = I , T ˜ T T ¯ ˙ T T ¯ ˙ T ˜ = ϕ ¯ ˙ 2 1 0 0 0 1 0 0 0 1 = ϕ ¯ ˙ 2 I , T ˜ ˙ T T ¯ T T ¯ T ˜ ˙ = ϕ ˜ ˙ 2 1 0 0 0 1 0 0 0 1 = ϕ ˜ ˙ 2 I , T ˜ T T ¯ T T ¯ ˙ T ˜ = T ˜ T T ¯ ˙ T T ¯ T ˜ = ϕ ¯ ˙ 0 1 0 1 0 0 0 0 1 = ϕ ¯ ˙ I ^ , T ˜ T T ¯ T T ¯ T ˜ ˙ = T ˜ ˙ T T ¯ T T ¯ T ˜ = ϕ ˜ ˙ 0 1 0 1 0 0 0 0 1 = ϕ ˜ ˙ I ^ , T ˜ T T ¯ T T ¯ ˙ T ˜ ˙ = T ˜ T T ¯ ˙ T T ¯ T ˜ ˙ = T ˜ ˙ T T ¯ T T ¯ ˙ T ˜ = ϕ ¯ ˙ ϕ ˜ ˙ 1 0 0 0 1 0 0 0 1 = ϕ ¯ ˙ ϕ ˜ ˙ I , T ˜ T T ¯ T T ¯ ¨ T ˜ = ϕ ¯ ˙ 2 1 0 0 0 1 0 0 0 1 + ϕ ¯ ¨ 0 1 0 1 0 0 0 0 1 = ϕ ¯ ˙ 2 I ϕ ¯ ¨ I ^ , T ˜ T T ¯ T T ¯ T ˜ ¨ = ϕ ˜ ˙ 2 1 0 0 0 1 0 0 0 1 + ϕ ˜ ¨ 0 1 0 1 0 0 0 0 1 = ϕ ˜ ˙ 2 I ϕ ˜ ¨ I ^ , T ˜ T T ¯ T T ¯ ¨ = ϕ ¯ ˙ 2 cos ϕ ˜ sin ϕ ˜ 0 sin ϕ ˜ cos ϕ ˜ 0 0 0 1 + ϕ ¯ ¨ sin ϕ ˜ cos ϕ ˜ 0 cos ϕ ˜ sin ϕ ˜ 0 0 0 1 = ϕ ¯ ˙ 2 T ˜ T + ϕ ¯ ¨ T ˜ ^ , T ˜ T T ¯ ˙ T T ¯ ˙ = ϕ ¯ ˙ 2 cos ϕ ˜ sin ϕ ˜ 0 sin ϕ ˜ cos ϕ ˜ 0 0 0 1 = ϕ ¯ ˙ 2 T ˜ T , T ˜ ˙ T T ¯ T T ¯ ˙ = ϕ ˜ ˙ ϕ ¯ ˙ sin ϕ ˜ cos ϕ ˜ 0 cos ϕ ˜ sin ϕ ˜ 0 0 0 1 = ϕ ˜ ˙ ϕ ¯ ˙ T ˜ ^
Substituting Equations (A1)–(A3) into Equation (16), Equations (18)–(20) can be observed.

Appendix B

According to Equations (29)–(33), it can be observed that
ε ^ T q T = ε 11 q T ε 22 q T ε 33 q T ε 23 q T ε 13 q T ε 12 q T T ε ^ q = ε 11 q ε 22 q ε 33 q ε 23 q ε 13 q ε 12 q 2 ε ^ T q T q = 2 ε 11 q T q 2 ε 22 q T q 2 ε 33 q T q 2 ε 33 q T q 2 ε 13 q T q 2 ε 12 q T q
and
ε 11 q T = C 11 2 Λ 1 q T + 2 C 11 C 21 Λ 4 q T + C 21 2 Λ 2 q T 1 ε 22 q T = C 12 2 Λ 1 q T + 2 C 12 C 22 Λ 4 q T + C 22 2 Λ 2 q T 1 ε 33 q T = Λ 3 q T 1 ε 13 q T = ε 31 q T = C 11 Λ 6 q T + C 21 Λ 5 q T ε 23 q T = ε 32 q T = C 12 Λ 6 q T + C 22 Λ 5 q T ε 12 q T = ε 21 q T = C 11 C 12 Λ 1 q T + C 12 C 21 + C 11 C 22 Λ 4 q T + C 21 C 22 Λ 2 q T
ε 11 q = C 11 2 Λ 1 q + 2 C 11 C 21 Λ 4 q + C 21 2 Λ 2 q 1 ε 22 q = C 12 2 Λ 1 q + 2 C 12 C 22 Λ 4 q + C 22 2 Λ 2 q 1 ε 33 q = Λ 3 q 1 ε 13 q = ε 31 q = C 11 Λ 6 q + C 21 Λ 5 q ε 23 q = ε 32 q = C 12 Λ 6 q + C 22 Λ 5 q ε 12 q = ε 21 q = C 11 C 12 Λ 1 q + C 12 C 21 + C 11 C 22 Λ 4 q + C 21 C 22 Λ 2 q
2 ε 11 q T q = C 11 2 2 Λ 1 q T q + 2 C 11 C 21 2 Λ 4 q T q + C 21 2 2 Λ 2 q T q 1 2 ε 22 q T q = C 12 2 2 Λ 1 q T q + 2 C 12 C 22 2 Λ 4 q T q + C 22 2 2 Λ 2 q T q 1 2 ε 33 q T q = 2 Λ 3 q T q 1 2 ε 13 q T q = 2 ε 31 q T q = C 11 2 Λ 6 q T q + C 21 2 Λ 5 q T q 2 ε 23 q T q = 2 ε 32 q T q = C 12 2 Λ 6 q T q + C 22 2 Λ 5 q T q ε 12 = ε 21 = C 11 C 12 2 Λ 1 q T q + C 12 C 21 + C 11 C 22 2 Λ 4 q T q + C 21 C 22 2 Λ 2 q T q
where
Λ 1 q T = 2 B i T S 1 ξ T S 1 ξ B i q + 2 B i T S 2 ξ T S 2 ξ B i q + 2 B i T S 3 ξ T S 3 ξ B i q Λ 2 q T = 2 B i T S 1 ζ T S 1 ζ B i q + 2 B i T S 2 ζ T S 2 ζ B i q + 2 B i T S 3 ζ T S 3 ζ B i q Λ 3 q T = 2 B i T S 1 η T S 1 η B i q + 2 B i T S 2 η T S 2 η B i q + 2 B i T S 3 η T S 3 η B i q Λ 4 q T = 2 B i T S 1 ξ T S 1 ζ B i q + 2 B i T S 2 ξ T S 2 ζ B i q + 2 B i T S 3 ξ T S 3 ζ B i q Λ 5 q T = 2 B i T S 1 ζ T S 1 η B i q + 2 B i T S 2 ζ T S 2 η B i q + 2 B i T S 3 ζ T S 3 η B i q Λ 6 q T = 2 B i T S 1 ξ T S 1 η B i q + 2 B i T S 2 ξ T S 2 η B i q + 2 B i T S 3 ξ T S 3 η B i q
2 Λ 1 q T q = 2 B i T S 1 ξ T S 1 ξ B i + 2 B i T S 2 ξ T S 2 ξ B i + 2 B i T S 3 ξ T S 3 ξ B i 2 Λ 2 q T q = 2 B i T S 1 ζ T S 1 ζ B i + 2 B i T S 2 ζ T S 2 ζ B i + 2 B i T S 3 ζ T S 3 ζ B i 2 Λ 3 q T q = 2 B i T S 1 η T S 1 η B i + 2 B i T S 2 η T S 2 η B i + 2 B i T S 3 η T S 3 η B i 2 Λ 4 q T q = 2 B i T S 1 ξ T S 1 ζ B i + 2 B i T S 2 ξ T S 2 ζ B i + 2 B i T S 3 ξ T S 3 ζ B i 2 Λ 5 q T q = 2 B i T S 1 ζ T S 1 η B i + 2 B i T S 2 ζ T S 2 η B i + 2 B i T S 3 ζ T S 3 η B i 2 Λ 6 q T q = 2 B i T S 1 ξ T S 1 η B i + 2 B i T S 2 ξ T S 2 η B i + 2 B i T S 3 ξ T S 3 η B i
According to Equation (37), it can be observed that
Q E q = U E q T = i = 1 n j = 1 m 0 1 0 1 A i j h ε ^ T q T E ε ^ d ξ d η
K E q = 2 U E q T q = i = 1 n j = 1 m 0 1 0 1 A i j h 2 ε ^ T q T q E ε ^ + ε ^ T q T E 2 ε ^ q E ε ^ d ξ d η

References

  1. Farag, N.H.; Pan, J. Modal characteristics of in-plane vibration of circular plates clamped at the outer edge. J. Acoust. Soc. Am. 2003, 113, 1935–1946. [Google Scholar] [CrossRef]
  2. Kim, C.B.; Cho, H.S.; Beom, H.G. Exact solutions of in-plane natural vibration of a circular plate with outer edge restrained elastically. J. Sound Vib. 2012, 331, 2173–2189. [Google Scholar] [CrossRef]
  3. Wang, Z.M.; Wang, Z.; Zhang, R. Transverse vibration analysis of spinning circular plate based on differential quadrature method. J. Vib. Shock 2014, 33, 125–129. (In Chinese) [Google Scholar] [CrossRef]
  4. Żur, K.K. Free vibration analysis of elastically supported functionally graded annular plates via quasi-Green’s function method. Compos. Part B Eng. 2018, 144, 37–55. [Google Scholar] [CrossRef]
  5. Tian, W.M.; Shi, M.J.; Tan, H.Y.; Wu, J.L.; Hao, B.Z. Structure and Development of Rubber Tree Bark, 1st ed.; Science Press: Beijing, China, 2015. (In Chinese) [Google Scholar]
  6. Zhang, X.R.; Wen, Z.T.; Zhang, Z.F.; Sun, Z.J.; Zhang, H.; Liu, J.X. Design and test of automatic rubber-tapping device with spiral movement. Trans. Chin. Soc. Agric. Mach. 2023, 54, 169–179. (In Chinese) [Google Scholar] [CrossRef]
  7. Soumahin, E.F.; Chantal, E.K.M.; N’dri, A.A.E.; Kouadio, Y.J.; Obouayeba, S. Influence of tapping time on rubber yield and the physiological status of rubber trees southeastern Côte d’Ivoire in a context of climate change. Int. J. Agric. Pol. Res. 2022, 10, 134–146. [Google Scholar] [CrossRef]
  8. Fang, J.; Shi, Y.; Cao, J.; Sun, Y.; Zhang, W. Active navigation system for a rubber-tapping robot based on trunk detection. Remote. Sens. 2023, 15, 3717. [Google Scholar] [CrossRef]
  9. Yang, H.; Sun, Z.; Liu, J.; Zhang, Z.; Zhang, X. The development of rubber tapping machines in intelligent agriculture: A review. Appl. Sci. 2022, 12, 9304. [Google Scholar] [CrossRef]
  10. Zhou, H.; Zhang, S.; Zhang, J.; Zhang, C.; Wang, S.; Zhai, Y.; Li, W. Design, development, and field evaluation of a rubber tapping robot. J. Field Robot. 2021, 39, 28–54. [Google Scholar] [CrossRef]
  11. Nie, F.; Zhang, W.; Wang, Y.; Shi, Y.; Huang, Q. A forest 3D Lidar SLAM system for rubber-tapping robot based on trunk center atlas. IEEE/ASME Trans. Mechatron. 2022, 27, 2623–2633. [Google Scholar] [CrossRef]
  12. Zhang, C.; Yong, L.; Chen, Y.; Zhang, S.; Ge, L.; Wang, S.; Li, W. A rubber-tapping robot forest navigation and information collection system based on 2D LiDAR and a gyroscope. Sensors 2019, 19, 2136. [Google Scholar] [CrossRef] [PubMed]
  13. Wang, S.; Zhou, H.; Zhang, C.; Ge, L.; Li, W.; Yuan, T.; Zhang, W.; Zhang, J. Design, development and evaluation of latex harvesting robot based on flexible toggle. Robot. Auton. Syst. 2022, 147, 103906. [Google Scholar] [CrossRef]
  14. Zhang, X.R.; Cao, C.; Zhang, L.N.; Xing, J.J.; Liu, J.X.; Dong, X.H. Design and Test of Profiling Progressive Natural Rubber Automatic Tapping Machine. Trans. Chin. Soc. Agric. Mach. 2022, 53, 99–108. (In Chinese) [Google Scholar] [CrossRef]
  15. Gao, K.K.; Sun, J.H.; Gao, F.; Jiao, J. Tapping error analysis and precision control of fixed tapping robot. Trans. Chin. Soc. Agric. Mach. 2021, 37, 44–50. (In Chinese) [Google Scholar] [CrossRef]
  16. Wen, X.Y.; Hong, Y.M.; Xiong, J.M.; Yu, Y.G.; Tan, J.H.; Meng, X.B. Design and experiment of fixed automatic natural rubber tapping machine. Trans. Chin. Soc. Agric. Mach. 2023, 54 (Suppl. S2), 128–135. (In Chinese) [Google Scholar] [CrossRef]
  17. Kamil, M.F.M.; Zakaria, W.N.W.; Tomari, M.R.M.; Sek, T.K.; Zainal, N. Design of automated rubber tapping mechanism. IOP Conf. Ser. Mater. Sci. Eng. 2020, 917, 012016. [Google Scholar] [CrossRef]
  18. Sun, Z.J.; Xing, J.J.; Hu, H.N.; Zhang, X.R.; Dong, X.H.; Deng, Y.G. Research on recognition and planning of tapping trajectory of natural rubber tree based on machine vision. J. Chin. Agric. Mech. 2022, 43, 102–108. (In Chinese) [Google Scholar] [CrossRef]
  19. Wongtanawijit, R.; Khaorapapong, T. Rubber tapping line detection in near-range images via customized YOLO and U-Net branches with parallel aggregation heads convolutional neural network. Neural Comput. Appl. 2022, 34, 20611–20627. [Google Scholar] [CrossRef]
  20. Wongtanawijit, R.; Khaorapapong, T. Nighttime rubber tapping line detection in near-range images. Multimed. Tools Appl. 2021, 80, 29401–29422. [Google Scholar] [CrossRef]
  21. Zhang, C.L.; Li, D.C.; Zhang, S.L.; Shui, Y.; Tan, Y.Z.; Li, W. Design and test of three-coordinate linkage natural rubber tapping device based on laser ranging. Trans. Chin. Soc. Agric. Mach. 2019, 50, 121–127. (In Chinese) [Google Scholar] [CrossRef]
  22. Yang, H.; Chen, Y.; Liu, J.; Zhang, Z.; Zhang, X. A 3D lidar SLAM system based on semantic segmentation for rubber-tapping robot. Forests 2023, 14, 1856. [Google Scholar] [CrossRef]
  23. Xu, W.X.; Zhang, L.C. Ultrasonic vibration-assisted machining: Principle, design and application. Adv. Manuf. 2015, 3, 173–192. [Google Scholar] [CrossRef]
  24. Bai, W.; Roy, A.; Guo, L.; Xu, J.; Silberschmidt, V.V. Silberschmidt. Analytical prediction of frictional behaviour and shear angle in vibration-assisted cutting. J. Manuf. Process. 2021, 62, 37–46. [Google Scholar] [CrossRef]
  25. Liang, X.L.; Zhang, C.B.; Cheung, C.F.; Wang, C.; Li, K.; Bulla, B. Micro/nano incremental material removal mechanisms in high-frequency ultrasonic vibration-assisted cutting of 316L stainless steel. Int. J. Mach. Tool Manu. 2023, 191, 104064. [Google Scholar] [CrossRef]
  26. Zhang, X.Y.; Peng, Z.L.; Liu, L.B. A transient cutting temperature prediction model for high-speed ultrasonic vibration turning. J. Manuf. Process. 2022, 83, 257–269. [Google Scholar] [CrossRef]
  27. Zou, L.; Yan, D.; Niu, Z.; Yuan, J.; Cheng, H.; Zheng, H. Parametric analysis and numerical optimisation of spinach root vibration shovel cutting using discrete element method. Comput. Electron. Agric. 2023, 212, 108138. [Google Scholar] [CrossRef]
  28. Shu, L.M.; Sugita, N. Analysis of fracture, force, and temperature in orthogonal elliptical vibration-assisted bone cutting. J. Comput. Electron. Agric. 2020, 103, 103599. [Google Scholar] [CrossRef]
  29. Chen, Y.; Pan, L.; Yin, Z.; Wu, Y. Effects of ultrasonic vibration-assisted machining methods on the surface polishing of silicon carbide. J. Mater. Sci. 2024, 59, 7700–7715. [Google Scholar] [CrossRef]
  30. Peng, Y.; Liang, Z.; Wu, Y.; Guo, Y.; Wang, C. Effect of vibration on surface and tool wear in ultrasonic vibration-assisted scratching of brittle materials. Int. J. Adv. Manuf. Technol. 2012, 59, 67–72. [Google Scholar] [CrossRef]
  31. Chinta, M.; Palazzolo, A.B. Stability and bifurcation of rotor motion in a magnetic bearing. J. Sound Vib. 1998, 214, 793–803. [Google Scholar] [CrossRef]
  32. Tiwari, M.; Gupta, K.; Prakash, O. Effect of radial internal clearance of a ball bearing on the dynamics of a balanced horizontal rotor. J. Sound Vib. 2000, 238, 723–756. [Google Scholar] [CrossRef]
  33. Rafsanjani, A.; Abbasion, S.; Farshidianfar, A.; Moeenfard, H. Nonlinear dynamic modeling of surface defects in rolling element bearing systems. J. Sound Vib. 2009, 319, 1150–1174. [Google Scholar] [CrossRef]
  34. Tu, W.B.; Liang, J.; Yu, W.N. Motion stability analysis of cage of rolling bearing under the variable-speed condition. Nonlinear Dyn. 2023, 111, 11045–11063. [Google Scholar] [CrossRef]
  35. Cooley, C.G.; Parker, R.G. Vibration of high-speed rotating rings coupled to space-fixed stiffnesses. J. Sound Vib. 2014, 333, 2631–2648. [Google Scholar] [CrossRef]
  36. Liu, Z.; Gao, Q. Development and parameter identification of the flexible beam on elastic continuous tire model for a heavy-loaded radial tire. J. Vib. Control 2018, 24, 5233–5248. [Google Scholar] [CrossRef]
  37. Liu, J.; Tang, C.K.; Shao, Y.M. An innovative dynamic model for vibration analysis of a flexible roller bearing. Mech. Mach. Theory 2019, 135, 27–39. [Google Scholar] [CrossRef]
  38. Viitala, R. Minimizing the bearing inner ring roundness error with installation shaft 3D grinding to reduce rotor subcritical response. CIRP J. Manuf. Sci. Technol. 2020, 30, 140–148. [Google Scholar] [CrossRef]
  39. Fan, B.; Wang, Z.M. Vibration analysis of radial tire using the 3D rotating hyperelastic composite REF based on ANCF. Appl. Math. Model. 2024, 126, 206–231. [Google Scholar] [CrossRef]
  40. Chen, X.Z.; Zhang, D.G.; Li, L. Dynamic analysis of rotating curved beams by using absolute nodal coordinate formulation based on radial point interpolation method. J. Sound Vib. 2019, 441, 63–83. [Google Scholar] [CrossRef]
  41. Kindt, P.; Sas, P.; Desmet, W. Development and validation of a three-dimensional ring-based structural tyre model. J. Sound Vib. 2009, 326, 852–869. [Google Scholar] [CrossRef]
  42. Wei, Y.T.; Liu, Z.; Zhou, F.Q.; Zhao, C.L. Three-dimensional REF model of tire including the out-of-plane vibration. J. Vib. Eng. 2016, 29, 795–803. (In Chinese) [Google Scholar] [CrossRef]
  43. Sopanen, J.; Mikkola, A. Dynamic model of a deep-groove ball bearing including localized and distributed defects. Part 1: Theory. Proc. Inst. Mech. Eng. Part K J. Multi-Body Dyn. 2003, 217, 201–211. [Google Scholar] [CrossRef]
  44. Sopanen, J.; Mikkola, A. Dynamic model of a deep-groove ball bearing including localized and distributed defects. Part 2: Implementation and results. Proc. Inst. Mech. Eng. Part K J. Multi-Body Dyn. 2003, 217, 213–223. [Google Scholar] [CrossRef]
  45. Vu, T.D.; Duhamel, D.; Abbadi, Z.; Yin, H.-P.; Gaudin, A. A nonlinear circular ring model with rotating effects for tire vibrations. J. Sound Vib. 2016, 388, 245–271. [Google Scholar] [CrossRef]
  46. Chung, J.T.; Heo, J.W.; Han, C.S. Natural frequencies of a spinning disk misaligned with the axis of rotation. J. Sound Vib. 2003, 260, 763–775. [Google Scholar] [CrossRef]
Figure 1. The traditional rubber trapping.
Figure 1. The traditional rubber trapping.
Forests 15 01071 g001
Figure 2. The vibration tool holder. (a) The vibration tool holder in 4GXJ-2 portable electric rubber cutter, (b) the vibration tool holder of automatic rubber tapping robot, (c) the partial internal components of vibration tool holder, (d) the kinematic sketch.
Figure 2. The vibration tool holder. (a) The vibration tool holder in 4GXJ-2 portable electric rubber cutter, (b) the vibration tool holder of automatic rubber tapping robot, (c) the partial internal components of vibration tool holder, (d) the kinematic sketch.
Forests 15 01071 g002
Figure 3. The eccentric 3D-REF model and element. (a) Deep groove ball bearing, (b) eccentric 3D-REF of planetary rotor, (c) vertical view of eccentric 3D-REF, (d) ANCF 3D-beam element.
Figure 3. The eccentric 3D-REF model and element. (a) Deep groove ball bearing, (b) eccentric 3D-REF of planetary rotor, (c) vertical view of eccentric 3D-REF, (d) ANCF 3D-beam element.
Forests 15 01071 g003
Figure 4. The two contact situations between the planetary rotor and the drive fork.
Figure 4. The two contact situations between the planetary rotor and the drive fork.
Forests 15 01071 g004
Figure 5. The output current test of brushless motor.
Figure 5. The output current test of brushless motor.
Forests 15 01071 g005
Figure 6. The forward and backward traveling wave circular frequencies ( ω ˜ f and ω ˜ b ) of rotating circular ring model [45] and proposed model.
Figure 6. The forward and backward traveling wave circular frequencies ( ω ˜ f and ω ˜ b ) of rotating circular ring model [45] and proposed model.
Forests 15 01071 g006
Figure 7. The comparison of calculation results of the rigid-body vibration frequency of eccentric revolving annular plate [46] and proposed model.
Figure 7. The comparison of calculation results of the rigid-body vibration frequency of eccentric revolving annular plate [46] and proposed model.
Forests 15 01071 g007
Figure 8. The unstable regions under different revolving speeds Ω ¯ .
Figure 8. The unstable regions under different revolving speeds Ω ¯ .
Forests 15 01071 g008
Figure 9. Under different revolving speeds, the ten-cycle nodal trajectories of planetary rotor in O ¯ X ¯ Y ¯ Z ¯ . (a) Ω ¯ = 850   rad / s . (b) Ω ¯ = 950   rad / s . (c) Ω ¯ = 1050   rad / s . (d) Ω ¯ = 1200   rad / s . (e) Ω ¯ = 1290   rad / s . (f) Ω ¯ = 1500   rad / s .
Figure 9. Under different revolving speeds, the ten-cycle nodal trajectories of planetary rotor in O ¯ X ¯ Y ¯ Z ¯ . (a) Ω ¯ = 850   rad / s . (b) Ω ¯ = 950   rad / s . (c) Ω ¯ = 1050   rad / s . (d) Ω ¯ = 1200   rad / s . (e) Ω ¯ = 1290   rad / s . (f) Ω ¯ = 1500   rad / s .
Forests 15 01071 g009
Figure 10. Rubber tapping experiments. (a) Ω ¯ = 750   rad / s . (b) Ω ¯ = 850   rad / s . (c) Ω ¯ = 950   rad / s . (d) Ω ¯ = 1050   rad / s . (e) Ω ¯ = 1150   rad / s .
Figure 10. Rubber tapping experiments. (a) Ω ¯ = 750   rad / s . (b) Ω ¯ = 850   rad / s . (c) Ω ¯ = 950   rad / s . (d) Ω ¯ = 1050   rad / s . (e) Ω ¯ = 1150   rad / s .
Forests 15 01071 g010
Table 1. Parameters of the planetary rotor.
Table 1. Parameters of the planetary rotor.
ParameterSymbolValueUnit
Radius R 6.8mm
Radial thickness a 2.5mm
Transverse thickness h 1.2mm
Elastic modulus of outer ring E 207GPa
Poisson’s ratio μ 0.3-
Density ρ 7810kg/m3
Eccentricity d 0.6mm
Table 2. Parameters of rotating circular ring model [45].
Table 2. Parameters of rotating circular ring model [45].
ParameterSymbolValueUnit
Radius R 0.285m
Radial thickness a 0.01mm
Transverse thickness h 0.085mm
Elastic modulus E 2.6611GPa
Poisson’s ratio μ 0.3-
Density ρ 1160.7kg/m3
Spring stiffness k ξ 4.35MPa
k ζ 0.319MPa
k η 0MPa
Table 3. Parameters of eccentric revolving annular plate [46].
Table 3. Parameters of eccentric revolving annular plate [46].
ParameterSymbolValueUnit
Radius R 40mm
Radial thickness a 50mm
Transverse thickness h 1.2mm
Elastic modulus E 6.55MPa
Poisson’s ratio μ 0.3-
Density ρ 1200kg/m3
Eccentricity d 13mm
Revolving speed Ω ˜ 1000rad/s
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Cao, J.; Fan, B.; Xiao, S.; Su, X. Stability Analysis of Planetary Rotor with Variable Speed Self Rotation and Uniform Eccentric Revolution in the Rubber Tapping Machinery. Forests 2024, 15, 1071. https://doi.org/10.3390/f15061071

AMA Style

Cao J, Fan B, Xiao S, Su X. Stability Analysis of Planetary Rotor with Variable Speed Self Rotation and Uniform Eccentric Revolution in the Rubber Tapping Machinery. Forests. 2024; 15(6):1071. https://doi.org/10.3390/f15061071

Chicago/Turabian Style

Cao, Jianhua, Bo Fan, Suwei Xiao, and Xin Su. 2024. "Stability Analysis of Planetary Rotor with Variable Speed Self Rotation and Uniform Eccentric Revolution in the Rubber Tapping Machinery" Forests 15, no. 6: 1071. https://doi.org/10.3390/f15061071

APA Style

Cao, J., Fan, B., Xiao, S., & Su, X. (2024). Stability Analysis of Planetary Rotor with Variable Speed Self Rotation and Uniform Eccentric Revolution in the Rubber Tapping Machinery. Forests, 15(6), 1071. https://doi.org/10.3390/f15061071

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop