Digital twin based condition monitoring of a knuckle boom crane: An experimental study

This paper presents a novel approach for implementation of a digital twin for condition monitoring of a small-scale knuckle boom crane. The digital twin of the crane is simulated real-time in a nonlinear finite element (FE) program, where the estimated payload weight is used as an input. We implement an inverse method for estimation of the weight as well as its force vector direction based on physical strain gauge measurements. Additional strain gauges were utilized for validation of accuracy of the digital twin and inverse method. Based on a few physical sensor outputs, the digital twin allows for real-time determination of stresses, strains and loads at an unlimited number of hot spots. Therefore a digital twin can be an effective tool for predictive maintenance and product life-cycle management. In addition, condition monitoring of cranes during heavy-lift operations increases safety and reliability. The presented approach is described in a general manner and is applicable for various robotic manipulators used in the industry.


Introduction
Cranes are indispensable in many industrial operations. Offshore cranes are both used for installation operations and are an important part of supply chains of offshore installations. Cranes operate in demanding environmental conditions and are subjected to dynamic wave loads. In addition, cranes are often used for heavy-lift operations. When a crane operates, its configuration changes, which makes it even more difficult to determine the effect of the load combination on various components of the crane. Application of digital twins for condition monitoring, predictive maintenance and life-cycle management can be a reasonable solution, as, based on data from only few sensors, it provides real-time stresses in the entire structure and loads in the components such as bearings, shafts, cables, hydraulic cylinders, etc. Such data can be recorded and stored throughout the life-time of the machine. This is a cost-efficient solution since no additional physical instrumentation is needed.
In general, digital twin is a term used to describe a model or several models of a physical structure, mechanism, process or system, where each model is specifically derived for collection and generation of relevant data [1][2][3]. The model should also accurately describe properties and behavior of a physical system, such that the data obtained from the model is reliable. The concept of a digital twin for a mechanical system was demonstrated by a simple example of a fixed beam in bending in [1], where the CAD (Computer-Aided Design) model of the beam and stress analysis tools were used in combination with sensors attached to the beam. The authors suggested that the next step of the research was to implement the concept for a moving mechanism. Some research on implementation of digital twins for mechanisms is already available. An interesting case of a digital twin for a helicopter rotor system was presented in [4], where the loads and displacements at the joints of the mechanism were calculated based on sensor data and dynamic multibody simulation. The authors suggested that knowledge of exact loads applied on mechanical components in real operation would improve prediction of maintenance. Some other examples of application of digital twins for condition monitoring and predictive maintenance within aerospace industry were given in [5][6][7]. In addition to digital twins of mechanisms in operation, similar concepts and experimental studies can be found for manufacturing machines and tools. For example, digital twins for CNC machines were discussed in [8][9][10][11]. Importance of digital twins for autonomous systems in manufacturing was discussed in [12]. A review of scientific literature reveals that many methods for implementation of digital twins for mechanisms, like cranes, can be applicable for autonomous systems used in manufacturing. This is especially relevant for robotic manipulator arms used in manufacturing, since a crane is mechanically very similar to a robotic manipulator. Digital twins for robotic manipulators were discussed in [13][14][15][16].
The main contribution of this work is a novel approach for implementation of a digital twin monitoring solution for a knuckle boom crane, which is also applicable for other robotic manipulators used in the industrial processes. This work is an extension of the short conference contribution [17]. The proposed digital twin solution was implemented for the case of a small-scale knuckle boom crane and verified by experiments. The proposed approach can be efficiently used for predictive maintenance and life-cycle management of cranes and their components, by continuously collecting data about motion of a crane and loads a crane and its components are exposed to. In addition, the possibility of real-time data acquisition allows for increasing reliability and reducing risk under heavy-lift operations. Industrial cranes as well as robotic manipulators are often an important part of the manufacturing process, and the proposed solution can be integrated into digital twin of an entire production line, which would contribute to efficiency and optimisation.

Preliminaries on inverse methods
A problem in structural dynamics where the response of the structure is found from applied known loads is referred to as a forward problem. An inverse problem for digital twins refers to the process of estimating the excitation the structure is subjected to from physical sensor measurements of the structural response. The method used to do this is referred to as an inverse method.

The inverse problem
The inverse problem is ill-posed [18]. In a well-posed problem, there is a unique solution to any reasonable measured data. For the inverse problem, this is not necessarily the case. Some problems are: • The solution does not have to be unique. • The measured data is most likely incomplete (all deformations are not measured). • Noise and drifting in the sensor data can lead to undefined solutions. • The location of the external forces may be unknown.
These problems must be addressed when developing inverse methods. For the case of digital twins of mechanisms like a crane, there are additional requirements for the method: • It should be able to run continuously with real-time data. • It should be computationally efficient for real-time calculations of large systems. • It should be stable over long time periods. • It should be able to handle large displacements of mechanisms.

Method basis
We propose a novel inverse method to meet the demands of digital twins. It is based on a known equal number of applied forces and strain gauges at known locations and utilizes finite element simulation. Rølvåg et al. [32] described how a virtual strain rosette can be formulated using a zero-stiffness membrane element connected to nodes around the area where strain is to be determined. For this membrane element, there exists an element strain-displacement matrix which gives a vector of strains as where B e is the strain-displacement matrix of the element, and v e is the strain rosette element local displacements. These strains must be transformed into Cartesian coordinates with the transformation matrix T r , and the strains in one direction (leg) of the virtual rosette isolated with t l . Expanding the equation further to relate strains to global displacements v, gives  (2) with N e being a matrix of ones and zeroes for selecting the strain rosette element displacements. T e transforms these displacements into the element local coordinate system. Finally, L expands the system to include linearly constrained displacements. This means that strains are linearly related to the displacements of the systems, where it is assumed that the matrices in Eq. (2) are constant.
In a linear static system, the displacements are linearly related to the external load Q by the stiffness matrix K = Kv Q. ( Assuming that only applied forces F contribute to Q, Eq. (3) can be expressed using F and a distribution matrix q = Kv qF. (4) For a system with n applied forces and n strain gauges, the relationship between them can be expressed by Solving this equation for the applied forces gives: The matrix S is square due to the equal number of strain gauges and applied forces, and invertible if strain gauge locations are chosen correctly. Eq. (6) gives the basis for the inverse method. Depending on the model reduction used, K might not be invertible. The matrix S can then be found by introducing unit forces, and recording the strains The recorded strains are the columns of the inverse compliance matrix S 1 n n n

Mechanisms
Problems with the relations above arise when considering multibody mechanisms, because large displacements will be induced by actuators. If strain gauges and applied forces are located on the same component, and the forces are defined in the component local coordinate system, only gravitational effects need to be accounted for. On the other hand, if they are spread across different components, the compliance matrix becomes a function of actuator values a (lengths and angles) Large displacements will also cause changes in strain readings due to gravitational forces. These need to be accounted for when using the measured strains in Eq.
where A a ( ) is a vector function of actuator values representing strains caused by changes in gravitational forces from the position where the strain gauges were zeroed. Internal dynamic effects may also contribute to strain readings, as well as high frequency measurement noise. These effects can be removed by means of a low-pass filter, eliminating frequencies above the highest frequency of the applied load.

Cable and pulley system
The effect of a cable and pulley system can be added to the method by considering static force equilibrium at the end pulley (see Fig. 1). Three mutually perpendicular forces representing the reaction forces of the pulley are used as F , and the equilibrium equations give the cable tension. This cable tension F c is applied to the digital twin as forces on each pulley and the attachment point of the cable. The effect of these forces on strain gauge readings must also be accounted for when calculating the forces. This is done in a similar manner as for gravitational forces in Eq. (10), using a vector function W a ( ).

Drifting sensitivity
Measurements tend to be prone to drifting (low frequency) and signal noise (high frequency). Both causing deviation from the desired measurement value. While signal noise can be effectively removed by a low-pass filter, the drifting can have frequencies coinciding with load frequencies, and thus can be difficult to filter out. When forces are calculated from a drifted sensor signal with the proposed inverse method, the drift will be read as actual deformation; the method has no way of distinguishing it from the actual deformation. Thus, we get artificial forces that the physical twin was not subjected to. The effect could be reduced by compensating for the drifting, e.g. compensating for temperature changes. However, if the setup used is very sensitive to drifting, meaning a very small drift causes large artificial forces, this might not be sufficient. An effort should therefore be made to minimize the effect of drifting. For comparing different setups, we can define a sensitiveness index s as where drift is some typical vector of drifted signal values, and the physical interpretation of s being the total resulting artificial force. The sensitiveness thus depends on S. By optimizing the sensor distribution, s can be minimized.

Implementation
The digital twin was implemented using Fedem, a nonlinear finite element software made for fast digital twin computations, utilizing efficient CMS model reduction techniques [33]. Fedem also embeds efficient virtual strain gauge calculations enabling automated setup of the inverse method for any model. The inverse method itself was implemented in a Python script, which received sensor data, estimated forces and controlled the real-time simulation.
The physical and digital twin setup is given in Fig. 1. The payload was hanging in a cable running over pulleys to the base of the crane. Strain gauges , UT UR and LT were used for Eq. (10) to estimate the forces. The forces on the end pulley were implemented as three mutually perpendicular forces in the upper boom local coordinate system. Two independent gauges UB and LL were used for validation.

Finite element model
The finite element models for each of the crane components were created with NX Nastran. 10-noded CTETRA elements were used for all parts [34]. The Fedem software imported the Nastran files, and connected the components by node to node connections. Rigid elements of type RBE2 were used at joints to facilitate this.

Actuators
A kinematic diagram of the crane showing the actuators is shown in Fig. 2. As can be seen, the motion of the crane and payload is spatial, and not restricted to a linear or planar motion. Encoders on the actuator motors were used to record the motion of the crane. Their digital twin counterparts were modeled using an approach very similar to the penalty spring constraint method. The slewing actuator was modeled as a revolute joint with high rotational stiffness, and controlled by introducing a stress-free angle change. The joint was connected to the base and ground, which meant that the stress-free angle change forced the whole crane model to rotate.
The lower and upper actuators were modeled as linear springs and dampers. A stress-free length change was used to control the length of each actuator. The stiffness was chosen to be 1.0 × 10 8 N m for both actuators. This is approximately equivalent to the stiffness of the ball screw inside the physical actuators. With no damping, the stiff spring elements exhibit unwanted vibrations. The damping coefficient was found iteratively by comparing measurements and simulation results to be 2.0 × 10 7 N s m . The mass of the actuators was added as concentrated masses at the spring end points. The mass of the side drive, motor and fixed cylinder was added to the lower end node, and the mass of the moving piston was added to the upper end node.

Strain rosettes
Virtual strain rosettes are implemented in Fedem with the formulation described in Section 2.2. A visualization of such a rosette is shown in Fig. 3b. The zero-stiffness membrane element is connected to four nodes in the same location as the physical rosette shown in Fig. 3a. If the type of the virtual rosette is set correctly, and the directions matched with the physical gauge, they become true twins. Filtering of the measured strains was achieved through a digital Butterworth low-pass filter. A cutoff frequency of 6 Hz and filter order 6 gave the best results. This filtering was applied only when using the measurements for force estimation, and thus the measurements shown as results are unfiltered.
The distribution of strain gauges proved to be important. The sensitivity index s was used to evaluate different distributions. A typical drifting observed in the physical measurements was 1 µm m across all gauges. This was most likely due to small temperature changes, as no efforts were made to counter this effect. If strain gauges , UT UR and UB were used, it resulted in a sensitivity index s = 121.7 N. Such an artificial force is very high, considering that a typical payload was 220 N. Replacing UB with LT however, reduced the sensitiveness to s = 6.8 N. The sensitiveness of the compliance matrix could potentially be improved even further by altering sensor placements, but for the purpose of this paper, the combination of , UT UR and LT was satisfactory.

Lookup tables
Some matrices used in the proposed inverse method are dependent on actuator values when considering mechanisms. By choosing an appropriate number of points in the range of each actuator where the inverse method matrices S and A are to be calculated, a hyperdimensional array of matrices was constructed (the number of dimensions depending on the number of actuators). In simple terms, each entry to a matrix or vector was saved as a lookup table. This is illustrated in Fig. 4 for an × n n compliance matrix with two actuators. Actuator 1 has k interpolation points, while actuator 2 has m points. The value of any S ij can be approximated by linear interpolation in the table to the right in the figure. Any number of actuators can be added by increasing the number of dimensions of the lookup table.
The calculation itself was automated by iterating through all possible combinations of interpolations points of the actuators and calculating the required matrices in each position. The entries of the gravitational force compensation matrix A were found from strain readings at zero load for the given position. The compliance matrix S was found by the method give by Eqs. (7) and (8) using unit forces. Finally, the cable force compensation matrix W was found in a similar manner through applying unit cable forces, though no matrix inversion was necessary. The required number of interpolation points depends on the nonlinearity of the matrices. For the crane, five interpolation points for each of the two linear actuators proved to be sufficient.

Results and discussion
Several tests were performed in order to benchmark the digital twin solution, and two representative runs with different payloads are presented in this paper. Each run started by lifting the payload with the crane, manually inducing abrupt movement in all three actuators to initiate spatial payload sway, and then finished by lowering the payload back to the ground. The actuator movement for each run in shown in Fig. 5.
The supplementary video shows a real-time simulation of the crane, with some result plots. The accuracy of the inverse method and model was evaluated in two ways. First, the independent strain gauge measurements were used for the verification and evaluation of accuracy of the virtual strain gauges. Second, the estimated cable tension was compared to an expected static tension (no payload movement). Figs. 6a and 7 show these results from the run with the heaviest payload ( = m 32.2 kg). In the strain measurement comparison there is good correlation between strain measurements, but the digital twin had some trouble replicating high frequency oscillations. This was expected as the inverse method assumes quazi-static behaviour. Some offset in mean values was also observed, and can most likely be attributed to measurement drifting due to e.g.

Fig. 3.
A strain rosette UR on the physical and digital twins. The physical rosette is covered by a protective substance, but is illustrated with overlaid graphics. Fig. 4. A compliance matrix with two actuators as a hyperdimensional array. Each entry in the matrix (left) is a lookup table (right).

T. Moi, et al.
Engineering Failure Analysis 112 (2020) 104517 temperature changes. Looking at the force plot in Fig. 6a, large oscillations can be observed. This is due to inertia forces from the sway of the payload, and thus only the mean value of these oscillations should be compared to the expected static F c . The effect of sensor drifting in strain measurements is apparent here also, with some low frequency variations in mean F c . Though inertia forces are included for the payload in the inverse method, the inertia forces of the crane itself is not. The Fedem software includes inertia forces of the crane components in the simulation, but the strain caused by this measured on the physical twin will be falsely converted to equivalent payload forces by the inverse method. The authors did not consider this to be a significant source of error in this particular case, but if the inertia forces in a different case are significant, they should be compensated for.
In order to evaluate if the results are sufficiently accurate, their use case must be considered. If the strain gauge measurements are used for fatigue damage calculations, a peak valley extraction is usually performed to reduce the number of data points. In this process, oscillations with an amplitude below a specified gate value are ignored. These small oscillations are considered to have no effect on the fatigue of the component. Such a processing of data for the independent gauge LL is shown in Fig. 7c. It is clear that the inaccurate representation of high frequency oscillations does not affect the accuracy of fatigue damage calculations when data is processed in this manner.
In fact, the strain results from the digital twin can produce more accurate input for fatigue damage calculations. In the results of this paper, both the strain gauge measurements and simulation results are zeroed at the start of each run, but the crane components will actually be pre-strained due to the cranes own weight, giving a higher mean strain. While this pre-strain can be hard to identify using physical measurements on an assembled crane, they can be found easily with the digital twin by zeroing the virtual strain gauges before applying gravitational load to the simulation.
Some additional results of interest for condition monitoring are shown in Fig. 6b-d. These are all quantities that could be hard to measure, at least without modifications to the crane. Joint forces and moments are shown as totals for each joint, but are available as individual components in global or local coordinate systems as well. Fig. 8 shows the main results from the run with the lighter payload ( = m 13.5 kg). Here we observe a more pronounced effect of  sensor drifting in both the force and strain plots. This can be explained by the fact that the magnitude of the strain signal is lower when using a lighter payload, making the drifting relatively larger compared to the actual signal. Thus, if higher accuracy at this load level is required, efforts should be made to compensate for the sensor drifting.

Error estimates
Quantifying the error between physical and virtual strain gauges is not straightforward. A simple time step by time step comparison would yield an unreasonable error due to the small offsets in time between the twins. Since mainly peak values are of interest for structural monitoring, these were chosen as basis for the error calculation. These could be paired between physical and virtual gauges, and thus removing the time offset problem. However, the number of peaks between the physical and virtual gauges were not the same. This was due to edge cases close to the gate value in the peak valley extraction, where the peak in one twin is barely  included while the other is not. Therefore, peaks where a pair between physical and virtual measurements could not be found within 0.5 s time difference were ignored. Depending on the use case of the digital twin, different measures of error would be of interest. If crane overload is to be monitored, then the magnitude of the strain is of interest. Since the value of a peak could potentially be zero, giving a division error, a mean absolute percentage error in relation to the mean absolute value was used for strain magnitude error: Where m is average magnitude percentage error, n is the number of measurements and p t , and s t , are the physical and virtual strain measurements for a given pair. p denotes the mean absolute value of the physical measured strain peaks.
Overloading can also be monitored by the cable tension. The error was calculated by comparing the mean estimated cable tension F c between 5 s and 50 s, where the payload was fully lifted, to the expected static tension mg: If the measurements are used for fatigue monitoring, the strain amplitude is the most important. The error in amplitude, a , was evaluated by comparing strain difference between peaks: The error results are displayed in Table 1. In general, the run with the heavier payload measured the least error. Strain gauge UB was the most accurate, which is expected since it is located close to the gauges used in the inverse method. LL , located further away, showed larger errors. In relation to the proposed sources of error, the equivalent temperature difference of the strain magnitude errors were between 0.14°C and 1.36°C (assuming a thermal expansion coefficient of 10 × 10 −6°1 C ). It should also be noted that the error calculated with this procedure may be highly dependent on the chosen gate value, depending on the data.

Conclusion
In this work we presented an implementation of a digital twin with an inverse method using strain gauges as load sensors. A formulation for virtual strain gauges was used to find a linear relation between applied forces and measured strains. This relation provided the basis for an efficient nonlinear method suitable for digital twins with large displacements.
Experiments were conducted using a small-scale knuckle boom crane. A payload was lifted with a cable and moved to induce sway. Three strain gauges and encoders on the actuators provided the required input data for the digital twin. Two independent strain gauges were used to benchmark the accuracy of the digital twin solution.
The physical measurements showed good correlation with simulated strain time histories proving that both the crane model and inverse method are accurate. Cable tension was estimated within 1.4-6.7% of expected static tension, strain peak magnitude within 1.6-16.8% and strain amplitude within 6.0-13.5% of physical measurements, depending on payload weight and strain gauge location. This means that virtual sensor outputs from anywhere on the crane model can provide strain and load time histories for structural monitoring and fatigue damage estimation without further physical instrumentation.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgment
This work was partially funded by the Norwegian Research Council, SFI Offshore Mechatronics, project number 237896. The authors have also received support from SAP Norway Engineering Center of Excellence.