Stability Analysis of Nonlinear Glue Flow Control System With Delay

In this paper, we improve a new mathematical model associated with glue ﬂow control system for glue applying of particleboard. Firstly, we study the existence and stability of the equilibria and the existence of fold, Hopf and Bogdanov-Takens bifurcations in above system. Next, the normal forms of Hopf bifurcation and Bogdanov-Takens bifurcation are derived, and the classiﬁcations of local dynamics near above bifurcation critical values are analyzed. Then, numerical simulation results show that the ﬂow control system associated with glue applying of particleborad exists stable equilibrium, stable periodic-1, periodic-2, and periodic-4 solutions, and chaotic attractor phenomenon from a sequence of period-doubling bifurcations. Final-ly, we compare the dynamical phenomena of ﬂow control system with and without cubic terms, showing that cubic terms can eﬀect the dynamical behaviors of ﬂow control system.


Introduction
With the progress of science and technology and the improvement of people's demand for wood quality in production and life, the technology of particleboard production continues to develop. Particleboard, also known as particle-board or chipboard, is a wood-based board material made of wood fragments, such as chips or shavings, that are mechanically pressed into sheet form and bonded together with resin. It is mainly used in furniture manufacturing, construction industry and carriage of train and car manufacturing. Compared with conventional wood and plywood, the advantages of particleboard is that it is cheaper, denser, and more stable and uniform. Thus, it can be processed into large-format plate with different styles and specifications as needed [1,2]. China is a country short of wood, and the development of particleboard production is an important way to improve the utilization rate of wood.
The flow chart of production process of particleboard is shown in Fig. 1. The furnish for particleboard normally consists of wood chips and logs are debarked, chipped and screened prior to washing. Then the chippers and wood particles are moved to the drying system, which is commonly used to reduce the moisture content of the particles to the desired levels. The wood particles are then mixed with an adhesive to bind them together, which is called glue applying process. This particle mix is then formed and cross cut, and the resulting mat is sent to the pre-press and hot-press devices and then to board finishing to create the end product particleboard. After pressing, the boards are unload, cooled, edge trimmed and sanded, to final dimensions. Finally, the finished product is packaged for shipment.
The process of glue applying of particleboard is one of the most important part in the production of parti-  cleboard, since the amount of glue directly effects the static strength, elastic mode, internal binding strength and other mechanical properties of particleboard, even influences the effect of hot-pressing and emission of volatile organic compound (VOC). Thus, it has important significant that studying the glue applying system of particleboard. Figure 2 shows the structure of particleboard glue applying system, and it consists of raw glue supply device (upper left part), particle supply device (upper right part) and mixing device (lower half part). The process of particleboard glue applying is as follows. Firstly, the raw glue and hardener, whose weight are measured by weighting sensor A, enter into mixing device through pipeline A. The rotating speed of pump decides the flow velocity of raw glue. Simultaneously, particle is transported on weighting belt through screw bin, whose weight and speed are measured by weighting sensor B and gyral coder, respectively. Then, particle is transported in the entrance for particle by belt and screw conveyor. Next, glue and particle are mixed in mixing device. Finally, mixed particle with glue is output from discharge chute. There exist various bifurcations in actual engineering problems, and delay differential equation may exhibit much more complicated behavior than ordinary differential equation (ODE) [3][4][5]. The automatic control process of glue applying system is to set the ratio of particle and glue according to the value of particle flow. The glue is pumped out from glue mixing tank by gear pump, then sprayed to particle. The amount of glue applying is determined by the speed of motor, and the speed of motor is measured and controlled by electromagnetic flow transmitter. The basic open loop structure of flow control system is given as Fig. 3. In this figure, squirrel cage asynchronous motor is powered by frequency converter, and we can adjust the speed of motor easily by changing the frequency of frequency converter. Flow sensor measures the flow in real time. Actually, we only study the inside part of dotted line for Fig. 3 in this paper. The glue comes out from the outlet of pump, and is piped to jet, which is a delayed process. The rest of the paper is organized as follows. In Sect. 2, we improve a new nonlinear delayed mathematical model associated with the glue flow control system of particleboard. In Sect. 3, we study the existence and stability of the equilibria and the existence of fold, Hopf, Bogdanov-Takens bifurcations in the delayed nonlinear system. In Sect. 4, we derive the normal forms of Hopf bifurcation and Bogdanov-Takens bifurcation, and also present the bifurcation analysis. In Sect. 5, we draw the numerical simulations. In Sect. 6, we compare the dynamical phenomena of flow control system with and without cubic terms (that is, frequency water supply system and glue flow control system of glue applying of particleboard). Finally, conclusion is shown in Sect. 7.

Mathematical Modeling
Consider some asynchronous motor, the torque T e produced by the interaction of field and armature currents in the machine will result in the angular velocity ω and will be opposed by the load torque T L on the shaft. If the motion is rotational, which is usually the case with electrical drives. Under dynamic conditions when the speed is changing, we obtain the following motion equation of electrical drive [6,7]: with T e being the driving-and T L the load torque. ω is the angular velocity, J is the moment of inertia of the rotating mass about the axis of rotation. Jω is the angular momentum. The term ω dJ dt is significance with variable inertia drives such as centrifuges or reeling drives, where the geometry of the load depends on speed or time, or industrial robots with changing geometry. In most cases however, the inertia can be assumed to be constant, hence, we have It should be noted that T e is the internal or electrical motor torque, not identical with the torque available at the motor shaft. The difference between internal torque and shaft torque is the torque required for accelerating the inertia of the motor itself and overcoming the internal friction torque of the motor. When the speed of pump is constant, hydraulic pressure will be increased with the decrease of flow, and will be decreased as the flow increasing. Comparing with constant pressure or constant speed water supplying systems, pump head will be saved in variable pressure water supplying system. For pump load, when the ratio of power supply voltage and frequency is a constant, the linearized expressions for torque T e and load torque T L are given as follows: T e = pK 0 ω 0 − (D + pK 0 )ω, where ω 0 is the deviations of power frequency for the stator, (Hz); D > 0 is friction coefficient (dimensionless); m is the deviation of load torque, (N · m); K 0 = p r2 ( V10 ω10 ) 2 , where V 10 and ω 10 are the values on static work point of voltage (V ) and frequency (Hz) of the stator power supply, respectively; p > 0 is pole pairs (dimensionless); r 2 > 0 is rotor resistance converted onto stator side, (Ω).
For the pump, and linearized by the small deviation, we can derive the simple and practical mathematical model of variable frequency pressure regulation supplying glue system in engineering, and when the voltage and frequency ratio of the power supply of squirrel-cage motors remains unchanged, its small deviation linear model is that: Liu et al. [8,9] considered a linearized relationship between torque and rotational speed, that is, m = K m ω, then, above system is written as following linear system: where ω 0 and ω are the deviations of stator power frequency (Hz) and rotor electric angular velocity (rad/s), respectively; D > 0 is friction coefficient (dimensionless); m is the deviation of load torque (N · m); J > 0 is moment of inertia (Kg · m 2 ); K 0 = p r2 ( V10 ω10 ) 2 , where V 10 and ω 10 are the values on static work point of voltage (V ) and frequency (Hz) of the stator power supply, respectively; p > 0 is pole pair number (dimensionless); r 2 > 0 is rotor resistance converted onto stator side (Ω).
By ignoring the dependence of driving and load torques on the angle of rotation, the corresponding interaction vanishes and so does the effect on the static and transient behaviour of the drive. If we neglect the electrical transients and the dynamics of the load, the remaining mechanical system is described by a first order nonlinear differential equation. A change in speed will result from a change in the voltage applied to the motor or from a change in load conditions when the machine is operating under stable conditions. However, asynchronous machine is a multi-variable, strong coupling and nonlinear controlled object. We also consider the internal transients and interactions in the electrical machine and the load. Actually, the rate of change of angular velocity is nonlinear function associated with angular velocity due to the effect of wind resistance, the change of load, load switching and the change of the liquid's viscosity. Here, the greater angular velocity is, the greater change of resistance will be, and both the angular velocity of motor and gear pump effect load.
In pump load, the flow of pump is related to its speed. Basing on the torque characteristics of pump load, we have following two facts: (1) the resistance of liquid in a certain speed range is roughly proportional to the square of the speed. As the speed decreased, the torque decreases as the quadratic power of the speed, that is, the "quadratic decreasing load torque", m = K 1 ω 2 . Moreover, there exists time delay during the process of pipeline transportation of fluid, thus, Ding [10] studied the Hopf bifurcation and Bogdanov-Takens bifurcation for variable frequency water supply system with delay as follows: where a and b are feedback strengthes, and τ is the time delay of pipeline transportation of glue, (s).
In pump load, the second fact is as follows: (2) Load power p is proportional to the third power of the speed, that is, p = K 2 ω 3 , where K 2 is power constant of load. Note that both load torque and load power can effect the angular velocity, thus, we consider the effects from the two nonlinear terms K 1 ω 2 and K 2 ω 3 on the change of angular velocity. We know that there exists time delay during the process of flow transportation in pipeline. To control the flow is essentially to control the speed of motor, and both torque and power can effect the change of speed. Combining with the facts (1) and (2), namely, torque is proportional to the square of the angular velocity, and load power p is proportional to the third power of the speed, and we consider the above effects as a nonlinear delay function expressed by a cubic polynomial approximatively, that is, the terms aω(t − τ ) + bω 2 (t − τ ) + cω 3 (t − τ ). Then, we obtain the improved model associated with glue flow control system of particleboard glue applying as follows: where a, b and c are effect strengthes, K 0 , K 1 and K 2 are proportionality coefficients, and τ is the time delay of pipeline transportation of glue, (s). Actually, Eq. (5) has two more cubic terms K 2 ω 3 and cω 3 (t − τ ) than Eq. (4).

Existence and stability of equilibria and existence of bifurcations
In this section, we consider system (5) associated with the glue flow control system of particleboard. First of all, we determine the equilibria of above system. Equilibrium solutions of (5) are obtained by simply settinġ ω = 0, and note that ω(t) = ω(t − τ ) for equilibrium solution, thus, the equilibrium of system (5) meets the following equation: Eq. (6) at least has one real root. Let then we have the following lemma associated with the equilibrium of system (5).
Without loss of generality, we assume that E k (k = 1, 2, 3) is the equilibrium of system (5). Transferring the equilibrium E k (k = 1, 2, 3) to the origin, we makê ω = ω−E k , substituting it into system (5), then take off the hat for convenience, we have the following equation: Eq. (8) has trivial equilibrium ω = 0.
We obtain the characteristic equation of the linearized system of Eq. (8) at trivial equilibrium as follows: where When τ = 0, Eq. (9) becomes , and h(ω) is given in (6), then we have following lemma. Lemma 2 Considering system (5).
Case 2: Hopf bifurcation. Next, we consider the existence of local Hopf bifurcation for system (5). Let λ = iβ k (i 2 = −1, β k > 0, k = 1, 2, 3) be a root of Eq. (9). Substituting λ = iβ k into Eq. (9) and separating the real and imaginary parts, we have Adding the square of two equations (10), we can obtain If |B k | ≤ |A k |, system (5) does not occur Hopf bifurcation. When |B k | > |A k |, we can obtain where Then, we deduce the following transversality conditions: Note that J > 0, and combining with Lemmas 3.1 and 3.2, we present a theorem associated with the stability of equilibria and the existence of Hopf bifurcation for system (5). (5), we have following results.
Combining with the results associated with fold bifurcation, we can obtain the following theorem on B-T bifurcation.

Hopf bifurcations
When τ = τ (j) k , k = 1, 2, 3; j = 0, 1, 2, · · · , the characteristic equation (9) has one pair of pure imaginary roots λ = ±iβ k , where τ (j) k is given by (11). We treat the time delay τ as the bifurcation parameter and τ ϵ as perturbation parameter, then system (5) undergoes a Hopf bifurcation near equilibrium E k (k = 1, 2, 3) at the critical point: τ = τ c . Due to the common routine for deriving normal form of Hopf bifurcation, here we only state the normal form results and omit the computation details, and the similar regular process can refer to ref. [10]. The Hopf bifurcation normal form of system (5) derived by using the multiple time scales method [11] is given as follows: where Under the polar coordinates transforming: G = re iΘ , we finally obtain the amplitude and phase equations of (13) on the center manifold for equilibrium E k (k = 1, 2, 3) as follows: Theorem 4 If Re(f k )Re(g k )τ ϵ < 0, the first equation of (14) has nontrivial positive equilibrium The trivial equilibrium of the first equation of (14) remains equilibrium of system (5), while non-trivial equilibrium corresponds periodic solution of system (5). Thus, if Re(f k )Re(g k )τ ϵ < 0 and Re(f k )τ ϵ > 0 (< 0), the bifurcating periodic solutions of system (5) are stable (unstable) on the center manifold.

Bogdanov-Takens bifurcations
When A 2 + B 2 = 0 and B 2 τ + 1 = 0, system (5) undergoes a Bogdanov-Takens bifurcation at , where x, q 1 , f 1 , f 2 are given by (7), Rescaling the time by t → t/τ to normalize the delay, then Eq. (8) yieldṡ Stability analysis of nonlinear glue flow control system with delay 7 (15) Next, we introduce two bifurcation parameters as τ = τ c + µ 1 and a = a c + µ 2 , and the Bogdanov-Takens bifurcation normal form of system (15) is derived by using the center manifold reduction and normal form theory due to Faria and Magalhães [12,13] as follows, where . We also omit the computation details for deriving normal form of system (5) associated with Bogdanov-Takens bifurcation, and the similar regular process can refer to ref. [10].
Note that d 1 d 2 ̸ = 0, then we have Theorem 5 If A 2 + B 2 = 0 and B 2 τ c + 1 = 0, system (5) undergoes a Bogdanov-Takens bifurcation at equilibrium E 2 , and system (15) is locally topologically equivalent near the origin to the system (5) near the equilibrium E 2 .
There exist two essentially different types of bifurcation phenomenon for the normal form (16) near Bogdanov-Takens bifurcation according to the signs of d 1 and d 2 .
Note that under suitable rescaling and the change of variables, for any d 1 d 2 ̸ = 0, Eq. (16) can be reduced to a normal form as follows, with s = ±1. The complete bifurcation phenomena of system (17) near trivial equilibrium can be found in Ref. [14]. Of course, the normal form (17) is not unique, we can also obtain another normal form with s = ±1. The complete bifurcation phenomena of system (18) can be found in Ref. [15]. For example, if d 1 d 2 < 0 in Eq. (16), by introducing the following change of variables and rescaling of time

system (16) becomes the normal form (17) with
In the plane of the parameters u 1 and u 2 , the small neighbourhood of the origin is divided into four regions by three types of bifurcation curves. The detailed explanation can be found in Ref. [14], and here we just briefly list some main results.
Theorem 6 (Ref. [14]) For sufficiently small u 1 , u 2 , Another example, if d 1 d 2 > 0 in Eq. (16), by introducing the change of variables and rescaling of time In the plane of the parameters u 1 and u 2 , the small neighbourhood of the origin is divided into four regions by three types of bifurcation curves. The detailed explanation can be found in Ref. [15], and here we just briefly list some main results. where
Remark 3: In the view of mathematics, there exists a negative equilibrium point E 1 = ω = −1.2166, while it makes sense only for positive value since ω means the angular velocity. We show this example only to verify the results of Theorem 3.2(1), actually, this negative equilibrium point E 1 is unstable, which is not contradict actual process.
Remark 4: System (5) is a scale differential equation. In order to show the periodic solutions clearly, we simulate state variablesω with respective to ω (see Figs. 6-8), and Fig. 8(b) gives the Poincaré map of the chaotic attractor with Poincaré section ω(t − τ ) = 0.  Next, we show a bifurcation diagram of system (5) for τ ∈ [0, 0.0675] (see Fig. 9), which describes the maximum and minimum values of solutions of system (5) as the unfolding parameter τ increases. Figure 9 shows that system (5) undergoes a series of period-doubling bifurcations, then occurs chaotic behavior as the time delay τ increases. Remark 5: Although we can not prove the dynamics phenomenon of period-doubling bifurcations theoretically, we can avoid chaos by limiting time delay according to numerical simulation analysis. We know that the flow of glue tracks the flow of particle in time, and we can adjust the speed of motor in order to control the flow of glue. Time delay is larger means the time of glue transportation is longer, then larger time delay leads to lower time efficiency for controlling the flow of glue, thus, the stability of system is poorer, even lead to chaotic phenomena, and we need to avoid large time delay, namely, limit the length of pipeline. Actually, we can prove the region of delay in which the system is to be stable states (stable equilibrium or stable periodic states), then we can design the proper length of pipeline, that's the reason that we consider time delay as a control parameters.
Example 4: Choosing a = 25, b = 8, p = 2, D = 1, J = 7, K 0 = 4, K 1 = 2, ω 0 = 2, K 2 = 4 27 , c = 4 27 , then we obtain f 1 = 8 27 , f 2 = 4, f 3 = 16, f 4 = 16, q 1 = −6.75, p 1 = −6.75 and ∆ = 0 from (7). Then from Theorem 3.2(2) and Eq. (9), we have E 2 = −6, A 2 = 55 7 and B 2 = − 55 7 . Let τ c = 7 55 ≈ 0.1273, then A 2 + B 2 = 0, 1 + B 2 τ c = 0, and system (5)  = 0} corresponds to a fold bifurcation, along this curve system (17) has an equilibrium with a zero eigenvalue. Entering from region 1 into region 2 through the component S − of the fold curve yields two equilibria: a saddle and a stable node. Then the node turns into a focus and loses stability as we cross the Hopf bifurcation boundary H = {(u 1 , u 2 ) : u 1 = 0, u 2 < 0}. A stable limit cycle is present for close parameter values to the left of H. As we move clockwise, the stable cycle "grows" and approaches the saddle, turning into a homoclinic orbit at Then, an unstable node and a saddle, existing for the parameter values in region 4, collide and disappear at the fold curve S + . The above bifurcation curves are given in Fig. 10, which describes the bifurcation phenomena near the Bogdanov-Takens bifurcation critical point for Example 4.  Hopf bifurcation to an unstable periodic orbit below u 2 = √ −u 1 , u 1 < 0. u 1 = − 49 25 u 2 2 is homoclinic bifurcation curve. Homoclinic bifurcation is described by the periodic orbit created in the Hopf bifurcation growing in amplitude as u 2 is decreased until it collides with the saddle point, creating a homoclinic orbit. As u 2 is further decreased, the homoclinic orbit breaks. The above bifurcation curves are given in Fig. 11, which describes the bifurcation phenomena near the Bogdanov-Takens bifurcation critical point for Example 5. Remark 6: Under the same group of parameters, a = 25, b = 8, p = 2, D = 1, J = 7, K 0 = 4, K 1 = 2, ω 0 = 2 in Examples 4 and 5, there exist different dynamics phenomena near Bogdanov-Takens bifurcation critical point for different cubic terms coefficients K 2 and c. That is, there exist stable periodic solutions and stable equilibrium near Bogdanov-Takens bifurcation critical point in Example 4 due to d 1 d 2 < 0 (see Fig. 10), while does not exist stable periodic solutions near Bogdanov-Takens bifurcation critical point in Example 5 due to d 1 d 2 > 0 (see Fig. 11).
Remark 7: Eq. (5) describes the glue flow control problem associated with glue applying of particleboard. ω describes the electrical angular velocity of rotor, which is preferred to vary in a stable state (stable equilibrium or stable periodic solution). Stable electrical angular velocity generates stable control signal, then leads to stable flow rate of glue. We study the critical conditions under which system undergoes some types of bifurcation (see Theorem 3.2 and Theorems 4. 1-4.4), and analysis the dynamics phenomena near above bifurcation critical values. Thus, the flow control system (5) can be controlled into some types of stable states by adjusting control parameters, and we give the regions of parameters under which system (5) exists stable equilibrium solutions or stable periodic solutions. The method presented in this paper may provide a very good tool in the study of real problems.
6 Comparison of Eqs. (4) and (5) Eq. (4) describes nonlinear variable frequency water supply system associated with glue applying of particleboard, and Eq. (5) describes glue flow control system associated with glue applying of particleboard. Ref. [10] studied the dynamical behaviors of Eq. (4), while in this paper, we consider Eq. (5). Although the two models illustrate different questions associated with glue applying of particleboard, Eq. (5) has only two more cubic terms K 2 ω 3 and cω 3 (t − τ ) than Eq. (4) looking at expressions of the two equations. Therefore, we compare their dynamical properties on the view of mathematics in this section.
Firstly, both of Eq. (4) and Eq. (5) can undergo fold bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation, and both of their periodic solutions can exist in some region of parameter, and then occur chaotic phenomena from a sequence of period-doubling bifurcations. These similar results are quite normal, after all, the two systems have most same terms and forms.
Secondly, under the group of parameters of Example 1 in Section 5, that is, a = 100, b = 3, p = 2, D = 0.8, J = 1, ω 0 = 5, K 0 = 10, K 1 = 2, K 2 = 0.2, c = 1, system (5) only has one equilibrium E 1 = −1.2166, which is unstable for ∀τ ≥ 0. While also under this group of parameters with K 2 = 0, c = 0, that is, we do not consider the effect of cubic terms, then Eq. (4) has two equilibria, one is local asymptotically stable for ∀τ ≥ 0, and the other one is unstable (see the Example 1 of numerical simulations in Ref. [10], note that K m in Ref. [10] corresponds to K 1 in this paper). Obviously, the cubic terms of system (5) may effect the number of equilibria.
Remark 8: Eq. (4) and Eq. (5) describe nonlinear variable frequency water supply system and glue flow control system associated with glue applying of particleboard, respectively. Eq. (5) has only two more cubic  Fig. 13 The Bogdanov-Takens bifurcation lines of system (5).
terms K 2 ω 3 and cω 3 (t − τ ) than Eq. (4), and the cubic terms of system (5) not only may change the number of equilibria, but also may effect the dynamical properties near Bogdanov-Takens bifurcation critical point.

Conclusion
In this paper, we have constructed a glue flow control system associated with glue applying of particleborad, and considered the dynamical properties of this system. We have shown the existence conditions of some types of equilibria, analysed the stability of equilibria and the existence of several types of bifurcations. Further, we have derived the normal forms of Hopf bifurcation and Bogdanov-Takens bifurcation, and shown the bifurcation analysis. Five numerical simulation examples are carried out to demonstrate the applications of the theoretical results. Next, under proper parameters, we have shown that system (5) may exist stable equilibrium, stable periodic-1, period-2 or period-4 solutions, and a complex chaotic attractor from a sequence of perioddoubling bifurcations. Finally, we compare the dynamical phenomena of Eqs. (4) and (5), showing that cubic terms in Eq. (5) can effect its dynamical properties.

Conflict of Interest
The authors declare that they have no conflict of interest.
Availability of data The datasets used or analyzed during the current study are available from the corresponding author on reasonable request.
Funding This study was funded by the Heilongjiang Provincial Natural Science Foundation of China (Grant No. LH2019A001).