Co-Simulation Platform for Simulating Heavy Mobile Machinery With Hydraulic Actuators and Various Hybrid Electric Powertrains

Computer simulations are important tools for evaluating hybrid electric powertrain concepts. In heavy mobile machinery, the powertrain undergoes traction and often hydraulic loads. The dynamic loads on hydraulic actuators of heavy mobile machinery can be studied using detailed physics-based simulation models. The effect of coupling hydraulic actuators to various configurations of hybrid electric powertrains is studied using a co-simulation approach in this study. The case example of a tractor modeled using multibody dynamic approach combined with four powertrain configurations with different topology and hydraulic actuator coupling is investigated. A farm tractor was simulated in a virtual environment with deformable terrain in a soil digging and dumping work cycle. The coupling of hydraulics to the electric energy source in a series-parallel topology achieved fuel savings of 10.5% and total energy savings of 12.3%. The presented modeling method can be used to simulate various heavy mobile machinery to evaluate hybrid electric powertrains using work cycles created by human operators.


I. INTRODUCTION
Computer simulations of heavy mobile machinery has become a central part of analyzing the performance and energy efficiency of hybrid electric powertrain concepts. In practice, it is possible to perform simulations that can describe the actual dynamics of such machinery using multibody dynamics [1]. Multibody machinery models are often coupled with hydraulic actuators to capture the realistic response of the system [2]. In modeling hydraulic circuits, a constant pressure source is typically considered for simplicity [3], [4]. However, in reality the hydraulic pump is driven by an engine or motor, the energy consumption of which accounts for a significant portion of the demanded power of such machines. To accurately simulate heavy mobile machinery, detailed hydraulic circuits that are coupled with an energy source must be modeled, so that all the dynamic The associate editor coordinating the review of this manuscript and approving it for publication was Wonhee Kim . loads are accounted for. Therefore, it is essential to consider the effects of hydraulic actuators on the powertrain, especially for energy consumption investigation. Due to the high power requirements of such machines, hybrid electric powertrains must be simulated and investigated to decrease the fuel consumption of hydraulic heavy mobile machinery [5].
Hybrid electric architectures can improve the fuel economy of heavy mobile machinery [6], [7]. Novel driveline concepts for heavy mobile machinery were proposed to improve the fuel economy and to achieve compactness in [8] and [9], respectively. Many studies have investigated control strategies for hybrid excavators to improve fuel economy [10], [11], [12], [13], [14]. A comparison of control strategies made in [15] showed that suitable controllers can achieve over 10% improvement in fuel economy for hybrid loaders. The efficiency of a hybrid excavator was improved by 8.8% compared to a conventional excavator by optimizing the system parameters in [16]. Hybrid electric energy storage systems have also been studied to bring improvement in performance and efficiency of heavy mobile machinery [17], [18]. Hybrid electric tractor powertrain architectures were investigated in [19] to improve fuel efficiency and peak power performance. These studies have focused on hybrid electrification of heavy mobile machinery while accounting for the hydraulic energy consumption. However, instead of detailed models of the hydraulic circuits and the machinery, the hydraulic loads have been directly used from measured experimental data or typical duty cycles.
Some studies have developed detailed hydraulic circuits for heavy mobile machinery when proposing novel powertrain component designs [20], [21]. A hydraulic drive consisting of a servomotor directly driving a pump was proposed in [22]. An integrated electro-hydraulic energy converter was developed in [23] for long booms in heavy working vehicles. Novel electro-hydraulic hybrid excavator powertrain designs have also been proposed to achieve improved control and energy regeneration [24], [25], [26]. For example, the energy-saving system developed in [27] used a hydraulic motor and an electric generator to achieve pressure compensation and energy regeneration in hybrid excavators. The powertrain designs developed in the above studies use data measured from test benches or prototypes to validate the designs or run control algorithms. In addition to experimental results, simulations have also been used to improve the control [28], energy regeneration [29] and hybridization factor [30] of hybrid electric hydraulic machinery. However, the simulation models of the machinery used in studies [28], [29], [30] are simple and do not account for the detailed model of traction and hydraulic loads.
Complex multibody simulation models of heavy mobile machinery comprising variable displacement pumps coupled with the engine have also been used to account for the detailed dynamics of the entire machine [31], [32]. Mechanical models accounting for the forces acting on the hydraulic systems have been coupled with energy recovery systems proposed in [33] and [34] using load sensing pumps, valves and accumulators. These models captured the dynamic behaviour of the machinery. However, it failed to include the two-way coupling interaction of the hydraulic pumps and the electrified powertrains. The coupling of a multibody micro excavator model with an electrically actuated hydraulic system is presented in [35]. But the electric drives in this study are not connected to an energy source and the work cycle considers a stationary excavator, thereby ignoring the powertrain dynamics as well as the influence of the working environment on the dynamic performance of the excavator. To evaluate the real-time control strategy proposed in [36] a wheel loader model was developed using multibody dynamics approach with a load sensing hydraulic system and pressure compensating control valve. The dynamic loads were taken from measurements and duty cycles.
As evident from the studied literature, powertrain architecture, hydraulic circuits and mechanical models of heavy mobile machinery are all important aspects in the hybridization of heavy mobile machinery. Most studies on hydraulic heavy mobile machinery coupled with hybrid electric powertrains have taken the energy consumption of the hydraulic pump from experimental data or typical duty cycles. Such data is not readily available and often the work cycles of heavy mobile machinery are not known beforehand. Other studies that have considered detailed hydraulic circuits to simulate the pump loads have simplified the machinery models. Relatively less literature exists where detailed models of complex machinery have been used to capture the loads on complete hydraulic circuit models. There is a lack of multi-domain simulation models accounting for all these aspects in the literature owing to the complexity of heavy mobile machinery applications. It is challenging, and at the same time, essential to consider the interactions between components of various physical behaviour and computation time scales in heavy mobile machinery simulations.
An important approach to couple the dynamics of such multiphysics systems is co-simulation, which has been used in some studies to combine electrical, hydraulic and mechanical components of vehicles [37], [38]. Mechanical, control and electronic disciplines have been incorporated in a functional mock-up interface in the development of a power-shift mechanism co-simulation model for a tractor in [39]. A high fidelity vehicle simulation model was developed in [40] integrating the powertrain, vehicle body, tire, driver and terrain subsystems to evaluate the overall dynamics of a military vehicle as a case example. The soil-tool interaction of the multibody model of a hydraulically actuated wheel-loader implement was also studied in a co-simulation approach in [41]. Mechanical-hydraulic co-simulation models for mobile machinery have been used in [42] and [43] to evaluate control strategies. Methods to improve the stability of multibody-hydraulic co-simulation systems have been proposed in [44]. As one of the case examples an off-road vehicle with soil-wheel interactions was investigated [44]. The hydraulic system of an excavator was coupled with a thermal model of an internal combustion engine in a real-time co-simulation model in [45].
Despite the research efforts explained above, there are limitations in the existing literature. Hybrid electric powertrains, hydraulic actuators and detailed machinery models have been investigated separately. Coupling between them have not included all three together in the reviewed literature. Using co-simulation models, several aspects of heavy mobile machinery have been evaluated, like electro-hydraulic components, interaction of machinery with terrain, control strategies and coupling of hydraulics with the engine. However, even in those studies, the combined effect of coupling the hydraulic actuators with different topologies of hybrid electric powertrains for detailed and complex heavy mobile machinery on deformable terrain has not been investigated. To comprehensively evaluate the performance and powertrain efficiency resulting from the hydridization of heavy mobile machinery it is necessary to develop a model that can comprehensively describe the dynamics of such machinery. This study aims to address this gap in the literature by proposing VOLUME 10, 2022 a co-simulation model that combines various hybrid electric powertrains with detailed heavy mobile machinery mechanical model while capturing the dynamic loads on the machine through interaction with terrain in a virtual environment.
The objective of this paper is to introduce a co-simulation platform for hybrid electric heavy mobile machinery and investigate the effects of hydraulic actuators on the various coupling configurations with the powertrain. The case example of a farm tractor modeled using a semi-recursive multibody formulation [1] is considered in this study. The hydraulic actuation of the front-loader of the tractor is modeled using the lumped fluid theory. In a three-dimensional maneuver of digging and dumping sand particles from one place to another, the hydraulic and traction loads are calculated from the interaction with a deformable terrain. In this study, the series, parallel and series-parallel topologies of the hybrid electric powertrain are considered for coupling FIGURE 1. Co-simulation between multibody dynamics (including hydraulic actuators and deformable terrain) and hybrid powertrain subsystems.
the hydraulic pump with the diesel engine or electric motor. The coupling effect of the hydraulic actuators with the electromechanical interaction of the hybrid electric powertrain configurations are investigated. This paper contains five sections. The rest of the paper is organized such that Sect. II describes co-simulation setup and modeling methods of the subsystems. Sect. III describes the case example of a hybrid electric tractor. Results and discussion are presented in Sect. IV and a conclusion is provided in Sect. V.

II. MODELING METHODS
The multibody dynamic model of the machinery and the hybrid powertrain model form the two subsystems solved separately in the co-simulation setup shown in Fig. 1. The multibody dynamics subsystem consists of the detailed mechanical model of the heavy mobile machinery developed using multibody modeling approach, the hydraulic actuators and the deformable terrain. The hybrid powertrain subsystem includes the user inputs and the hybrid electric powertrain comprising the diesel-electric drive, transmission and controller. Each subsystem is solved separately in its own simulation tool. The subsystems are combined through the update and exchange of variables orchestrated by the communication manager in a discrete event-based co-simulation at designated time scales [46]. The communication manager is contained in the master simulation tool which is the hybrid powertrain subsystem in this co-simulation setup.
The co-simulation begins with the initialization of the TCP/IP socket connection through a successful exchange of handshake data. The dynamics of the system are driven by first providing the user inputs through steering wheel, pedals and joystick. These inputs include the reference torque (T ref ) to drive the machine which is sent to the hybrid powertrain controller, the braking torque (T brake ), steering angle (θ steer ) and the reference voltage signal (U ref ) which are sent to the multibody dynamics subsystem. These user inputs can be provided by a human operator for short as well as large computation time of the co-simulation setup. However, the user experience is enhanced when the computation time is low and synchronized as closely as possible to real-time. T ref is used by the hybrid powertrain to compute the powertrain torque (T pt ) which is sent to the multibody mechanism.
The multibody dynamics subsystem then uses T brake , θ steer and T pt to drive the multibody mechanism model. The hydraulic actuators use U ref along with the piston position (x) and velocity (ẋ) computed by the multibody mechanism to send back the cylinder force (F cyl ). The multibody mechanism also exchanges the friction force (F µ ) and the normal force (F n ) of the tires and hydraulic implements with the deformable terrain to generate soil particles. In return, the hybrid powertrain receives the vehicle velocity (v veh ) and the powertrain rotational speed (ω pt ) from the multibody mechanism, and the pump speed (ω pump ) and pump torque (T pump ) from the hydraulic actuators. In this way the user inputs are supplied by a human operator to produce a customized work cycle in the hybrid powertrain subsystem. Using these inputs the dynamical loads are computed in the mutibody dynamics subsystem through interactions between the mechanical model, hydraulic actuators and the deformable terrain. These loads are sent to the hybrid electric powertrain to investigate the effects on the powertrain components.
A fixed-step discrete solver is used in the master simulation tool to solve the hybrid powertrain subsystem at a time-step of 0.5 ms and to handle the data transfer through the communication manager at a time-step of 1.2 ms [47]. A fourth order Runge-kutta solver is used to solve the multibody dynamics subsystem at a time-step of 1.2 ms [3], [4]. These solver settings are aimed at achieving the lowest possible computation time so that human operators can use the developed model to create customized work cycles for specific applications.

A. MULTIBODY DYNAMICS SUBSYSTEM
The detailed dynamics of heavy mobile machinery, such as underground mine loaders [48], wheel loaders [49], excavators [50], tree harvesters [51] and forklifts [49], can be modeled using the multibody approach. In this approach, the modeling of a generic heavy mobile machinery is usually divided into modeling of mechanical and hydraulic subsystems. The dynamics of the mechanical subsystem can be described using a multibody method such as the augmented Lagrangian-based semi-recursive method [1], whereas the dynamics of the hydraulic subsystem can be described using the lumped fluid theory [52]. The scope of this study is limited to linear hydraulic actuators, where double-acting hydraulic cylinders are considered. It should be noted that the mechanical and hydraulic subsystems are monolithically coupled in this study. Furthermore, a realistic environment such as deformable terrain can be modeled using a combination of the grid-and particle-based methods [4].

1) MULTIBODY MECHANISM MODEL
The equations of motion in semi-recursive multibody methods are formulated first in the Cartesian coordinates and then switched to relative joint coordinates using a velocity transformation matrix. In defining the Cartesian coordinates, the reference point used in this study is at the center of mass of the body. In this method, closed-loop mechanisms are opened using temporary cut-joints and, later on, the corresponding loop-closure constraints are incorporated. Here, the constraints are incorporated using the augmented Lagrangian method [1].
Consider a tree-structure open-loop mechanism of N b rigid bodies with any number of branches as shown in Fig. 2. The open-loop kinematics can be recursively expressed from the base to the leaves or vice-versa using the classical kinematic relations [1], [53]. The Cartesian velocities Y j and accelera-tionsẎ j of the center of mass of body j can be written as where g,ġ, andg are the respective position, velocity, and acceleration of the center of mass, and ω andω are the respective angular velocity and acceleration of the body.
In an open-loop mechanism, Y j andẎ j can be recursively expressed as [54] where B is a transformation matrix, the scalars z,ż, andz are the respective relative joint position, velocity, and acceleration between bodies j − 1 and j, and b and d are the joint type dependent vectors [53]. This method assumes only revolute and prismatic joints because other joints can be represented as their combination using massless intermediate bodies.
Furthermore, the Cartesian coordinates can be mapped onto a set of relative joint coordinates using a velocity trans- where an asterisk (*) denotes the kinematically admissible virtual velocities, 0 is a zero matrix, and the mass matrix M j and the external force vector Q j can be written as [55] VOLUME 10,2022 where m is the body mass, I is an identity matrix, J is the inertia tensor, f is the external force vector, and τ is the external moment vector. Using Eqs. (4) and (5) in Eq. (6), the equations of motion for the open-loop mechanism can be written as where M ∈ R 6N b ×6N b and Q ∈ R 6N b are the composite diagonal mass matrix and external force vector constituting M j and Q j , respectively. For a closed-loop mechanism, a set of N m loop-closure constraints (z) = 0 are incorporated in Eq. (8) using the augmented Lagrangian method [1]. The constraints here are assumed to be holonomic. Accordingly, the equations of motion for the closed-loop mechanism can be written as where α, , and µ are the diagonal matrices containing the respective penalty factors, natural frequencies, and damping ratios, z is the Jacobian matrix of , λ is the vector of iterated Lagrange multipliers, and the terms˙ and¨ can be written as˙ where t is the time derivative of , and˙ z and˙ t are the respective time derivatives of z and t . To solve Eq. (9) without using the algebraic constraint equations (z) = 0, the unknown vector of Lagrange multipliers is computed iteratively as where h = 0, 1, 2, . . . is the iteration step, and λ (0) = 0 for the first iteration. This iterative procedure can be embedded into Eq. (9) to obtain an iterative solution for the accelerations as where M z (0) = Q for the initial iteration.

2) HYDRAULIC ACTUATORS
The pressures in a hydraulic circuit can be modeled using the lumped fluid theory [52]. In this theory, a hydraulic circuit is divided into discrete volumes, where the pressures are assumed to be equally distributed. Thus, the effects of acoustic waves are assumed to be negligible. The modeling concept is shown in Fig. 3, where a hydraulic control volume V s is associated with an incoming flow Q s 1 and outgoing flows Q s 2 and Q s 3 . The hydraulic pressure p s within this hydraulic control volume V s can be computed aṡ whereṗ s is the pressure build-up, B e s is the effective bulk modulus, dV s dt is the change of control volume with respect to time t, Q s k is the incoming/outgoing volume flow rate, and n f is the total number of volume flow rates. The effective bulk modulus B e s in Eq. (14) can be written as where B oil is the bulk modulus of oil, V k is a sub-volume forming the control volume V s , B k is the bulk modulus of the sub-volume V k , and n c is the total number of sub-volumes. A directional control valve can be modeled using a semi-empirical modeling method [56]. The volume flow rate Q d through a directional control valve can be written as where C v d is the semi-empirical flow rate constant of the valve, sgn (·) is the signum function, p is the pressure difference, and U is the relative spool position (measured in volts) that can be computed aṡ whereU is the time derivative of U , U ref is the reference voltage signal for the reference spool position, and τ is the time constant. The value of C v d can be obtained from the manufacturer's catalogue and the value of τ can be obtained from the Bode-diagram of the valve describing the valve spool dynamics. In a double-acting hydraulic cylinder as shown in Fig. 4, the volume flow rates can be expressed as where Q in and Q out are the respective incoming and outgoing volume flow rates, A a and A b are the respective areas on the piston and piston-rod sides, and x andẋ are the respective piston position and piston velocity. The force F cyl produced by a double-acting hydraulic cylinder can be written as where p a and p b are the respective chamber pressures on the piston and piston-rod sides, and F µ is the total friction caused by sealing [57]. Hydraulic pumps are typically used to convert the mechanical energy from diesel engines or electric motors into hydraulic energy. The volumetric flow rate of the pump can be expressed as [58]: where α is the ratio between the current and maximum displacement of the pump, ω p is the rotational speed of the pump, D max is the maximum displacement of the pump, and Q l is the leakage flow rate. The load torque required by the pump shaft can be calculated as [58]: where p p is the pressure difference across the pump. Depending on the coupling the diesel engine or the electric motor has to provide this torque to operate the hydraulics. The power of the pump is: The load-sensing controller setting pressure for the pump output is determined as: where p LS is the pressure in the hydraulic control volume that the load-sensing pump controller is connected to, p ref is the reference pressure that must be kept above the load-sensing pressure, and p sb is the stand-by pressure. The derivative of the pump displacement ratio is then calculated as: where k p is the pump controller gain, p p is the working pressure of the pump and τ p is the pump controller time constant.

3) DEFORMABLE TERRAIN
Deformable terrain can be modeled using a combination of the grid-and particle-based methods [4]. In this hybrid method, the vertical component of the applied force compresses and displaces land. The horizontal force component generates sand particles as shown in Fig. 5. Consider a cellular automata on an L × L square grid such that h(i, j) is the system height at the center of a cell (i, j) and h f (i, j) is the height shift because of an applied vertical force [59], [60]. Here, the terrain material is assumed to be a cohesion-less material with low compressibility such as dry sand. The total height h (i, j) of a column can be written as The displacement of sand field is triggered when the gradient angle of h (i, j) is greater than the pile slope of 35 • in this study. Accordingly, the automata state can be updated as where z + is the velocity of flowing matter, and h x and h y are the respective partial derivative of h (x, y). Avalanches are also considered in this method. Sand particles are generated when the horizontal component of the applied force is greater than the shear impulse limit. Here, sand particles are simulated as rigid spherical bodies with six degrees of freedom containing both mass and inertia properties. The dynamics of the sand particles are described using the augmented Lagrangian-based semi-recursive multibody method presented in Sect. II-A1, whereas their rotations are described using Euler parameters. The size of the granular sand particles is determined based on the limit of 1000 particles. In sand particle representation, the degree of sand compaction is incorporated by considering a void ratio for all the free particles.
Sand particles in contact with the sand heightfield merge with the sand grid when their linear and angular velocities are below a threshold value [61]. The merging and generation of sand particles are performed using a volume update, where sand particles are replaced by the corresponding static heightfield and vice versa, respectively. More details about the modeling of deformable terrain can be followed from the previous publication [4] of the authors.

B. HYBRID POWERTRAIN SUBSYSTEM
In a conventional heavy mobile machinery powertrain with hydraulics the diesel engine provides the power for both the traction as well as the hydraulic loads. When the powertrain is electrified, the power for these loads can be provided by the VOLUME 10, 2022 diesel engine, electric energy storage or both. There are three basic topologies based on power flow -series, parallel and series-parallel. In the series architecture shown in Fig. 6a, the diesel engine is decoupled from the traction load and is connected directly to the hydraulics. In the parallel architecture shown in Fig. 6b, the diesel engine is completely decoupled from the traction load and the hydraulics which are both driven by individual electric motors. This allows the engine to be operated at its optimal efficiency. However, when the engine is decoupled from one or more loads there is a double energy conversion required from mechanical to electric and then back to mechanical which can cause additional losses. In a series-parallel architecture, double energy conversion are avoided by using a power-split device so that the engine can support the traction demand also. In this architecture the hydraulics may be coupled directly to the engine as shown in Fig. 6c to avoid double energy conversion or it may be coupled to a second motor as shown in Fig. 6d to have better control of the engine operating points. [5]

1) DIESEL-ELECTRIC DRIVE
A diesel engine is used as the primary energy source of the powertrains in this study. A maximum torque curve limits the torque of the engine while its efficiency is looked up from an efficiency map. A first-order transfer function is used to model the dynamics of the engine. The fuel consumption of the engine can be expressed as [62]: where T d is the engine torque, E J is the fuel heat value, η d is the efficiency of the engine, and ω d is the engine speed which can be computed as: where J d is the rotational inertia of the engine shaft and T load is the load torque on the engine shaft. The secondary energy source used in this study is a generic rechargeable battery described by the Shepherd equations [63]. The battery voltage at the output terminals can be expressed as: where E 0 is the constant voltage value, R is the internal resistance, I bat is the battery current, E bat is the voltage at the output terminals, K is the polarization constant/polarization resistance, Q max is the maximum battery capacity, i bat* is the filtered current, A is the exponential zone amplitude, and B is the exponential zone time constant inverse. These parameters are extracted from the manufacturer-provided discharge curve. The actual charge, Q, can be expressed as: The state of charge (SOC) of the battery can be estimated using the modified Coulomb counting method as: where SOC init is the initial SOC.
A quasi-static model of the electric machine has been used in this study to replicate the losses and to reduce computation time while achieving the required accuracy [64]. Similar to the diesel engine a maximum torque curve combined with a first-order transfer function limits the electric machine torque and its changes while the efficiency is looked up from an efficiency map. The power of the electric machine flowing to or from the battery can be calculated as [65]: where η EM is the efficiency of the electric machine, ω EM the rotational speed of the electric machine and T EM is the torque of the electric machine.

2) TRANSMISSION
The gearbox can consist of clutches and simple gears controlled such that the gears are changed based on the vehicle speed. The frictional torque of a disc friction clutch can be written as [66]: where µ cl is the coefficient of friction of the clutch plate, n is the number of friction surfaces, r eff is the effective torque radius, P is the engagement pressure, and A is the engagement surface area. The clutch is unlocked when the clutch torque is zero, slipping when the clutch torque is above zero but less than the input torque and locked when the clutch torque reaches the input torque [66]. The transmission ratio of the gears allow transfer of torque and speed such that [65]: where ω in is the speed input to the gear, i is the gear ratio, ω pt is the speed output from the gear or powertrain, η g is the transmission efficiency, T pt is the torque output from the gear or powertrain and T in is the torque input to the gear. Planetary gearsets can be used in power-split devices to direct the energy coming from the different diesel engine and battery to the wheels and the hydraulic pump as needed. The relations between the torque and speed of the planet carrier, sun and ring gears of the planetary gearsets can be written as [65]: where ρ is the gear ratio between the sun and the ring, ω d is the angular velocity of the carrier which is the angular velocity of the diesel engine connected to it, ω S is the angular velocity of the sun, ω R is the angular velocity of the ring, T C is the torque of the carrier, T S is the torque of the sun, T R is the torque of the ring, η RC is the efficiency of power transfer between the ring and the carrier, and η RS is the efficiency of power transfer between the ring and the sun.

3) CONTROLLER
The co-simulation platform requires the powertrain controller to be able to perform with low computational load and without prior knowledge of the work cycle. A rule-based supervisory control (Fig. 7) based on heuristics has been used because it is effective in real-time implementation without the knowledge of a pre-determined work cycle [65]. The rules of the controller are set using both thermostat and power follower heuristic design principles [67], leading to the following modes: 1) Mode 1 -The machinery starts its operation with mode 1 where the engine provides a constant minimum power (P d_min ). When the battery state of charge (SOC) drops below its lower threshold (SOC low ) or the total power demand (P load ) exceeds the maximum allowable battery power (P b_max ) the control is switched to mode 2. 2) Mode 2 -When SOC < SOC low or P load > P b_max the controller uses a thermostat strategy to level the load. The engine is commanded to run at a constant optimal power (P d_opt ) to charge the battery while the battery acts as the equalizer. When the SOC increases above its upper threshold (SOC upp ) and P load < P d_min the mode changes back to 1. 3) Mode 3 -When P load > P b_max the controller uses a power follower strategy so that the engine can follow the load to avoid harsh battery usage. In this mode the engine is commanded to provide P load . The mode changes back to 2 only when P load < P b_max . The values of P d_min , P d_opt , SOC low , SOC upp and P b_max can be tuned to achieve the desired output. In this study, the value of P d_min is 25 kW, P d_opt is 78.5 kW, SOC low is 50%, SOC upp is 80%, and P b_max is 158 kW. In addition to the reference power dictated by the controller (P d_ref ), the engine also provides the power consumed by the hydraulic pump if it is directly coupled to the diesel engine.

III. CASE STUDY OF A TRACTOR MODEL
In this study, the case example of a tractor modeled using the semi-recursive multibody formulation explained in Sect. II-A1 has been combined with the four powertrain configurations shown in Fig. 6. The specifications of the mechanical model of the tractor are based on the Valtra N174 tractor. The specifications of the hybrid electric powertrain components used in the four models are provided in Table 1.  The front loader is operated by four hydraulic cylinders connected to two directional control vales as shown in Fig. 8, modeled as explained in Sect. II-A2. The piston side of each cylinder is denoted as port A while the rod side is denoted as port B. The valves are connected to a hydraulic pump which is coupled with the diesel engine or electric motor according to the configurations presented in Fig. 6.
A digging and dumping operation of the tractor is simulated in a virtual environment shown in Fig. 9. In the beginning of the work cycle, the front loader is first aligned as required for the digging operation. The tractor then begins to move forward towards the mound of soil, performs the digging operation and reverses back as shown in Fig. 9. It then proceeds to transport the sand in an uphill maneuver to the dumping location. Once the dumping operation is completed  and the front loader is lowered, the cycle ends. The speed of the tractor and the mass of the load in the bucket is shown in Fig. 10. The mass of soil picked up by each model is not identical due to the chaotic movement of the soil particles in the deformable terrain. However, the mass varies within an acceptable range with minor differences. The specifications of the computer used for the simulation are provided in Table 2.

IV. RESULTS AND DISCUSSION
This section presents hydraulic actuator performance and energy efficiency results of the hybrid electric tractor model under the digging and dumping work cycle. The hydraulic actuator forces in the cylinders are shown in Fig. 11 and the pressures are shown in Fig. 12. The fuel consumption, diesel engine power, battery SOC, traction motor power and hydraulic pump power are shown in Fig. 13. The diesel engine, traction motor and hydraulic pump energy consumption are shown in Fig. 14. The integration time taken by the simulation models is shown in Fig. 15.

A. HYDRAULIC ACTUATOR FORCES AND PRESSURE
The actuator force plot (Fig. 11) indicates identical responses for both lift and tilt cylinders. However, after 15 seconds, the actuator forces begin to rise, denoting the initiation of the tractor digging operation. The multibody dynamic system, the hydraulic system and the deformable terrain model start interacting here. A high fluctuation and peak force of about 50 kN for the lift cylinder and 20 kN for the tilt cylinder can be observed here due to the harsh collision between the bucket and the soil particles. A resulting fluctuation in the pressure of the lift and tilt cylinders can also be seen. After the digging operation the actuator forces lower to about 40 kN for the lift cylinder and 10kN for the tilt cylinder. This force is required to support the weight of the soil particles being carried by the bucket. A slight drop in actuator forces is observed at around 35 seconds which is a result of some soil falling off from the bucket during the reverse maneuver (Fig. 10). At around 45 seconds, a sudden increase in pressure in all sides of the cylinders is observed, followed by steep drops in the piston side pressures of all the cylinders at about 50 seconds. This  variation in pressure is brought about by the lowering of the front loader for transportation. The actuator forces fluctuate further during the uphill transportation due to the gradient angle. At around 70 seconds, the actuator forces are observed to decline as the soil is dumped. Variation in pressure is observed again in Fig. 12 as the front loader is lowered at the end of the work cycle.

B. FUEL CONSUMPTION, POWER AND SOC
As seen in Fig. 13(a), the highest amount of fuel is consumed by the series model, followed by the series-parallel, parallel and series-parallel electric models, in that order. There is a 1.5% fuel-saving in the series-parallel model compared to the series as listed in Table 3. This comes from the marginally lower diesel power of the series-parallel model seen in the diesel power plot as shown in Fig. 13(b). As listed in Table 3, the parallel and series-parallel models accrue fuel savings of 7.3% and 10.5%, respectively, compared to the series due to further lower diesel power consumption. This shows that, some fuel savings can be made by using the series-parallel topology which distributes the traction load between the motor and the diesel engine. These savings in the series-parallel topology are possibly achieved by reducing the double energy conversion that occur in the series topology. Furthermore, it is seen from Fig. 13(a) that coupling the hydraulics directly to a second motor (parallel and series-parallel electric) results in lower fuel consumption than coupling it to the engine. In heavy mobile machinery applications which involve long durations of intensive work, such improvements can amount to significant fuel-savings.
The diesel power plot in Fig. 13(b) reflects the control strategy behavior as well as the coupling of the hydraulics with the energy sources. The work cycle begins with low engine power because the controller operates at mode 1 where the engine tries to provide the minimum engine power. In the case of the series and series-parallel models, the engine power is higher than the other models in this region because the engine also provides the pump power due to direct coupling of the hydraulics with the engine. At around 12 seconds, the engine is commanded by the controller to run at the total load power because the traction motor power increases beyond the maximum battery power resulting in mode 3. The engine power steeply increases here to follow the load according to the mode 3 rule. Although the diesel engine attempts to run at the load power, it must be noted that the engine does not individually supply the load demands. The engine power is directed to the battery, hydraulic pump and to the transmission depending on the powertrain configuration. At around 17 seconds, the diesel power of the parallel and the series-parallel electric models start to decline as the traction motor power subsides and the controller shifts to mode 2. On the other hand, at the same time, the diesel power of the series and series-parallel models continue to rise because the engine provides the power for the coupled hydraulics. Following this peak, the diesel power comes down to the optimal diesel power for all the models as determined by the control logic (Fig. 7) due to lower traction and hydraulic loads. Similar patterns are seen in the rest of the work cycle where the diesel power can be co-related to the traction motor power and the hydraulic pump power depending on the coupling. Therefore, it can be summarized that relatively low diesel power is observed when the hydraulic load is high for parallel and series-parallel electric models because the hydraulic pump is uncoupled with the engine. Relatively low diesel power is also observed when the traction load is high for the series-parallel models because the traction load is distributed between the diesel engine and the traction motor.
The battery SOC is seen to increase and decrease in a similar manner for all four models (Fig. 13(c)). In mode 1 the SOC stays almost at the same level of 80% for all models because both the electric loads and the charging from diesel power are minimal. Mode 3 is triggered by high load demands which are primarily provided by the electric motor. Although the engine follows the load demand to charge the battery at a higher rate during this mode, the power consumption of the electric motor is higher than the charging power provided by the engine, as seen in Fig. 12. Therefore, the SOC declines steeply whenever mode 3 occurs because of the high motor power. During mode 2 the SOC replenishes as the electric loads are low and the diesel power can charge the battery. For the first half of the work cycle, the series-parallel models maintain a relatively higher SOC level than the other topologies due to lower traction motor power. However, when mode 2 occurs for a second time, the SOC of these models declines steeply as the diesel power is not high enough to charge the battery more. The SOC of the series-parallel model is eventually replenished to 73.7% which is higher than the other models due to lower electric loads for both traction and hydraulic demands. Whereas, the SOC for the series-parallel electric model acquires the lowest final value of 71.0% because the hydraulic power is supplied by a second motor, depleting the SOC relatively more. Also, the diesel power for this model is lower which causes relatively slow charging. The series and parallel models end up in a middle range with 73.4% and 72.6% final SOC, respectively. The SOC of the parallel model is lower than that of the series model because of the hydraulic power extracted from a second motor. All four models attain final SOC in the range of 70-75% which is appropriate for the simulated case. The charging or discharging characteristics can be altered as per the requirements of the machinery by tuning the parameters of the controller. For example, a charge depleting profile may be desirable for machines that have plug-in hybrid options, while a charge sustaining profile may be desirable if it is not a plug-in hybrid. Such profiles can be attained by adjusting the parameters.
The traction motor power shown in Fig. 13(d) can be related to the tractor speed (Fig. 10). High motor power is observed in three instances of the work cycle. These peaks correspond to the forward, reverse and uphill maneuvers as the user increases the reference torque input through the accelerator pedal. As the pedal is fully pressed the traction motor power reaches the peak value of 600 W. When the pedal is fully released the motor power drops to zero as the user intends the tractor to be stationary. The traction motor power behavior of the series model is identical to that of the parallel, and that of the series-parallel model is identical to that of the series-parallel electric. This is because these models have similar topology in terms of fulfillment of the traction demand. The motor power is lower for the series-parallel models compared to the series and the parallel because in this topology the traction load is distributed among the diesel engine and the electric motor. Whenever the motor power increases, the increased diesel power supports the traction demand, thus causing the motor power to decline. This displays the advantage of the series parallel topology where double energy conversion are minimized.
The pump power (Fig. 12(e)) is seen to rise whenever the front loader is actuated. This includes the alignment of the loader and the digging and dumping operations. In the cases of the series and series-parallel models, the increase in pump power leads to increased diesel power because of the direct coupling of the pump to the engine. In the cases of the parallel and series-parallel electric models, the pump power is extracted from the battery through the second motor to which the pump is coupled. This depletes the SOC by a small amount. However, the battery charging power simultaneously provided by the engine is higher, resulting in an overall increase in SOC.

C. ENERGY CONSUMPTION
The energy consumption of the four models is shown in Fig. 14. The series and the series-parallel models have the highest diesel energy consumption (Fig. 14(a)) because of the direct coupling of hydraulics to the engine. However, the diesel energy consumption of the parallel and the series-parallel electric models are the lowest because the hydraulic energy comes from electric sources. In the traction motor energy plot in Fig. 14(b), steep increases are observed for all models in mode 3 where the traction demand of the work cycle is high. The series-parallel topology shows lower traction motor energy consumption than the other models owing to the lower motor power consumption seen in Fig. 14(b). As explained in Sec. IV-B, the motor power consumption is lower in the series-parallel topology due to the distribution of the traction load between the engine and the electric motor. Therefore, corresponding lower motor energy consumption seen in Fig. 14(b) substantiates the advantage of the series-parallel topology. The pump energy (Fig. 14(c)) increases during modes 1 and 2 as the hydraulic actuation takes place when the tractor is stationary in this work cycle. The pump energy (Fig. 14(c)) is reflected in the higher diesel energy consumption of the series and series-parallel models due to the direct coupling of the engine and pump.
Even though the fuel consumption of the parallel model is 7.3% lower than that of the series model, the relative savings in total energy consumption is only 2.9% as listed in Table 3. This is because the energy consumed by the second motor to support the hydraulic load causes double energy conversion lowering the overall energy savings. It can be seen in Table 3 that the series-parallel model reduces the total energy consumption by 8.4% compared to the series model. In both of these models the hydraulic pump is coupled to the engine directly, eliminating the possibility of double energy conversion. However, the series-parallel topology saves energy by distributing the traction load between the diesel engine and electric motor. The series-parallel electric model lowers the total energy consumption by 12.3% (Table 3) compared to the series model. In this model there is a possibility of similar energy losses due to double energy conversion as in the case of the parallel model. In the simulated work cycle the energy consumed by the pump constitutes a smaller portion of the entire load because the traction demand is high. However, depending on the application the proportion of traction and hydraulic loads may vary largely in heavy mobile machinery.

D. REAL-TIME CAPABILITY
The real-time capability of the four models is shown in the integration time plot in Fig. 15. The integration time of all the models is below the time-step of the multibody subsystem which is 1.2 ms, for the parts of the work cycle where the tractor performs only transportation work. However, when the soil is dug at around 17 seconds, the integration time of all the models peak to a value of about 1.6 ms. This occurs due to the generation of soil particles which requires high computation time. During this period the models are not real-time capable because the integration time is higher than 1.2 ms and therefore, the simulation time is not synchronized to realtime. When the soil particles are returned to the ground during the dumping at the end of the work cycle, the integration time is, once again, low. This implies that the assimilation of the soil particles back into the ground is computationally efficient and synchronized to real-time.
The series model is the simplest model of the four and therefore takes marginally lower computation time. The series-parallel topology-based models have an additional power-split device, while the parallel and series-parallel electric models have an additional electric motor. Therefore, these models have slightly higher complexity compared to the series model. However, the patterns of integration time taken by all the models do not vary much even with varying complexities of the models. While the models can be simulated in real-time for the most part of the work cycle, the soil particle generation can be improved to reduce the computation time to make the model fully real-time capable. Simulations closer to real-time result in a more realistic driving experience for the human operator while creating work cycles specific to the working conditions of the machinery. However, the efficiency of the simulation achieved in this study is adequate for human operators to drive the tractor model.

E. DISCUSSION
The dynamic performance of hydraulic actuators caused by the loads computed by the detailed physics-based model of the multibody tractor and deformable terrain is shown in Sec. IV-A. Using the presented simulation models, these loads can be simulated by human operators to take into account the diversity in heavy mobile machinery work cycles instead of using experimental data or typical duty cycles. The interactions between the detailed machinery model, hydraulic actuators, deformable terrain and hybrid electric powertrains were captured in the developed models through a cosimulation approach. The co-simulation setup was designed to maintain complexity and simplicity where required while handling the time-scales of the components based on the subsystem requirements.
The effects of coupling the hydraulic actuators with the different power sources were investigated using the four powertrain configurations shown in Fig. 6. Higher fuel savings were achieved by the models in which the hydraulics are coupled to a second electric motor. However, the electric coupling of hydraulics alone in the parallel model did not result in significant savings in total energy consumption because of double energy conversions. The series-parallel models showed improvements in total energy consumption by distributing the traction demand between the diesel engine and the traction motor. The highest fuel and total energy savings were achieved by the electric coupling of hydraulics in the series-parallel topology.
To determine the suitability of the configuration and hydraulics coupling, various factors must be considered. If higher fuel savings is the goal, then the series-parallel electric model displays the best results. The losses resulting from double energy conversion can be considered minimal and may be agreeable if, for example, the electric energy comes from clean sources in the case of plug-in hybrids.
However, if all the electric energy in the battery comes from the diesel engine charging, then it may be better to couple the hydraulics directly to the engine. The manufacturability aspect also needs to be considered because the gains from these configurations come at the cost of adding additional components.

V. CONCLUSION
This study introduced a co-simulation platform for a detailed multibody model of hydraulically actuated heavy mobile machinery in a deformable terrain environment, coupled with different power sources of hybrid electric powertrains. The case study of a farm tractor equipped with a hydraulically actuated front loader with four powertrain configurations were considered -a series topology with the hydraulics directly coupled to the engine, a parallel topology with hydraulics coupled to an electric motor, a series-parallel topology with hydraulics directly coupled to the engine and another series-parallel topology with hydraulics coupled to an electric motor.
In a user-defined work cycle of soil digging and dumping operation the dynamic performance of the hydraulic actuators was studied from the cylinder forces and pressures. The pressure and forces in the hydraulic actuator cylinders showed how the front loader was operated by the user using a control voltage signal and dynamic hydraulic loads were simulated through interaction with the soil. Fuel savings of 7.3% and 10.5% were achieved by the parallel and the seriesparallel models, through the coupling of the hydraulics to the electric motor instead of the engine. The models with series-parallel topology showed improvements in total energy consumption of 8.4% and 12.3% compared to the series topology by sharing the traction demand between the diesel engine and the electric motor. The highest fuel and total energy savings were observed in the series-parallel model with electric coupling of hydraulics. This demonstrated the benefits of the series-parallel topology and coupling of the hydraulics to an electric motor in the simulated case study. The study shows that the developed simulation method can be used to couple multi-domain simulation models of heavy mobile machinery capturing complexities where required. Such simulation models can be effectively used to compare and select powertrain concepts where both hydraulic and traction loads are considered from user-defined work cycles without the need for any measurement data.
The results obtained in this study have not been validated because of a lack of experimental measurements. However, the behavior of the system and the individual components are logical. The energy consumption characteristics shown in the results have a good agreement with those observed in existing literature. It should be noted that the purpose of the co-simulation platform developed in this study is to compare various coupling configurations of the powertrain energy sources and the hydraulics. The developed models do not generate exact values of energy consumption, but they aid the relative comparison between various configurations. However, in the future, the models can be improved further to include detailed modeling of the components and working environment to produce more realistic results.
In future work, the soil particle modeling can be improved to reduce computation time so that the model can be real-time capable during digging operations. This will help in improving the user experience such that the effect of different human operators on the energy efficiency of hybrid powertrains can be studied. The influence of the deformable terrain on the operation and energy consumption characteristics of the machine can also be studied. Global optimal solution can be implemented in the controller when the work cycle is fixed. Data from conventional machinery can also be utilized in defining the controller rules. Coupling the hydraulics to electric energy sources presents opportunities for energy recovery which can be studied further. The multibody modeling approach also allows the developed models to be utilized to analyze various product processes such as operator training and failure analysis along with the possibility to model other complex hybrid heavy machinery.