Analysis of Underlapped Symmetrically Ported Valve Controlled Asymmetric Cylinder Drive

1. Introduction

Hydraulic power systems [3] are widely used in industrial manufacture, building construction, port transportation, funfair, etc. They use fluid to achieve energy transmission without a rigid joint. Hydraulic power is able to produce massive force and torque within a compact and lighter package, and the conventional rigid components like gears and joints can be avoided.

Hydraulic linear actuators are mainly symmetric cylinder or asymmetric cylinder. The symmetric cylinder has equal piston area and the asymmetric type (also named as single rod) cylinder comes with different piston areas. The symmetric cylinder is usually utilised for high-performance requirements, due to its linear performance with symmetric ported control valve. This type of cylinder is widely applied in aerospace and military for high performance and dynamic response [4]. However, the asymmetric cylinder [5] drive occupies 80% of the hydraulic applications [6], such as cranes, injection moulding systems, and hydraulic press machines, etc. The asymmetric cylinder drive is usually controlled by a symmetric ported valve, unequal flow rate will be delivered into cylinder chambers with same valve opening in the extending or retracting state, which means system parameters do not remain the same after the valve crosses its null position [7]. Leaney [8] did research in these nonlinearities, but how the system performs when the valve in the underlap region is not explicitly explained. Viersma [1] indicates that pressure jumps will occur when the valve spool is around its null position, causing ‘implosion’ or ‘explosion’ of oil owing to its compressibility, so that smooth operation around valve null position is impossible.

A test rig of a four-way servovalve controlled asymmetric cylinder hydraulic system is constructed by Leaney [2], the system is a typical symmetric ported valve controlled asymmetric cylinder system, which is depicted in Figure 1, the Moog series 76 model 102 four-way valve is used to control the asymmetric cylinder system to achieve the desired output.

Leaney [8] simulated his test rig with ‘Component Linking’ method, which models each of the components such as valve, actuator/cylinder and load, then combines them to perform the whole system simulation. This method aims at simplifying the ‘Small Perturbation’ model [9], which utilised a single transfer function to describe a whole hydraulic system. However, as not all parameters of a hydraulic system remain constant during operation and a single transfer function is mainly designed for a linear time-invariant (LTI) system, a single transfer function [10] is difficult to capture the behaviours of a hydraulic system. Modelling the components separately is easier for the analysis and adjusting of the system parameters.

The components of the asymmetric cylinder system are connected by the concept of ‘power bonding’, in other words, the components are connected by the power transfer, such concept is still utilised nowadays for a hydraulic system of the 50 wheel loader [11]. In Leaney’s paper [2], three components are the major concerns, which are valve, actuator and load model. Their power transfers are indicated in Figure 2, the voltage is input to the valve model to control the valve opening, the flow passes through the valve and drops to a certain value. Pressure and flow are delivered to the actuator model, which outputs the force/torque and velocity. Finally, the load model outputs the velocity and force.

1.1. Valve model

The control valve (electrohydraulic valve) is a Moog 76 model 102 flow control four-way servovalve, its schematic diagram is depicted in Figure 3. If current is applied in the armature coil in the torque motor, the T-shape flapper rotates slightly. Assuming the flapper rotates clockwise, the offset between the left nozzle and flapper decreases, which leads to the decrease of flow passing through the left nozzle and causes the pressure increasing in the left control line. A similar process occurs in the right control line, leading to a decrease of pressure.

The pressure difference in the control line moves the spool, therefore bends the feedback spring to balance the torque generated by the torque motor, so that the spool displacement is dependent on the rotation of the torque motor; furthermore, the spool displacement is dependent on the current.

Appropriate transfer functions, equations 1 and 2, are given by Moog [13]. The expressions are linear, empirical relationships that approximate the responses of the actual servovalve. For the frequency range below 50 Hz, the valve dynamic can be approximated by a first-order transfer function below [13]:

$$\frac{\mathrm{Q}\left(\mathrm{s}\right)}{\mathrm{I}\left(\mathrm{s}\right)}=\mathrm{K}\frac{1}{1+{\mathsf{\tau }}_{\mathrm{s}}\mathrm{s}}$$

(1)

where $\mathrm{K}$ is the flow gain $\left[\raisebox{1ex}{${\mathrm{m}}^3$}\!\left/ \!\raisebox{-1ex}{$\mathrm{s}·\mathrm{a}\mathrm{m}\mathrm{p}$}\right.\right]$ and ${\mathsf{\tau }}_{\mathrm{s}}$ is the servovalve time constant $\left[\mathrm{s}\mathrm{e}\mathrm{c}\right]$.

A second-order transfer function is necessary to represent wider frequency range response as given below [13]:

$$\frac{\mathrm{Q}\left(\mathrm{s}\right)}{\mathrm{I}\left(\mathrm{s}\right)}=\mathrm{K}\frac{1}{1+\left(\frac{2\mathsf{\zeta }}{{\mathsf{\omega }}_{\mathrm{n}}}\right)\mathrm{s}+{\left(\frac{\mathrm{s}}{{\mathsf{\omega }}_{\mathrm{n}}}\right)}^2}$$

(2)

where ${\mathsf{\omega }}_{\mathrm{n}}$ is the natural frequency [$\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}$], and $\mathsf{\zeta }$ is the damping ratio. In a limit cycle oscillation study of a low-frequency phenomenon, a higher order valve model does not have significant difference compared to a lower order model [14]. Hisa [15] stated that a reduced order transfer function can be used to approximate a higher order transfer function, and a first order model is used as a reference model by Persson et al [16].

As all of the valve operational frequencies in this research are much lower than 50 Hz, the first order transfer function equation 1 is chosen for the valve model. As the displacement of the spool, $\mathrm{X}\left(\mathrm{s}\right)$, controls the flow rate of the valve, equation 1 can be modified as:

$$\frac{X_s\left(s\right)}{I\left(s\right)}=K_s\frac{1}{1+{\tau }_{s}s}$$

(3)

where ${\mathrm{K}}_{\mathrm{s}}$ is the gain of the system [$\raisebox{1ex}{$\mathrm{m}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{a}\mathrm{m}\mathrm{p}$}\right.$].

The current is generated by an amplifier. Assuming the amplifier is ideal, then current is proportional to the input voltage as given by equation 4:

$$\frac{I\left(s\right)}{V_v\left(s\right)}=K_A$$

(4)

where ${\mathrm{K}}_{\mathrm{A}}$ is the gain of voltage to current [$\raisebox{1ex}{$\mathrm{a}\mathrm{m}\mathrm{p}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{t}$}\right.$]. Hence the valve dynamics are represented by:

$$\frac{X_s\left(s\right)}{V_v\left(s\right)}=\frac{K_A{K}_s}{1+{\tau }_{s}s}$$

(5)

A block diagram for this valve model is depicted in Figure 4.

1.2. Flow gain

The flow gain is the gain that depicts the relationship between the valve opening and the oil flow across it. In order to interpret it clearly, a physical layout of the four-way valve-controlled [17] asymmetric cylinder is listed in Figure 5:

${\mathrm{Q}}_1$ is the flow rate that delivered into the piston side chamber, ${\mathrm{V}}_1$ is the volume of fluid that trapped in the piston chamber, ${\mathrm{A}}_1$ is the piston side area of the cylinder, the subscript 2 indicates corresponding parameters in rod side chamber. ${\mathrm{P}}_{\mathrm{s}}$ is the supply pressure, ${\mathrm{P}}_{\mathrm{T}}$ is the return line pressure of oil tank, which is zero. ${\mathrm{A}}_{\mathrm{V}}$ is the opening area of the corresponding port of the valve, subscripts 1, 2, 3 and 4 indicate different ports of the valve which are depicted in Figure 5. ${\mathrm{x}}_{\mathrm{P}}$ is the cylinder displacement. ${\mathrm{Q}}_{\mathrm{P}1}$ is the flow pass through the ${\mathrm{A}}_{\mathrm{V}1}$ port, ${\mathrm{Q}}_{1\mathrm{T}}$ is the flow pass through ${\mathrm{A}}_{\mathrm{V}2}$ port, ${\mathrm{Q}}_{2\mathrm{T}}$ is the flow pass through the ${\mathrm{A}}_{\mathrm{V}3}$ port and ${\mathrm{Q}}_{\mathrm{P}2}$ is the flow pass though the ${\mathrm{A}}_{\mathrm{V}4}$ port.

The flow deliver into the chambers can be described as:

$$Q_1=Q_{P1}-Q_{1T}=C_D\sqrt{\frac{2}{\rho }}\left[A_{V1}{\left(P_s-P_1\right)}^{\frac{1}{2}}-A_{V2}{\left(P_1\right)}^{\frac{1}{2}}\right]$$

(6)

$$Q_2=Q_{2T}-Q_{P2}=C_D\sqrt{\frac{2}{\rho }}\left[A_{V3}{\left(P_2\right)}^{\frac{1}{2}}-A_{V4}{\left(P_s-P_2\right)}^{\frac{1}{2}}\right]$$

(7)

For a symmetric ported valve, ${\mathrm{A}}_{\mathrm{V}1}={\mathrm{A}}_{\mathrm{V}3}$ and ${\mathrm{A}}_{\mathrm{V}2}={\mathrm{A}}_{\mathrm{V}4}$, so that equations 6 and 7 can be simplified as:

$$Q_1=Q_{P1}-Q_{1T}=C_D\sqrt{\frac{2}{\rho }}\left[A_{V1}{\left(P_s-P_1\right)}^{\frac{1}{2}}-A_{V2}{\left(P_1\right)}^{\frac{1}{2}}\right]$$

(8)

$$Q_2=Q_{2T}-Q_{P2}=C_D\sqrt{\frac{2}{\rho }}\left[A_{V1}{\left(P_2\right)}^{\frac{1}{2}}-A_{V2}{\left(P_s-P_2\right)}^{\frac{1}{2}}\right]$$

(9)

These two equations can be applied when the valve is in the extending, retracting and underlap region, the underlap region is the valve spool position where all ports opening areas are non-zero. Even for a zerolapped valve, due to the manufacture tolerance [18] , there will be a small underlap when the spool is placed in the null position. Figure 6 depicts the theoretical and actual situation for a nominally zero lapped valve.

The manufacture tolerances can be identified as:

i Radial clearance between spool and sleeve

ii Tolerance on lap

Such tolerance [19] will lead to variation of flow gain when the valve in its null region. The flow gain model block diagram is depicted in Figure 7.

1.3. The actuator model

The actuator model connects the valve and load model, which accepts the information of flow and pressure as input, then output force to the load model. The types of hydraulic actuators are mainly symmetric cylinder, asymmetric cylinder and hydraulic motor. Due to the target of this research, the asymmetric cylinder is the only one under analysis.

No matter the asymmetric cylinder in the extending or retracting state, the flow in its chamber consists of input flow, displacement flow and compressibility flow. But under different operational states, the flow balance equations in different chambers are not the same.

For the cylinder extending state:

$$Q_1=Q_{D1}+Q_{C1}$$

(10)

$$Q_2=Q_{D2}-Q_{C2}$$

(11)

For the cylinder retracting state:

$$Q_2=Q_{D2}+Q_{C2}$$

(12)

$$Q_1=Q_{D1}-Q_{C1}$$

(13)

where

$$Q_{D1}=v·A_1 Q_{D2}=v·A_2$$

(14)

$$Q_{C1}=\frac{s{P}_1·V_1}{B} Q_{C2}=\frac{s{P}_2·V_2}{B}$$

(15)

$\mathrm{A}$ is the area of cylinder piston side area or rod side area [${\mathrm{m}}^2$], $\mathrm{v}$ is the cylinder velocity [$\raisebox{1ex}{$\mathrm{m}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{s}$}\right.$], $\mathrm{V}$ is the cylinder chamber volume [${\mathrm{m}}^3$], $\mathrm{P}$ is the cylinder chamber pressure [$\raisebox{1ex}{$\mathrm{N}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^2$}\right.$] and $\mathrm{B}$ is the fluid Bulking Modulus [$\raisebox{1ex}{$\mathrm{N}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^2$}\right.$]. The subscript $1$ and $2$ indicate the piston side chamber and rod side chamber respectively.

A generalised definition of these equations is proposed by Leaney [2], which uses effective parameters to describe the system behaviours rather analyse them separately in different cylinder chambers. The generalised equations are depicted as:

$$\begin{gathered} P_{LE}=\frac{P_1{A}_1-P_2{A}_2}{A_E}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{Q}_{LE}=\frac{Q_1+Q_2}{2} \\ {A}_E=\frac{A_1+A_2}{2}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{Q}_{DE}=\frac{Q_{D1}+Q_{D2}}{2} \\ {Q}_{CE}=\frac{Q_{C1}+Q_{C2}}{2}=Q_{LE}-Q_{DE} \end{gathered}$$

(16)

Substitute these effective parameters into the flow and pressure equations reveals:

For cylinder extending (${\mathrm{A}}_{\mathrm{V}2}$=0 and ${\mathrm{A}}_{\mathrm{V}1}$>0), the service line pressures can be rewritten as [2]:

$${\left(P_1\right)}_{ext}=\frac{P_s+{\mathsf{\gamma }}^2\frac{1}{2}\left(1+\mathsf{\gamma }\right)P_{LE}}{1+{\mathsf{\gamma }}^3}$$

(17)

$${\left(P_2\right)}_{ext}=\frac{{\mathsf{\gamma }P}_s-\frac{1}{2}\left(1+\mathsf{\gamma }\right)P_{LE}}{1+{\mathsf{\gamma }}^3}$$

(18)

For cylinder retracting (${\mathrm{A}}_{\mathrm{V}2}$>0 and ${\mathrm{A}}_{\mathrm{V}1}$=0), the service line pressures can be rewritten as [2]:

$${\left(P_1\right)}_{ret}=\frac{{{\mathsf{\gamma }}^{2}P}_s+{\mathsf{\gamma }}^2\frac{1}{2}\left(1+\mathsf{\gamma }\right)P_{LE}}{1+{\mathsf{\gamma }}^3}$$

(19)

$${\left(P_2\right)}_{ret}=\frac{{{\mathsf{\gamma }}^{3}P}_s-\frac{1}{2}\left(1+\mathsf{\gamma }\right)P_{LE}}{1+{\mathsf{\gamma }}^3}$$

(20)

where $\mathsf{\gamma }$ is the area ratio of the cylinder, Leaney [2] introduced $\mathrm{X}$ factor to unify the expression for ${\mathrm{P}}_1$, ${\mathrm{P}}_2$ and utilised the generalised concepts in the component linking method. The $\mathrm{X}$ factor value depends on the region of the spool operation, $\mathrm{X}$=1 for cylinder extending state and $\mathrm{X}$=${\mathsf{\gamma }}^2$ for cylinder retracting state. And the service line pressure ${\mathrm{P}}_1$ and ${\mathrm{P}}_2$ can be written as:

$$P_1=\frac{{XP}_s+{\mathsf{\gamma }}^2\frac{1}{2}\left(1+\mathsf{\gamma }\right)P_{LE}}{1+{\mathsf{\gamma }}^3}$$

(21)

$$P_2=\frac{{\mathrm{X}\mathsf{\gamma }P}_s-\frac{1}{2}\left(1+\mathsf{\gamma }\right)P_{LE}}{1+{\mathsf{\gamma }}^3}$$

(22)

$$P_{LE}=\left(\frac{4B}{s{V}_{TE}}\right)Q_{CE}$$

(23)

where $\mathsf{\gamma }$ is the area ratio and ${\mathrm{V}}_{\mathrm{T}\mathrm{E}}$ is effective total volume [${\mathrm{m}}^3$].

$$V_{TE}=\frac{2A_E}{A_1}\left[\frac{V_1+{\left(\frac{A_2}{A_1}\right)}^2{V}_2}{1+{\left(\frac{A_2}{A_1}\right)}^3}\right]$$

(24)

The $\mathrm{X}$ factor value is known and clear when the valve is outside the underlap region. However, when the valve spool is in the underlap region, the $\mathrm{X}$ value is difficult to be found explicitly and it can be numerically solved by an implicit equation 25 [2] as below:

$$A_{v1}{\left(1+{\mathsf{\gamma }}^3-X\right)}^{\frac{1}{2}}-A_{v2}{\left(X\right)}^{\frac{1}{2}}=\mathsf{\gamma }\left[A_{v1}{\left(X\mathsf{\gamma }\right)}^{\frac{1}{2}}-A_{v2}{\left(1+{\mathsf{\gamma }}^3-X\mathsf{\gamma }\right)}^{\frac{1}{2}}\right]$$

(25)

Now the effective load flow can be revealed [2]:

$$\begin{gathered} Q_{LE}=\frac{C_D}{\sqrt{2\rho }}\sqrt{\frac{{\left(1+\gamma \right)}^2}{1+{\gamma }^3}}\left[A_{v1}{\left(X\gamma P_s-\frac{1}{2}\left(1+\gamma \right)P_{LE}\right)}^{\frac{1}{2}} \\ -A_{v2}{\left[\left(1+{\gamma }^3-X\gamma \right)P_s+\frac{1}{2}\left(1+\gamma \right)P_{LE}\right]}^{\frac{1}{2}}\right] \end{gathered}$$

(26)

The actuator model block diagram is depicted in Figure 8.

1.4. The load model

The load model is various based on different test rig set up, in order to simplify the analysis afterwards, the load set up is a simple asymmetric cylinder coupling with a load, which is placed on a surface, So that the load force balance can be depicted by:

$$P_1{A}_1-P_2{A}_2=P_{LE}{A}_E=m·a+friction force+F_{external}$$

(27)

where $\mathrm{m}$ is the load mass [$\mathrm{k}\mathrm{g}$] and $\mathrm{a}$ is the acceleration [$\mathrm{m}/{\mathrm{s}}^2$].

2. Materials and Methods

Leaney [2] stated that the cylinder chamber pressures ${\mathrm{P}}_1$ and ${\mathrm{P}}_2$ must be numerically solved when the valve in the underlap region, due to the derived implicit equation 25. The chamber pressures are the power source to drive the cylinder, so that the system behaviours are decided by chamber pressures. How the pressures perform when the valve is outside the underlap region is explicit. The chamber pressures behaviours in underlap region are solved numerically as shown in Figure 10.

An input signal is set to drive the valve moving across the underlap region in a small-time interval with a constant speed. The supply pressure is 70 bar and all the modelling parameters are the same as those in Moog valve controlled asymmetric cylinder system in the last section.

As the chamber pressures of the asymmetric cylinder when the valve in the underlap region is modelled and they are numerically solved. The key part is the $\mathrm{X}$ factor equation 25 from [2]. The solving command in Matlab is fsolve and it uses the rust-region algorithm [14] to find the root of factor $\mathrm{X}$ with given valve opening areas ${\mathrm{A}}_{\mathrm{v}1}$, ${\mathrm{A}}_{\mathrm{v}2}$ and area ratio $\mathsf{\gamma }$. Numerically solving equation 25 takes a lot of computing time, if the $\mathrm{X}$ value can be analytical solved, simulation time cost can be massively reduced. A trial is given in the rest of this section.

Square both sides of equation 25 leads to:

$$\begin{gathered} \frac{{A_{V1}}^2}{{\gamma }^2}\left(1+{\gamma }^3-X\right)+\frac{{A_{V2}}^{2}X}{{\gamma }^2}-\frac{2A_{V1}{A}_{V2}}{{\gamma }^2}{\left(X+{\gamma }^{3}X-X^2\right)}^{\frac{1}{2}} \\ ={A_{V1}}^{2}X\gamma +{A_{V2}}^2\left(1+{\gamma }^3-X\gamma \right)-2A_{V1}{A}_{V2}{\left(X\gamma +X{\gamma }^4-X^2{\gamma }^2\right)}^{\frac{1}{2}} \end{gathered}$$

(28)

Rearrange it:

$$\begin{gathered} \left[\frac{A_{V1}}{2A_{V2}}\left(1+{\gamma }^3\right)-\frac{A_{V2}}{2A_{V1}}\left({\gamma }^2+{\gamma }^5\right)\right]-\left[\frac{A_{V1}}{2A_{V2}}\left(1+{\gamma }^3\right)-\frac{A_{V2}}{2A_{V1}}\left(1+{\gamma }^3\right)\right]X \\ ={\left[\left(1+{\gamma }^3\right)X-X^2\right]}^{\frac{1}{2}}-{\gamma }^2{\left[\left(\gamma +{\gamma }^4\right)X-{{\gamma }^{2}X}^2\right]}^{\frac{1}{2}} \end{gathered}$$

(29)

Let $\mathrm{m}$ be given as:

$$m=\frac{A_{V1}}{2A_{V2}}\left(1+{\gamma }^3\right)-\frac{A_{V2}}{2A_{V1}}\left({\gamma }^2+{\gamma }^5\right)=\left(1+{\gamma }^3\right)\left(\frac{A_{V1}}{2A_{V2}}-\frac{A_{V2}{\gamma }^2}{2A_{V1}}\right)$$

(30)

Let $\mathrm{n}$ be given as:

$$n=\frac{A_{V1}}{2A_{V2}}\left(1+{\gamma }^3\right)-\frac{A_{V2}}{2A_{V1}}\left(1+{\gamma }^3\right)=\left(1+{\gamma }^3\right)\left(\frac{A_{V1}}{2A_{V2}}-\frac{A_{V2}}{2A_{V1}}\right)$$

(31)

Substitute $\mathrm{m}$ and $\mathrm{n}$ into equation 31 reveals:

$$m-nX={\left[\left(1+{\gamma }^3\right)X-X^2\right]}^{\frac{1}{2}}-{\gamma }^2{\left[\left(\gamma +{\gamma }^4\right)X-{{\gamma }^{2}X}^2\right]}^{\frac{1}{2}}$$

(32)

Square both sides of equation 32 leads to:

$$\begin{gathered} m^2+n^2{X}^2-2mnX \\ =\left(1+{\gamma }^3\right)X-X^2+{\gamma }^5\left(1+{\gamma }^3\right)X-{\gamma }^6{X}^2 \\ -2{\gamma }^2{\left[\gamma {\left(1+{\gamma }^3\right)}^2{X}^2-\gamma \left(1+{\gamma }^3\right)X^3-{\gamma }^2\left(1+{\gamma }^3\right)X^3+{{\gamma }^{2}X}^4\right]}^{\frac{1}{2}} \end{gathered}$$

(33)

Rearrange the equation 33 reveals:

$$\begin{gathered} \left(n^2+1+{\gamma }^6\right)X^2-\left[2mn+\left(1+{\gamma }^5\right)\left(1+{\gamma }^3\right)\right]X+m^2 \\ =2{\gamma }^2{\left[\gamma {\left(1+{\gamma }^3\right)}^2{X}^2-\gamma \left(1+{\gamma }^3\right)X^3-{\gamma }^2\left(1+{\gamma }^3\right)X^3+{{\gamma }^{2}X}^4\right]}^{\frac{1}{2}} \end{gathered}$$

(34)

Let

$$\begin{gathered} p=n^2+1+{\gamma }^6 \\ q=2mn+\left(1+{\gamma }^5\right)\left(1+{\gamma }^3\right) \\ h=m^2 \end{gathered}$$

(35)

Equation 34 becomes:

$$p{X}^2-qX+h=2{\gamma }^2{\left[\gamma {\left(1+{\gamma }^3\right)}^2{X}^2-\gamma \left(1+{\gamma }^3\right)X^3-{\gamma }^2\left(1+{\gamma }^3\right)X^3+{{\gamma }^{2}X}^4\right]}^{\frac{1}{2}}$$

(36)

Square both side of the equation 36, reveals:

$$\begin{gathered} p^2{X}^4+q^2{X}^2+h^2-pq{X}^3+ph{X}^2-pq{X}^3-qhX+ph{X}^2-qhX \\ =4{\gamma }^6{X}^4-4{\gamma }^5\left(1+\gamma \right)\left(1+{\gamma }^3\right)X^3+4{\gamma }^5{\left(1+{\gamma }^3\right)}^2{X}^2 \end{gathered}$$

(37)

Rearrange above equation reveals:

$$\begin{gathered} \left(p^2-4{\gamma }^6\right)X^4+\left[4{\gamma }^5\left(1+\gamma \right)\left(1+{\gamma }^3\right)-2pq\right]X^3+\left[\left(2ph+q^2\right)-4{\gamma }^5{\left(1+{\gamma }^3\right)}^2\right]X^2 \\ -2qhX+h^2=0 \end{gathered}$$

(38)

Let

$$\begin{gathered} a_4=p^2-4{\gamma }^6 \\ {a}_3=4{\gamma }^5\left(1+\gamma \right)\left(1+{\gamma }^3\right)-2pq \\ {a}_2=\left(2ph+q^2\right)-4{\gamma }^5{\left(1+{\gamma }^3\right)}^2 \\ {a}_1=-2qh \\ {a}_0=h^2 \end{gathered}$$

(39)

Thus equation 38 becomes:

$$a_4{X}^4+a_3{X}^3+a_2{X}^2+a_{1}X+a_0=0$$

(40)

Let

$$\begin{gathered} a=\frac{a_3}{a_4}, b=\frac{a_2}{a_4} \\ c=\frac{a_1}{a_4}, d=\frac{a_0}{a_4} \end{gathered}$$

(41)

As a result, equation 40 becomes:

$$X^4+a{X}^3+b{X}^2+cX+d=0$$

(42)

The value of $\mathrm{X}$ factor is a problem of find the root of a quartic equation 42. A general form of roots can be depicted as below equations 43:

$$\left\{\begin{array}{c}{X}_1=-\frac{a}{4}+\frac{{\left(U+V\right)}^{\frac{1}{2}}}{2}-\frac{{\left(2U-V-W\right)}^{\frac{1}{2}}}{2}\\ X_2=-\frac{a}{4}-\frac{{\left(U+V\right)}^{\frac{1}{2}}}{2}-\frac{{\left(2U-V+W\right)}^{\frac{1}{2}}}{2}\\ \begin{array}{c}{X}_3=-\frac{a}{4}+\frac{{\left(U+V\right)}^{\frac{1}{2}}}{2}+\frac{{\left(2U-V-W\right)}^{\frac{1}{2}}}{2}\\ X_4=-\frac{a}{4}-\frac{{\left(U+V\right)}^{\frac{1}{2}}}{2}+\frac{{\left(2U-V+W\right)}^{\frac{1}{2}}}{2}\end{array}\end{array}\right.$$

(43)

where

$$\begin{gathered} U=\frac{a^2}{4}-\frac{2b}{3} \\ V=\frac{L}{3}+\frac{Z_0}{3L} \\ W=\frac{\left(a^3-4ab+8c\right)}{4{\left(U+V\right)}^{\frac{1}{2}}} \end{gathered}$$

(44)

where

$$\begin{gathered} Z_0=b^2-3ac+12d \\ {Z}_1=2b^3-9abc+27a^{2}d+27c^2-72bd \\ L={\left[\frac{{Z_1+\left({Z_1}^2-4{Z_0}^3\right)}^{\frac{1}{2}}}{2}\right]}^{\frac{1}{3}} \end{gathered}$$

(45)

Actually, there is only one value of factor with a set of parameters. However, there are four roots from a quartic equation and the value of is one of them.

To identify the correct root, the method is depicted as following steps:

1)Numerically solving the factor in underlap region.

2)Identify one of the roots in a certain valve displacement in underlap region corresponding to the same value from numerically solving.

3)Repeat step 1 and 2 with different area ratio .

4)Repeat step 1, 2 and 3 with different underlap region.

5)Estimate the final analytical solution of factor .

Start with the parameter from Moog valve controlled asymmetric cylinder drive system, for the area ratio = 11.43/5.8=1.9707 and underlap region is from 0.0045mm to -0.0045mm, root identification is listed in Table 1.

A new phenomenon can be noticed, the $\mathrm{X}$ factor values in the underlap region contains two roots of equation 25 ${\mathrm{X}}_1$ and ${\mathrm{X}}_4$. When valve displacement is from -0.0045mm to -0.00075mm, the value of $\mathrm{X}$ factor is the root ${\mathrm{X}}_4$, when valve displacement is from -0.00074mm to 0.0045mm, the value of $\mathrm{X}$ factor is the root ${\mathrm{X}}_1$.

All the four roots values in underlap region are described as in Figure 11(a), different line styles are corresponding to different roots. The values of the roots in the underlap are compared with the numerically solved $\mathrm{X}$ values as in Figure 11(b). It can be noticed that a part of ${\mathrm{X}}_4$ and a part of ${\mathrm{X}}_1$ are merged into the numerically solved curve, they meet at some specific point.

Another interesting phenomenon can be found that some parts of the root curve can be combined with parts of other root curves to generate a smoothed curve. The point where ${\mathrm{X}}_4$ and ${\mathrm{X}}_1$ curves meet may change with some parameters. Change the form of valve opening area ${\mathrm{A}}_{\mathrm{v}1}$ and ${\mathrm{A}}_{\mathrm{v}2}$ in equation 25 into below equation.

$$\begin{gathered} \pi D(x_s+0.0000045){\left(1+{\mathsf{\gamma }}^3-X\right)}^{\frac{1}{2}}-\pi D(x_s-0.0000045){\left(X\right)}^{\frac{1}{2}} \\ =\mathsf{\gamma }\left[\pi D(x_s+0.0000045){\left(X\mathsf{\gamma }\right)}^{\frac{1}{2}}-\pi D(x_s \\ -0.0000045){\left(1+{\mathsf{\gamma }}^3-X\mathsf{\gamma }\right)}^{\frac{1}{2}}\right] \end{gathered}$$

(46)

where ${\mathrm{A}}_{\mathrm{v}1}=\mathsf{\pi }\mathrm{D}({\mathrm{x}}_{\mathrm{s}}+0.0000045)$, ${\mathrm{A}}_{\mathrm{v}2}=\mathsf{\pi }\mathrm{D}({\mathrm{x}}_{\mathrm{s}}-0.0000045)$ [${\mathrm{m}}^2$]

Cancel the $\mathsf{\pi }\mathrm{D}$ from the above equation, and divide both sides by $0.0000045$ reveals

$$\begin{gathered} \left(1+\frac{x_s}{0.0000045}\right){\left(1+{\mathsf{\gamma }}^3-X\right)}^{\frac{1}{2}}-\left(1-\frac{x_s}{0.0000045}\right){\left(X\right)}^{\frac{1}{2}} \\ =\mathsf{\gamma }\left[\left(1+\frac{x_s}{0.0000045}\right){\left(X\mathsf{\gamma }\right)}^{\frac{1}{2}}-\left(1-\frac{x_s}{0.0000045}\right){\left(1+{\mathsf{\gamma }}^3-X\mathsf{\gamma }\right)}^{\frac{1}{2}}\right] \end{gathered}$$

(47)

The $\frac{{\mathrm{x}}_{\mathrm{s}}}{0.0000045}$ in equation 47 can be regard as normalized spool travel from -1 to +1 in underlap. It can be observed that for a certain area ratio $\mathsf{\gamma }$, the factor $\mathrm{X}$ is only affected by the travel percentage. An assumption can be raised that the intersection point in Figure 11(a) is affected by the area ratio, different area ratios $\mathsf{\gamma }$ are applied into equation 47 to verify it. Runs the simulation with different area ratios $\mathsf{\gamma }$ in underlap region, the root ${\mathrm{X}}_4$ and ${\mathrm{X}}_1$ distribution is depicted in Table 2.

Input these data in the Matlab CFTOOL reveals the relationship for ${\mathrm{X}}_1$ percentage with different area ratio $\mathsf{\gamma }$ as equation 48.

$$X_1 Percentage=0.00999{\gamma }^3-0.09352{\gamma }^2+0.3228\gamma +0.2328$$

(48)

The percentage of root ${\mathrm{X}}_4$ is simply obtained by $1-{\mathrm{X}}_1\mathrm{P}\mathrm{e}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{g}\mathrm{e}$, so that the $\mathrm{X}$ factor roots distribution when the valve in the underlap region is obtained.

The overall component linking model block diagram for this Moog valve-controlled asymmetric cylinder drive is depicted in Figure 12.

The overall block model of the Moog valve controlled asymmetric cylinder drive system is simulated in Matlab Simulink as in Figure 13.

The Simulink model with analytical solution simulation results are compared with the numerical solution simulation results in Fortran and the experimental test results in Figure 14 and Figure 15.

The modelling results from Simulink shows consistency with the modelling results from Fortran and experimental results [2], which indicates that the analytical solving method is validated.

3. Results and discussion

The analytical solution for factor in the component linking method is able to save computing time during the simulation. The valve stroke keeps changing during operation, which indicates the computer must solving the equation 25 with different values of . If the step size of the simulation is small, for instance , the computer will need to solve thousands of implicit equations in one second, this process will consume a lot of computing power. This analytical solution in the second section improves this situation.

Simulation tests of the Moog valve controlled asymmetric cylinder drive system are carried out to compare the time consumes between the implicit solving and analytical solving method, including a square wave test, a sawtooth test and a sine wave test.

There is no observable difference between the two results, which indicates both solutions have the same performance. The numerical solve solution uses 5.49 seconds, while the analytical solution uses 1.07 seconds, the numerical solution consumes nearly five times time to finish the simulation. However, the valve in square wave test in Figure 16 is operating outside the underlap region most of the time, the advantage of the analytical solution is not highlighted here.

The sawtooth command test is to compare the time cost between the numerical solution and the analytical solution when the valve spool sweeps across its underlap region. 2.31% of the valve spool travelling in the underlap region in the sawtooth test and the simulation results of the sawtooth command are depicted in Figure 17.

The analytical solution in the sawtooth test produces the same results as the numerical solution, but the analytical solution costs 1.23 seconds to finish and the numerical method costs 8.26 seconds to finish. Nonlinear behaviours can be observed in the results in Figure 17, which include the gradient change of the velocity curve after the system switches its state and some observed oscillations when the valve in the underlap region around null position.

Now assuming the underlap region of the Moog four-way is 20%, a sine wave motion command is sent to the four-way valve so that it only operates in the underlap region. Its pressure behaviours of numerical and analytical simulation results are depicted in Figure 18.

The modelling time costs of calculating chamber pressures when valve travel in underlap region is compared in Figure 18, the numerical solved simulation used 1647.17 seconds to finish, while the analytical solution only uses 8.50 seconds which is a significant reduction by a factor of 200.

The time consumptions of the three tests are depicted as in Table 3.

4. Conclusions

The component linking method is validated in Matlab Simulink, its simulation results show consistency with the published results. So that it can be utilised to achieve further research when the valve of the asymmetric cylinder drive system in the underlap region. How the chamber pressures behave when the valve in underlap is revealed by Simulink modelling, the inevitable pressure jump process when valve switches its state is numerically solved with factor method. To reduce the computing power required, an analytical solution for the factor is developed. Equation 47 proved the factor value is only related to the percentage of valve spool traveling in the underlap region and offered a solution to identify the root distribution percentage in underlap region. This analytical solution is approximately maximum 200 times computationally efficient than the numerical solution method.