Fatigue life estimation of hydraulically actuated mobile working machines using Internet of Things and Digital Twin concepts

The method and its implementation are proposed that allow estimating the fatigue life of mechanical structure of hydraulically actuated mobile working machines. The method does not use stress or strain gauges and relies just on the data about the pressure and position of hydraulic cylinders actuating the machine components. These data are measured in real-time during the machine operation and are transmitted via the internet for the processing in the cloud. The simulation model of the machine (Digital Twin) is used for the simulation of machine dynamics on the basis of received data. The Digital twin provides input for the finite-element analysis of machine’s components in order to calculate the stress history of each component. The results of these calculations are used for fatigue life estimation of each component. Applicability of the method is tested by comparing its results with the similar simulation verified by experimental measurements.


Introduction
The mechanical structures of mobile working machines suffer from material fatigue phenomena caused by cyclic loading. A common practice is to estimate the fatigue life on the design phase of the machine lifecycle [1]. The dissemination of fault-tolerant design approach increases the importance of on-line fatigue monitoring and prediction. Internet of Things (IoT) and Digital Twin concepts that have become popular over the last years offer new opportunities for solving these tasks. With the help of technologies and methods used in these concepts it is possible to gather real-time sensor data from the machine and simulate different processes occurring in the machine. The simulation can be performed concurrently with the machine operation. This work introduces the way of using real-time sensor data about the pressure and position of hydraulic cylinders to predict the fatigue life of hydraulically actuated mobile working machines at the operational phase of their lifecycle.
Different methods for fatigue life estimation can be utilized at different stages of the machine lifecycle. As the fatigue life of heavy equipment is evaluated at the design stage, its estimation during the early period of machine life is of low importance. The stress-based methods [2] can be used at this stage for cumulative damage calculation to monitor and predict the failure caused by improper use of the machine. Such misuse can generate the stress levels that differ significantly from those considered at the design stage. As a result, the real fatigue life changes as well. For example, hydraulically actuated mining drills are reported to be misused by scaling, which means that instead of drilling, the boom is used to pry loose rocks from the wall [3]. In-service estimation of fatigue life can detect such kind of unexpected work cycles of the machine in order to warn the owner of the machine and to account for them in the design of new machines.
When the cumulative damage reaches the value of crack initiation, the transition can be made to using the fracture mechanics methods in order to calculate the speed of crack propagation. This is especially important when machines remain in operation after the expiration of the design fatigue life which is often the case for the large and expensive machines like mining excavators or gantry cranes [4].
There are several ways for application of fracture mechanics methods to the old machines, where the cracks already exist. One option is to enter the position and other parameters of the cracks discovered by the inspection of the machine performed during the maintenance activities. The speed of propagation of discovered cracks can be monitored automatically in order to define the time of the next maintenance. This approach should not replace the inspection intervals defined by standards but could compliment them in case where the large or uncommon loads lead to the need for more frequent inspections.
It is also possible to detect the crack initiation zones automatically by analyzing the points with maximum stress levels along the lifecycle. Combination of this approach with subsequent application of fracture mechanics methods can be useful for remanufacturing process as part of sustainable practices. Fatigue damage obtained by the component of a machine during its lifecycle can serve as a basis for the decision about the possibility of using this component in remanufacturing process [4]. Automatic detection of crack initiation zones and estimation of crack propagation speed can also define the locations on the component for repair.
For old machines that have been used for years and that are most vulnerable for fatigue damage (such as mining excavators) the following approach can be proposed. The hydraulic system of a machine is modified by adding pressure and position meters to the cylinders. The load sensor, if it is absent on the machine, is also added or simulated. The stress history of the machine is gathered during some period of time that is large enough to cover all different conditions of machine operation. The length of the time period depends on the intensity of machine use and the number of different conditions in which the machine can be used. The gathered stress history is extrapolated back in time for the full period of machine exploitation prior to measurements. The resulting "augmented" stress history is used for fatigue life estimation of the machine. The detected hot spots with highest stress levels can define the zones for thorough inspection and the calculation of crack size and propagation speed can define the methods of repair. As a result the remaining fatigue life of a machine can be extended. The practicability of the proposed method depends on the repetitiveness of operations that machine performs. Mining excavators are good examples of the machines that perform repetitive operations for years.
The similar approach is used in useful life estimation for prolonging the life of used cranes [5]. The difference is that in case of dockside cranes and gantry cranes all work cycles are nearly the same and vary in weight and size of the loads. In such conditions the load history is represented by distribution of weight of the loads and total number of lifts during the life of the crane. The work cycles of mobile working machines, on the other hand, differ from each other and gathering real-time data with the help of IoT technologies improves the accuracy of fatigue estimation.  Computing and storing the stress history for each mesh element of each machine component requires a lot of computational and storage resources. It seems more reasonable to consider just some subset of the mesh elements. For this purpose a time period with the typical load history of the machine is analyzed and the elements of the mesh with highest stress levels are chosen as a reference set for fatigue life estimation. As higher stress levels provide shorter fatigue life, the calculation of fatigue life for this set of elements only is sufficient for fatigue life estimation of the whole component. The nodes belonging to the elements from the reference set are called the hot-spots. These are usually the points where the holes, grooves, fillets and other geometrical features are located which cause high local stresses under the load, called stress concentrations.
Finite-element analysis is usually applied for static problems when the simulated component does not move. In order to make it suitable for investigating the moving components of mobile machines a technique called "Inertia Relief" is commonly used. This technique is provided by the FEA software. It calculates the acceleration obtained by the simulated body as a result of the boundary conditions application. In order to simulate the equilibrium and eliminate the rigid body motion, additional boundary condition is generated that applies the opposite acceleration to the center of mass of the body or distributes the opposite acceleration among several elements of the model.
Utilizing the "Inertia Relief" technique, the whole process of fatigue life estimation for mobile working machines can be constructed the following way: 1.
Real-time data from a working machine is used as an input for the digital-twin model to simulate dynamics of machine components. This simulation can assume the machine component as a rigid or flexible body. In this study the rigid body motion is assumed and the flexibility is taken into account at the FEA stage only. 2.
During the simulation the acceleration of each machine component as well as constraint forces are calculated at every time step and stored for later processing. 3.
Using the mesh constructed from CAD drawing of each component and setting constraint forces as boundary conditions, a finite-element analysis is performed to calculate the stress in hot spots of each component. Acceleration of the component is utilized in Inertia Relief technique to eliminate the rigid body motion and provide quasi-static analysis. The time step for FEA, i.e. the moments in time at which the stress is calculated, can differ from that one used for simulation. The main requirement is that it should be suitable for accounting all vibrations of the component. 4.
The stress history obtained at step 3 is used for fatigue life estimation. This work uses the rainflow counting method to extract stress cycles from continuous stress history but other counting methods can also be applied depending on the fatigue estimation software being used. Applying some cumulative damage theory to the stress cycles extracted from the stress history the cumulative damage can be obtained and fatigue life can be estimated. This work uses the Palmgren-Miner rule for cumulative damage calculation but other theories can also be applied [6].
The hot spots considered at step 3 can be defined manually by the manufacturer of the machine or by its user relying on the experience of machine utilization. They can also be defined automatically. For this purpose the CAD drawing of each machine component can be covered with the nodes positioned from each other with some user-defined distance. These nodes can also be placed automatically during the mesh construction process. At the beginning of machine lifecycle the stress history is calculated and stored for all the nodes. After some period of time, that is large enough to cover a common cycle of machine operation, for example a week or a month during which the machine was used in all modes of its operation, the nodes with highest levels of cumulative damage are identified and defined as hot spots.
Following this approach it is possible to account for different manners of machine operation used by different users. Gathering statistics along the lifecycle of several different machines it is also possible to define the safest manners of operation for every type of machine and every kind of work performed with them (mode of operation). It can be accomplished by first choosing the machines with lowest levels of cumulative damage and then finding the differences between the manners of operation of these machines and other machines performing the similar tasks. Comparing the stress histories of different models of the machines executing the similar tasks it is possible to choose the best suitable model for the particular mode of operation. It can be useful for companies owing large fleets of different machines and operating them the similar way throughout the lifecycle, for example, for the mining companies.

Application of the proposed method to the test model
This section considers the test application of the proposed method of fatigue life estimation to the model of the mobile crane described in [7]. The crane being simulated is the hydraulically actuated mobile log crane PATU655 presented in the figure 1(a). This work uses the CAD drawings of the crane developed in [8].
(a) (b) Figure 1. The simulated mobile log crane and its kinematic model.
The model uses iterative Newton-Euler dynamic formulation (INEF) to simulate dynamics of the crane. The simulated mobile crane consists of the pillar that rotates around the vertical axis, the lift boom and the jib boom that are actuated by the hydraulic cylinders and the extension boom that is actuated by the hydraulic cylinder located inside the jib boom. The INEF formulation provides the system of differential equations (1) that is solved to get the values of independent coordinates of the system: the angle of the lift boom θ 1 , the angle of the jib boom θ 2 , the angle of rotation of the pillar θ 0 and the length of the extension boom L ( figure 1(b)). The actuator forces are represented in equation (1) This system of differential equations is derived analytically using the types of joints between the booms, the geometrical parameters of the booms and their inertial properties as the input values. Inertia properties such as the mass, inertia tensor and position of the center of mass can be obtained using the FEA software and CAD drawings of the booms together with their material properties. The system of differential equations is solved numerically, for example, by the Runge-Kutta fourth-order method to get the values of independent coordinates at each time step.
The cylinder forces are calculated using the equation (2): where P i are the values of pressure and S i are the values of cross-section area of the cylinder chambers, F fr is the friction force. The torques provided by the cylinders in the hinge joints are calculated using the geometry of the crane. Equations (1) and (2) show that the dynamics of the crane is defined by the values of pressure in the cylinder chambers and also by the values of displacement and velocity of the cylinder rods. These data can be measured by the sensors located on the crane and can be transmitted to the computing facilities in the IoT environment.
To get the values of constraint forces in the joints at each time step the following system of equations is considered: where M is the mass matrix of the system, A is the vector of linear and angular accelerations of the bodies comprising the system, f(F ext ) is the vector function that defines for each body and each component of its linear and angular acceleration the sum of constraint and external forces and moments acting on that body. This system of equations can be solved for constraint forces using the accelerations and external forces such as actuator forces and gravitational forces as input values. The system of equations (3) is just another dynamic formulation of the multibody system. Using it this way eliminates the need of solving a differential-algebraic system of equations. If the differential equation solver is implemented as a custom program, the calculation of constraint forces can be performed at each time step after the values of independent coordinates at this step have been found. This The Code-Aster software contains also some capabilities for fatigue analysis. For example, the CALC_FATIGUE command allows (among its other capabilities) the calculation of cumulative damage using the stress history and the S-N curve (Wöhler curve) of the material. This command uses the rainflow counting method to extract the elementary load cycles from the stress history. The cumulative damage is calculated as the sum of the damage associated with the elementary cycles using the Palmgren-Miner rule.
As the Palmgren-Miner linear damage hypothesis assumes that the total damage is a sum of contributions from elementary load cycles over the time, it is possible to develop an iterative procedure for fatigue life estimation during the lifetime of the machine. For this purpose the real-time data gathered from the machine are being split into the parts suitable for processing on the available computing resources. Multibody dynamics simulation can be performed in real time without considerable requirements for resources [7]. The main limitations are the random-access memory available for FEA software and the amount of time needed for stress calculation at all time steps of the period being processed. Depending on the available resources and the complexity of the mesh the size of time period for stress history calculation should be defined. This definition could be made manually or automatically on the basis of measurements of time needed to process the first datasets acquired from the machine. The best solution should periodically adjust the size of the time period being processed depending on the load of the computing system.
The selected size of the time period for processing defines the time boundaries for multibody simulation and the size of CSV files with the calculated constraint forces and actuator forces. After the finite-element analysis of that time period the calculated cumulative damage is added to the total damage of the machine component that is stored in the database. The maximum value of damage among all hot spots should be monitored in order to produce an alarm when it exceeds the preset limit.

Computational experiments
To validate the proposed method the simulation results were compared with the results of the prior research performed by Professor Aki Mikkola in [10] that was used as a reference. In his research A.Mikkola considered also the mobile crane PATU 655. He measured the stress in two points of the crane boom during two work cycles of crane operation defined by the input voltage profile of control valves and positions of the booms. The measured results were compared with the simulated ones obtained with the help of ANSYS and ADAMS simulation software systems. Figure 2 shows the input voltage of the control valves considered in the research. The input voltage was the same for all hydraulic cylinders of the crane. The duration of work cycles was 8 seconds. As the hydraulic circuit considered in the reference work differs from that one used in simulation model described in [7], in order to compare the stress histories of the crane its motion was replicated. The initial position of each cylinder that was 307 mm for the lift cylinder and 71 mm for the jib cylinder and the load mass 370 kg were also replicated. The extension boom was considered to be fully extended.
In order to simulate the similar motion of the crane as in the reference research work, the input voltage profile of the control valves was adjusted to provide the similar positions of the lift cylinder and the jib cylinder. The values of input voltage applied to the control valves were changed but the moments in time when the change was taking place remained the same as in the reference work. The slewing actuator in this study was simulated by applying the torque to the pillar of the crane. The torque was adjusted in time to provide the similar variation of angle as in the reference research work. Figure 3 shows the simulated positions of the actuators in comparison with the reference work.  The red lines in figure 3 show the positions calculated in this work, the thick black lines represent the positions measured in the reference work, and the thin black lines represent the positions calculated in the reference work.
The differences between the model used in this work and the model used in the reference research work are summarized in the following The pillar is rotated by the slew mechanism consisting of two hydraulic cylinders The load is considered to be a point rigidly attached at the tip of the crane The load is considered to be a body connected to the crane via a rigid link and a universal joint The booms are modelled as volumes using the 3D mesh constructed from the CAD drawings of the booms The boom's structure is modelled using linear 6 DOF beam elements, the lift arm is modelled using 14 elements and the jib arm -using four elements The test environment used in the experiment is shown in figure 4. The real-time sensor data were produced by the simulation program that was calculating the pressure in the piston side (Pa) and the blind side (Pb) of the lift cylinder and the jib cylinder as well as the positions of those cylinders (s1 and s2). These data together with the timestamp and the value of the torque applied to the pillar were sent by UDP protocol to the server side. The server program was receiving the data and storing them to the database. The received values were used as the input for the digital twin simulation program that was reproducing the crane motion using the dynamic model described in the section 3. This program was calculating the constraint forces and acceleration of each boom that were used as boundary conditions for the FEA software. In the experiment the open-source FEA software "Code_Aster" was utilized. The values of the stress in hot spots at every time step were calculated as a result of the finiteelement analysis and were stored as the stress history to be used in the estimation of fatigue life. is taken from [10]). In these points the stress gauges were located in the reference work. The points were located on the upper flange of the profile, their distance from the edge of the profile was 1/4 of the total width of the profile. The strain gauges provide the scalar values of the composition stress, the different components of which cannot be separated. In order to compare these values with the simulation results that consist of individual stress components, a superposition of the simulated stress components onto each other was performed in the reference research work to produce the composition stress. As the reference work does not provide any equation for calculation of the composition stress, the equivalent stress was used as a common way for the calculation of the scalar value from the stress tensor. The equivalent stress was calculated the following way: The stress histories consisting of the composition stress measured and calculated in the reference work and the equivalent stress calculated in this study are presented in figure 6.

Conclusion
The differences between the hydraulic models mentioned in Table 1 lead to the higher frequency of stress oscillations in the model created in this work. The large amplitude of the oscillations during the last two seconds of the simulation is caused by the low damping of the model created in the current