A parameter optimization method for stress simulation of double horse head pumping unit dynamics model with cable

3.1. The force model of bushing joint

As shown in Fig. 2, bushing joint presents the spring and damper forces acting between two parts. It defines the force and torque magnitudes using six Cartesian components ($F_x$, $F_y$, $F_z$, $T_x$, $T_y$, and $T_z$), which are a linear function of the translational and rotational displacement between two coordinate systems moving with the two parts [8, 9].

The following constitutive equations define how a force and a torque were applied to the parts according to the displacement and velocity between marker $I$ and $J$ in a bushing joint:

The bushing element applies an equilibrating force and torque to the $J$ marker in the following way:

According to the calculation method of materials mechanics, the coefficients to calculate bushing joint force are defined as follows:

Fig. 2The spring and damper forces acting between two parts in bushing joint

3.2. Elastic-damping contact force model

To calculate the contact forces for the steel cables and horse head in the multibody dynamic model efficiently, the numerical description of contact force model needs to pay special attention. For the modeling of the contact under ADAMS, we have used IMPACT-Function. In this method, ADAMS/Solver can give a continuous flow responses of accelerations, velocities, positions and forces for all contact elements.

The normal contact force has two components: rigidity and viscous damping. The component of rigidity is a function of the penetration $\delta$, while the viscous damping is a function of the speed of penetration.

In this model, the normal force of contact is given by the following relation [13]:

The instantaneous damping coefficient $C$ is defined as:

4.1. The 3D model of pumping unit

To simulate the dynamic stress of DHH pumping unit, a multibody dynamic model was built with software ADAMS. A CYJS10-5-37HB type DHH pumping unit is used as the research case, and the working and structural parameters of this pumping unit were listed in Table 1. The 3D geometric outline of pumping unit was first constructed by Solidworks, and then was imported into ADAMS. To maintain the kinematics and dynamics constraint, a series of joints were added to the model. The three-phase asynchronous motor, gearbox and support were connected to the base in real working condition; therefore, they were fixed to ground. The walking beam can revolute around steel support bearing, so they were assembled to support pin with revolute joint. Cranks was linked to the connecting rod and the shaft of gearbox with revolute joints respectively. The two horse head were connected to the walking beam with fixed joints.

The area near two pin holes in back horse head often occurs fatigue cracks, and this paper aims to simulate the dynamic stress distribution of back horse head. Therefore, the back horse head connected to the walking beam was created as a flexible body. This unit was firstly meshed in ANSYS software by using the Solid45 element, and then was created with two external points in the center of each pin hole, finally was imported into ADAMS software through a modal neutral file (MNF). The flexible back horse head is shown in Fig. 3, and its details of meshing characteristics are listed in Table 2. The back horse head was connected to the walking beam by two fixed joints between the external points and pins in walking beam.

Fig. 3The flexible model of back horse head

4.2. The model of steel wire ropes

In this pumping unit, the crossbeam drives the back horse head by steel cables, and the front horse head also drives the polish rod with steel cables. However, ADAMS software have no parts to model cable, usually the cable was modeled by cylinders linked with bushing joints. According to the structural parameters of given type pumping unit, a cylinder 10 mm diameter and 50 mm long was created and connected to the crossbeam with revolute joint, and then a series of cylinders in the same size were built by ADAMS macro along the route of real steel cables until it reached the fixed point in the back horse head. To simulate the flexibility of steel cables, bushing joints were also applied between each two cylinders by ADAMS macro. In the bushing joint, the stiffness $K$ and damping coefficient $C$ in different directions used the recommended values.

In the real working condition, the steel cables drives the back horse head and walking beam rotating around the steel support. The contact points between steel cables and back horse head are changing with the angle of walking beam. Therefore, a series of contact forces between back horse head and cylinders in steel cables should be added besides the fixed joint. This model employs the IMPACT function available in the ADAMS function library to simulate the contact force. The stiffness $K$ and damping coefficient $C$ in the IMPACT function also used the recommended values.

The steel cables between the front horse head and the suspension point was modeled by the same method. Fig. 4 demonstrated the rigid-flexible model of DHH pumping unit built in ADAMS.

4.3. The load and simulation

Suspension point load and gearbox shaft motion are the two loads of this model. A 5 rpm gearbox shaft rotate motion was applied to the end of gearbox shaft near cranks, and an experimental tested and smoothed suspension point load (as shown in Fig. 5) was applied to the suspension point in this model.

The selection of solver and its integrator are key factors that affect the solution efficiency and accuracy. This simulation used Gear Stiff integrator (GSTIFF) and (SI1) integrator formulation. The integration tolerance was set as 0.00001 and the step size was 0.00001 second.

Fig. 4The DHH pumping unit dynamic model

Fig. 5The experimental suspension point load of pumping unit

4.4. Experimental test

To verify the simulated results, an experimental test was conducted for the CYJS10-5-37HB type DHH pumping unit to obtain the suspension point indicator diagram and dynamic stress process in the same working condition. In this test, a wireless indicator diagram collector, as shown in Fig. 6 left part, was employed to collect the load and displacement signal of suspension point, and then the tested signals were transformed to the indicator diagram illustrated in Fig. 5.

Fig. 6The test of indicator diagram and stress for a CYJS10-5-37HB type DHH pumping unit

To test the stress changing process of area near pin A and pin B, four vertical arranged strain gauges with full bridge circuit, as point 1-4 demonstrated in Fig. 6 right part, were bonded to the surface. The strain signal was recorded with the indicator diagram synchronously, and then was transformed to stress accompanied with crank rotation angle in computer.

4.5. Results and discussion

We choose one test point from each pin, and the model experimental tested dynamic stress of point 1 and point 3 were shown as blue line in Fig. 7(a) and Fig. 7(b) respectively. We can observe that dynamic stress in two points both exhibit periodic fluctuations, and the alternating stress amplitude is 12.87 MPa for point 1, and 21.87 MPa for point 3. Though they are both in the range of material allowable stress, point 3 appears tensile and compressive stress alternately, and this made the components easier to occur fatigue cracks after a long and high frequent cycle working period.

Fig. 7The comparison between simulated and experimental tested stress

a) Point 1

b) Point 3

The model simulated dynamic stress of point 1 and point 3 were also shown as red line in Fig. 7(a) and Fig. 7(b) respectively. We can observe that significant differences exist between simulated and experimental stress in amplitude, and there are also many vibrational peaks in the simulated stress due to the unreasonable cable flexibility and contact force, while the tested stress are smooth relatively. After serval trials, we found the differences in amplitude are mainly caused by unsuitable flexibility of simulated steel cables, and the vibrational peaks is probably due to the unstable contact force. While the stiffness $K$ and damping coefficient $C$ are key factors for performance in flexibility and contact force of cables, therefore, by adjusting the parameters of the stiffness and damping coefficient is an alternative way to improve the model.

5.1. The chromosomes and fitness function of GA-based MAs

In this model, the stiffness and damping coefficient of bushing joint and contact model were treated as variables, and there are six stiffness coefficients $K_{bi}$ and six damping coefficients $C_{bi}$ in bushing joint of steel cables and one stiffness coefficients $K_c$ and one damping coefficients $C_c$ in contact model between the steel cables and horse head. While the $K_{22}=K_{33}$, $K_{55}=K_{66}$ in Eq. (4), here we use one variable $K_{b23}$ to stand for the equal $K_{22}$ and $K_{33}$ to reduce the number of variables, and $K_{b56}$ stands for the equal $K_{55}$ and $K_{66}$. The damping coefficients $C_{bi}$ are defined in the same rule. Therefore, the chromosome is defined as the permutation of these parameters:

The objective of this optimization is to find the most similar simulated stress of area near pin A and pin B compare to the experiential one, that is, minimize the error between simulated and experiential stress in the same test point. This paper utilizes the mean square error (MSE) between simulated and experiential stress in each test point to evaluate their similarity, and we give the MSE of each test point an equal weight, therefore, the objective function is defined as:

Minimize $f(K,C)=MSE(Q_{As},Q_{Ae})+MSE(Q_{Bs},Q_{Be})$.

Subject to $K_l{<K<K}_u$, $C_l{<C<C}_u$.

Subscript $l$ and $u$ are the upper and lower value of the variable range.

5.2. Algorithm and flowchart of optimization method

This parameters optimization operation was implemented through the co-simulation between ADAMS and MATLAB. The calculations of MAs were conducted in MATLAB while the dynamics simulation of pumping unit was executed in ADAMS.

The parameters optimization process was implemented according to the following steps:

Step 1: The population, that is the parameters permutation including the stiffness and damping coefficient of bushing joint and contact model, were initialized randomly;

Step 2: The macro file of DHH pumping unit dynamics model was created by using the parameters in each individual of population, each individual corresponds to one macro file;

Step 3: ADAMS run all the macro files automatically, each individual obtained a pair of simulated dynamic stress near pin A and pin B in back horse head;

Step 4: MATLAB imported the simulated and experimental stress, and calculated the MSE as fitness function for each individual;

Step 5: If one of the fitness was equal to the given threshold value, or the generation reach the max number set previously, the calculation was finished, and the best results was obtained.

Step 6: If the operation continues, a new generation of offspring were selected by using the Roulette Wheel rule from parents;

Step 7: The new offspring conduct a greedy crossover based hill climbing to search some elite offspring. The elite individual is taken as one parent and another parent is selected from the current population using a roulette wheel. Then, a special uniform crossover is executed between these two parent individuals to generate an offspring. The offspring will replace the elite individual in the current population if it has a better fitness than the elite [14];

Step 8: Subsequently, the offspring were further refined by a steepest mutation-based hill climbing method. The chromosome only changes several bits randomly when executing one mutation operation on it. The elite individual is picked out from the current population and several random bits are changed. If the newly mutated individual has a better fitness, it will replace the elite individual [19];

Step 9: When the new generation of population was formed, this algorithm returned to step 2 to find a better fitness.

Based on the steps discussed above, the flowchart of model parameters optimized method was shown in Fig. 8.

5.3. Simulation

In the optimization process of model parameters, the steel cables of front and back horse head use the same stiffness coefficient $K$ and the coefficient $C$. The recommended values $K_r$ and $C_r$ were defined as the center of value for variables $K$ and $C$, and the variables $K$ and $C$ of each individual were created in the value interval from 10^-4 to 10⁴ multiples of value center $K_r$ and $C_r$ respectively.

In the optimization process of MA, the MA parameters plays an important role on the optimization results. Furthermore, due to the plentiful variables in this model, the variables may easily fall to local optimum values. Therefore, for searching a better global optimum values, orthogonal test method was employed to find the best MA parameters for this model. The best MA parameters are as follows: the size of population is 80, the maximum iterations are 60, the crossover probability is 0.9, and the mutation probability is 0.05.

5.4. Results and discussion

The optimized parameters, that is, stiffness coefficient $K$ and coefficient $C$, MSE and relative errors were listed in Table 3. To verify this method, the simulation results using recommended parameter values were also shown in Table 3. We can see that obvious differences of parameters exist between two models, and the MSE and relative errors of the improved model is much lower compared with that of recommended parameters model.

The simulated stresses of the improved model for point 1 and point 3 were illustrated in Fig. 9, and we can observe that, compared with the results simulated by the recommend parameters values under the same working condition in Fig. 7, the amplitude and distribution of simulated stress in Fig. 9 were closer to experimental stress in amplitude, and the peaks were suppressed significantly in the simulated stress. While there are still some differences because of several reasons: (1) The bushing joint and contact force model still exist theoretical shortage need to improve compared with real working condition; (2) The cylinders connected with bushing joints are a series of discrete rigid bodies, their dynamic performance are close to, but different from real continuous flexible steel cables; (3) Some structural details were simplified in the model, so the mass centers of components are some different. Despite the differences in the simulation results, it is a novel way to improve the dynamics response simulation accuracy for mechanism with steel cables.

For the purpose of verifying the applicability of optimized parameters in another way, we use the optimized parameters to simulate dynamic stress with the same pumping unit model in a different working condition, which the crank rotating speed is 4 rpm and an experimental tested suspension point load was applied to the suspension point. The simulated dynamic stress by the optimization parameters model was shown in Fig. 10 together with the experimental tested one, and the detailed comparison was given in Table 4. We can observe the amplitude and distribution of simulated dynamic stress is similar to that of the experimental tested dynamic stress, and though bigger than the MSE and relative error in Table 3, they are still in an acceptable range. So, it is feasible to simulate the dynamic stress of model for different working conditions by using the parameters optimized from one working condition, which can achieve an economic efficiency by the reduction of experimental tests. The model built in this paper can not only be used to simulate the dynamic stress of components for a safety assessment but also be employed to balance the crank torque for an energy saving.

Fig. 8Flowchart of model parameters optimized method

Fig. 9The comparison between experimental stress and simulated stress with optimized parameters

a) Point 1

b) Point 3

Fig. 10The comparison between experimental stress and simulated stress with optimized parameters in another working condition

a) Point 1

b) Point 3