On the cosimulation of multibody systems and hydraulic dynamics

Abstract

The simulation of mechanical devices using multibody system dynamics (MBS) algorithms frequently requires the consideration of their interaction with components of a different physical nature, such as electronics, hydraulics, or thermodynamics. An increasingly popular way to perform this task is through co-simulation, that is, assigning a tailored formulation and solver to each subsystem in the application under study and then coupling their integration processes via the discrete-time exchange of coupling variables during runtime. Co-simulation makes it possible to deal with complex engineering applications in a modular and effective way. On the other hand, subsystem coupling can be carried out in a wide variety of ways, which brings about the need to select appropriate coupling schemes and simulation options to ensure that the numerical integration remains stable and accurate. In this work, the co-simulation of hydraulically actuated mechanical systems via noniterative, Jacobi-scheme co-simulation is addressed. The effect of selecting different co-simulation configuration options and parameters on the accuracy and stability of the numerical integration was assessed by means of representative numerical examples.

1 Introduction

During the last decades, multibody-based simulation tools have established themselves in industry as they allow rapid testing of mechanical systems without costly prototyping. With the increase in computational power at affordable prices, these tools have also achieved the real-time limit outside the realm of academic discussions; see, for example, [1, 2]. However, in parallel with the maturity of this technology, the demand has risen for system-level simulation tools. Most engineering applications of industrial interest include not only mechanical components whose behavior can be described with conventional multibody system dynamics (MBS) formulations, but also elements with different physical properties, behavior, and time scale. These include control and electronics components, hydraulics circuits, and thermodynamic subsystems, among many others. When this is the case, the dynamics of the multibody components can no longer be considered on its own; instead, the interaction with the rest of elements in the engineering assembly must be appropriately described and taken into account in the simulation. For some of these applications, it is possible to develop monolithic formulations that are capable of describing the overall system dynamics with a common set of equations [3].

Monolithic methods are often accurate and efficient but require the modeling details of every subsystem to be accessible by the solution solver. This poses a problem in industrial applications in which at least some of the information about the system components must remain undisclosed to safeguard intellectual property rights. Moreover, in some applications the different physical behavior and time scales of the subsystems can make it challenging to develop a dynamics formulation able to represent all the phenomena involved. Assigning a particular solver to each subsystem, on the other hand, makes it possible to tailor its equations and implementation to the nature of the component that it is meant to represent. Although it would be, in theory, possible to build an all-encompassing simulation tool to suit all the engineering needs, in practice, this may not always be the case, and, in addition, the use of tailored modeling and solver tools for each engineering field is often advantageous. Co-simulation is an alternative to monolithic methods, which addresses these issues by assigning a solver to each subsystem in an application and then coupling their solution processes rather than developing a unified model of the whole system under study [4]. The subsystems exchange information through coupling variables, that is, input and output variables, at communication points in a discrete-time fashion. The result is a modular simulation environment, in which the implementation details of each subsystem remain hidden from the others. Additionally, co-simulation setups enable the distribution of the computational workload between several processing units [5, 6], even over network connections, and thus they can be used to improve the computational efficiency.

One of the main disadvantages of co-simulation with respect to its monolithic counterparts is the introduction of discontinuities and delays caused by the discrete-time nature of the communication interface [7]. These can eventually render the numerical integration process unstable in some cases, especially when noniterative co-simulation schemes are adopted [8]. Iterative co-simulation methods can be used to alleviate this problem; see, for example, [9, 10]; their use, however, becomes problematic in certain applications. This is the case of some real-time environments, in which the available time to carry out the computations is limited, and of systems that include components that cannot retake their integration steps, for example, physical components in hardware- or system-in-the-loop settings. In these cases, noniterative co-simulation schemes must be used.

The stability of noniterative co-simulation schemes can be improved in a number of ways. Polynomial extrapolation and approximation are among the most popular and simplest approaches [11,12,13,14]. The only information they require is the time history of the coupling variables exchanged at the co-simulation interface; moreover, they are relatively easy to implement. Further analysis of the evolution of the coupling variables has led to the development of methods to improve stability based on the energy exchanges at the interface [15, 16] and to the determination of instantaneous frequencies of the subsystems in the setup [17]. The availability of additional information about the behavior and properties of the subsystems enables the use of stabilization techniques that can be particularized to match the needs of particular problems [18, 19]. It is also possible to estimate these properties from the information contained in the coupling variables to develop correction algorithms [20].

Co-simulation has been applied successfully to a wide range of applications that involve mechanical systems in multiphysics contexts, such as automotive [12], railway [21, 22], biomechanics [23], and heavy machinery [24] simulation. Hydraulically actuated mechanisms are a particularly interesting case of these, as they are present in a large number of industry and research applications. Several monolithic approaches to solve together the dynamics of mechanical and hydraulic components can be found in the literature [25,26,27]. Hydraulics and mechanics have also been coupled by means of co-simulation schemes [24, 25] and publications exist that deal with particular applications of this approach [28, 29]. However, the systematic study of the coupling method and the effect of selecting a given co-simulation scheme and set of coupling parameters on the efficiency and accuracy of the numerical integration has received comparatively less attention [30].

The objective of the present paper is to gain insight into the co-simulation of multibody and hydraulic dynamics when noniterative Jacobi schemes are used to couple the subsystems. Although there exist other noniterative schemes, such as Gauss–Seidel, in this study, we have chosen the Jacobi scheme for its efficiency potential. Moreover, explicit Jacobi schemes are often the option of choice in industrial applications with demanding computational requirements, such as those that couple physical and virtual components, that is, Hardware-in-the-Loop and System-in-the-Loop setups. To this end, we carry out a benchmark study. We compare the influence on the solution and coupling accuracy of two sets of coupling variables, two approaches to extrapolate the coupling variables, two multibody dynamics formulations with different internal integrators, and three solvers for the hydraulic system. We take the largest achievable communication step-size as an indicator of the simulation efficiency; in this work, we adjust its value to match the integration step-size of the slow subsystem, that is, the mechanical one. It can be expected that the computational overhead associated with information exchanges at communication points in practical applications will be reduced when the size of the communication step increases. On the other hand, longer step-sizes decrease the computational load required to integrate the multibody system dynamics. The maximum achievable step-size required to integrate the hydraulic subsystem is also touched upon, since it can also be expected to have an effect on system-level efficiency. To limit the number of considered variables, here we consider only systems with constant integration step-sizes and matching time grids.

We use two representative examples of hydraulically actuated mechanical systems to illustrate the effect of the above-mentioned co-simulation configuration options on the overall system behavior and coupling accuracy. The results showed that the accuracy of Jacobi-scheme co-simulation of mechanical and hydraulic systems can be enhanced by a proper selection of coupling variables and integration schemes.

2 Simulation of multibody and hydraulic dynamics

For the purposes of this study, we implemented a generic co-simulation manager and relevant subsystems. Here we only briefly introduce the multibody and hydraulic dynamics, whereas the co-simulation scheme and coupling accuracy evaluation are covered in more detail.

2.1 Multibody formalisms

A mechanical system can be described with a set of $n$ generalized coordinates $\mathbf {q}$ subjected to $m$ kinematic constraints $\mathbf{\Phi }=\mathbf{0}$, which in this work are assumed to be holonomic for simplicity. Without loss of generality, we assume that the generalized velocities are the derivatives with respect to time of the coordinates $\mathbf {q}$, that is, $\mathbf{v} = \dot {\mathbf {q}}$. The equations of motion can then be expressed as

$$\begin{aligned} \mathbf {M}\ddot {\mathbf {q}}+ \mathbf {c}&= \mathbf {Q}+ \mathbf {Q}_{c}, \end{aligned}$$

(1a)

$$\mathbf{\Phi }(\mathbf{q},t)=\mathbf{0},$$

(1b)

where $\mathbf {M}$ is the $n \times n$ mass matrix, $\mathbf {c}$ contains the Coriolis and centrifugal forces, and $\mathbf {Q}$ and $\mathbf {Q}_{c}$ are the applied and constraint forces, respectively. The generalized constraint reactions can be expressed as ${\mathbf{Q}}_c=-{\mathbf{\Phi }}_{\mathbf{q}}^{\mathrm{T}}\mathit{\lambda }$, where ${\mathbf{\Phi }}_{\mathbf{q}}=\partial \mathbf{\Phi }/\partial \mathbf{q}$ is the $m \times n$ Jacobian matrix of the constraints, and $\mathit{\lambda }$ contains $m$ Lagrange multipliers.

The dynamics in (1a), (1b) are described by a system of differential-algebraic equations (DAEs). Several algorithms have been proposed to perform the numerical solution of such sets of equations. One of them is the index-1 augmented Lagrangian algorithm (ALI1), introduced in [31] and further developed in [32]. This method enforces the satisfaction of kinematic constraints via an iterative algorithm, in which both the Lagrange multipliers and the system accelerations $\ddot {\mathbf {q}}$ are consecutively updated:

$$\begin{array}{rl}& (\mathbf{M}+{\mathbf{\Phi }}_{\mathbf{q}}^{\mathrm{T}}\alpha {\mathbf{\Phi }}_{\mathbf{q}})\ddot{\mathbf{q}}+{\mathbf{\Phi }}_{\mathbf{q}}^{\mathrm{T}}{\mathit{\lambda }}^{\ast }\\ & =\mathbf{Q}-\mathbf{c}-{\mathbf{\Phi }}_{\mathbf{q}}^{\mathrm{T}}\alpha ({\dot{\mathbf{\Phi }}}_{\mathbf{q}}\dot{\mathbf{q}}+{\dot{\mathit{\Phi }}}_t+2\xi \omega \dot{\mathbf{\Phi }}+{\omega }^2\mathbf{\Phi }).\end{array}$$

(2a)

$${\mathit{\lambda }}_{j+1}^{\ast }={\mathit{\lambda }}_j^{\ast }+\alpha (\ddot{\mathbf{\Phi }}+2\xi \omega \dot{\mathbf{\Phi }}+{\omega }^2\mathbf{\Phi }),$$

(2b)

where ${\mathit{\lambda }}^{\ast }$ are the $m$ modified Lagrange multipliers, and the subscript $j$ stands for the iteration number. The modified Lagrange multipliers are made proportional via a penalty factor $\alpha$ to the violation of kinematic constraints at the configuration, velocity, and acceleration levels, following an approach similar to Baumgarte stabilization [33], in which $\omega$ and $\xi$ are scalar parameters.

The ALI1 formalism can be used together with a given numerical integrator, in this study a semi-implicit forward Euler formula, following a fixed-point integration procedure to conduct the forward-dynamics simulation of the system motion. If the state of a mechanism at time-step $n$, $\mathbf {q}_{n}$ and $\dot {\mathbf {q}}_{n}$, is known, then Eqs. (2a), (2b) can be used to obtain the accelerations $\ddot {\mathbf {q}}_{n}$. These can be later integrated to obtain the positions and velocities in time-step $n+1$.

On the other hand, replacing the fixed-point scheme with a Newton–Raphson solution approach can contribute to the improvement of the efficiency and stability of the solution process [34]. This requires the merging of the integration formulas, for example, the Newmark integrator family, and the dynamics equations, the latter of which now take the form

$$\mathbf{M}\ddot{\mathbf{q}}+{\mathbf{\Phi }}_{\mathbf{q}}^{\mathrm{T}}\alpha \mathbf{\Phi }+{\mathbf{\Phi }}_{\mathbf{q}}^{\mathrm{T}}{\mathit{\lambda }}^{\ast }=\mathbf{Q}.$$

(3)

Introducing the equations of the integrator in the dynamics equations (3) and establishing the equilibrium at time-step $n+1$, we obtain a system of nonlinear equations

$$ {\mathbf{g}}\left ( \mathbf {q}, \; \dot {\mathbf {q}}\right ) = \mathbf {0}. $$

(4)

The solution of this nonlinear system of equations can be found through Newton–Raphson iteration:

$$ \begin{aligned} &\left [ \dfrac{d \mathbf{g} ( \mathbf {q}, \dot {\mathbf {q}}) }{d \mathbf {q}} \right ]_{j} \Delta \mathbf {q}_{j+1} = - \left [\mathbf{g} (\mathbf {q}, \dot {\mathbf {q}})\right ]_{j}, \\ &\mathbf {q}_{j+1} = \mathbf {q}_{j} + \Delta \mathbf {q}_{j+1}. \end{aligned} $$

(5)

Mass-damping-stiffness orthogonal projections are used to ensure the satisfaction of kinematic constraints at the velocity and acceleration levels, thus improving the algorithm stability. Assuming that $\dot {\mathbf {q}}^{*}$ and $\ddot {\mathbf {q}}^{*}$ are the velocities and accelerations obtained after the Newton–Rhapson iteration has converged, these projections can be written for velocities as

$$\left[\frac{d\mathbf{g}(\mathbf{q},\dot{\mathbf{q}})}{d\mathbf{q}}\right]\dot{\mathbf{q}}=[\mathbf{M}+\frac{\mathrm{\Delta }t}{2}\mathbf{C}+\frac{\mathrm{\Delta }{t}^2}{4}\mathbf{K}]{\dot{\mathbf{q}}}^{\ast }-\frac{\mathrm{\Delta }t}{4}{\mathbf{\Phi }}_{\mathbf{q}}^{\mathrm{T}}\alpha {\mathbf{\Phi }}_t$$

(6)

and for accelerations, as

$$\left[\frac{d\mathbf{g}(\mathbf{q},\dot{\mathbf{q}})}{d\mathbf{q}}\right]\ddot{\mathbf{q}}=[\mathbf{M}+\frac{\mathrm{\Delta }t}{2}\mathbf{C}+\frac{\mathrm{\Delta }{t}^2}{4}\mathbf{K}]{\ddot{\mathbf{q}}}^{\ast }-\frac{\mathrm{\Delta }t}{4}{\mathbf{\Phi }}_{\mathbf{q}}^{\mathrm{T}}\alpha ({\dot{\mathbf{\Phi }}}_{\mathbf{q}}\dot{\mathbf{q}}+{\dot{\mathbf{\Phi }}}_t),$$

(7)

where $\mathbf{C}$ and $\mathbf{K}$ are the damping and stiffness contributions of the system. These steps represent the core of the index-3 augmented Lagrangian algorithm with projections (ALI3-P), described in [35] and extended to handle nonholonomic constraints in [36].

2.2 Hydraulic dynamics

The lumped fluid theory [37] is used in this study to model the dynamics of hydraulic actuators. According to this approach, the hydraulic circuit is divided into volumes $V$, in which the pressure $p$ is assumed to be equally distributed. Each volume has sets of incoming and outgoing volume flows, $Q_{\mathrm{in}}$ and $Q_{\mathrm{out}}$, respectively, and an effective bulk modulus $B_{e}$, which describes the flexibility of the oil and the containers, such as hoses and cylinders. The concept of the lumped fluid theory is illustrated in Fig. 1.

Fig. 1

A volume with equally distributed pressure and effective bulk modulus

The pressure $p$ can be solved from the first-order differential equation

$$ \dot{p} = \frac{B_{e}}{V} (Q_{\mathrm{in}} - Q_{\mathrm{out}} ) \text{,} $$

(8)

where $V$ is the volume under investigation. The effective bulk modulus, in turn, takes the form

$$ B_{e} = \Biggl(B_{o}^{-1} + \sum_{l=1}^{n_{c}} \frac{V_{l}}{V B_{l}} \Biggr)^{-1} \text{,} $$

(9)

where $B_{o}$ is the oil bulk modulus, $n_{c}$ is the number of chambers in volume $V$, and $V_{l}$ and $B_{l}$ are the volume and bulk moduli of each chamber, respectively. The volume flows $Q_{\mathrm{in}}$ and $Q_{\mathrm{out}}$, in turn, can originate from, for instance, cylinder movement or pressure difference over a valve, and therefore their equations are problem-dependent. Thus, these equations are presented with the case-example in Sect. 3.

2.3 Co-simulation setup

Figure 2 depicts a generic co-simulation scheme with, for simplicity, two subsystems and a generic co-simulation manager that couples the subsystems and orchestrates the co-simulation process. Each subsystem has its own internal set of states, which is integrated with method $\int _{1}$ and step-size $h_{1}$ in Subsystem 1 and integrator $\int _{2}$ and step-size $h_{2}$ in Subsystem 2. To couple their dynamics, the subsystems exchange information at communication steps of size $H$. Although, in general, all three step-sizes may be different, it is often practical to use the largest internal step-size as the communication step-size. It should also be noted that hereafter the notation of the subscripts $h_{1}$ and $h_{2}$ is for clarity particularized for each example.

Fig. 2

A generic co-simulation scheme

In the used noniterative Jacobi scheme co-simulation, at each communication step the subsystems send outputs $\mathbf{y}_{1}$ and $\mathbf{y}_{2}$ and receive inputs $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$, respectively. During the communication step, both subsystems integrate their states independently. Since the updated values of inputs $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ are not available during this time, data extrapolation may be needed. For the purposes of this study, we implement a polynomial extrapolation and a reduced interface model [24]. The former method uses Lagrange polynomials, that is, it is oblivious to the actual system dynamics, whereas the latter makes a physics-based assumption on the mechanical system dynamics for the update of subsystem inputs, as will be explained next.

2.3.1 Co-simulation via interface models

In co-simulation setups, multibody systems are often coupled to other subsystems with faster dynamics, which need to be integrated with shorter step-sizes. When a multibody system ℳ integrated with step-size $h_{\mathcal{M}}$ is coupled in a multirate Jacobi co-simulation setup to another subsystem $\mathcal{S}$ with faster dynamics, integrated with step-size $h_{\mathcal{S}} < h_{\mathcal{M}}$, the fast subsystem needs to anticipate its inputs between two consecutive communication points. This can be done by means of polynomial extrapolation. It is also possible to build an interface model (IM), to provide a prediction of the evolution of the dynamics of ℳ between communication points [24]. Such an interface model expresses the dynamics of the multibody system ℳ in terms of the variables that describe its interface to $\mathcal{S}$. Let us assume that the dynamics at the interface can be described with a set of $p$ generalized velocities $\mathbf{w}_{\mathrm{i}}$ such that

$$ {\mathbf{w}}_{\mathrm{i}} = \mathbf {A}_{\mathrm {i}} \dot {\mathbf {q}}, $$

(10)

where $\mathbf {A}_{\mathrm {i}}$ is a $p \times n$ Jacobian matrix that makes it possible to express the interface velocities $\mathbf{w}_{\mathrm{i}}$ in terms of the generalized velocities $\dot {\mathbf {q}}$ originally used to describe the multibody system. The matrix $\mathbf {A}_{\mathrm {i}}$ also enables writing the interaction forces $\mathbf {Q}_{\mathrm{i}}$ that couple the dynamics of ℳ and $\mathcal{S}$ in terms of the local coupling forces at the interface ${\mathit{\lambda }}_{\mathrm{i}}$ as

$${\mathbf{Q}}_{\mathrm{i}}={\mathbf{A}}_{\mathrm{i}}^{\mathrm{T}}{\mathit{\lambda }}_{\mathrm{i}}.$$

(11)

Then the multibody dynamics can be expressed in terms of the velocities $\mathbf{w}_{\mathrm{i}}$ that describe the interface as

$${\mathbf{M}}_{\mathrm{eff}}{\dot{\mathbf{w}}}_{\mathrm{i}}={\mathbf{Q}}_{\mathrm{eff}}+{\mathit{\lambda }}_{\mathrm{i}},$$

(12)

where

$$ \mathbf {M}_{\mathrm {eff}} = \left ( \mathbf {A}_{\mathrm {i}} \left (\mathbf {I}- \mathbf {P}_{c} \right ) \mathbf {M}{^{-1}}{\mathbf {A}}_{\mathrm {i}}^{\mathrm {T}}\right ) {^{-1}}$$

(13)

and

$$ \mathbf {Q}_{\mathrm {eff}}= \mathbf {M}_{\mathrm {eff}} \left (\mathbf {A}_{\mathrm {i}} \left (\mathbf {I}- \mathbf {P}_{c} \right ) \mathbf {M}{^{-1}}\left (\mathbf {Q}-\mathbf {c}\right ) + \dot {\mathbf {A}}_{\mathrm {i}} {\mathbf{w}}_{\mathrm{i}} + \mathbf {A}_{\mathrm {i}} \mathbf {P}_{c} \dot{\mathbf{w}}_{\mathrm{i}}\right ), $$

(14)

and where $\mathbf {P}_{c}$ is the projector matrix

$${\mathbf{P}}_c={\mathbf{M}}^{-1}{\mathbf{\Phi }}_{\mathbf{q}}^{\mathrm{T}}{\left({\mathbf{\Phi }}_{\mathbf{q}}{\mathbf{M}}^{-1}{\mathbf{\Phi }}_{\mathbf{q}}^{\mathrm{T}}\right)}^{-1}{\mathbf{\Phi }}_{\mathbf{q}}.$$

(15)

It should be noted that the present form of Eq. (15) relies on the mass matrix $\mathbf {M}$ being invertible, and therefore the method is limited to mechanical systems with a regular mass matrix.

The multibody dynamics provided by Eq. (12) can be used to obtain a physics-based prediction of the evolution of the inputs of the fast subsystem $\mathcal{S}$ between communication points in the way shown in Fig. 3.

Fig. 3

Scheme of the operation of an interface model (IM)

The terms $\mathbf {M}_{\mathrm {eff}}$ and $\mathbf {Q}_{\mathrm {eff}}$ are evaluated at each communication point and sent to the interface model, which keeps them constant during the macrostep. These are used to formulate the dynamics of the interface model, which is integrated in a computationally inexpensive way with the same integration step-size of $\mathcal{S}$. This integration process is updated with the exchange of inputs between the IM and the faster subsystem $\mathcal{S}$. This way the IM provides $\mathcal{S}$ with a prediction of its inputs between communication points, which is based on the dynamics of the multibody system ℳ.

2.3.2 Systematic evaluation of co-simulation accuracy

In the benchmark process that will follow, a systematic method to evaluate the accuracy of the co-simulation is needed. Although the total energy of the system provides useful information about the system-level solution, it gives indication if the energy error is due to the decoupling of the problem or due to the internals of each subsystem. On the other hand, indicators, such as the RMS-value of position errors with respect to a certain reference solution, may be a poor indicator in some cases. For instance, a small phase shift in a high-frequency signal can appear as a large error, even if the resulting solution can be considered otherwise acceptable. For these reasons, this work uses a method based on the work presented in [15], where the energy flows in the co-simulation interface are monitored via the concept of power bonds. This allows the error estimation at the co-simulation interface based on the input and output signals only, namely in terms of energy residuals. In contrast to [15], but similar to [16], the error indicator here can be considered a “generalized energy” and may not strictly represent energy in the physical sense.

To briefly introduce the key concepts related to the energy flow based error estimation, consider a Jacobi scheme co-simulation of two subsystems with scalar inputs and outputs between communication points $t_{n}$ and $t_{n+1}$, as presented in Fig. 4. In the figure, at $t_{n}$, each subsystem sends and receives its input $u_{k} (t_{n})$ and output $y_{k} (t_{n})$, respectively, $k$ being an index of power bond between the systems. After the time integration within each subsystem has taken place without additional communications, at time instant $t_{n+1}$, both subsystems send and receive $u_{k} (t_{n+1})$ and $y_{k} (t_{n+1})$, which are computed based on the previous inputs. With this information only, it is possible, as shown in [15], to estimate the introduced energy error in the co-simulation interface. It is important to note that, from hereafter, the generalized power and energy may not be power and energy in the physical sense if the product of inputs and outputs at each power bond does not correspond to an energy magnitude.

Fig. 4

An explicit Jacobi scheme co-simulation

The concept of power bonds can now be applied to monitor the energy flow between the subsystems [15]. Applied to Fig. 4, this concept states that the energy transfer between Subsystems 1 and 2 via power ports $k_{1}$ and $k_{2}$ is called a power bond $k$. At communication time $t$, Subsystem 1 reports the total transmitted generalized power $P_{k_{1}}$ as

$$ P_{k_{1}}(t) = \tilde{u}_{k_{1}}(t) y_{k_{1}}(t)\text{,} $$

(16)

where $\tilde{u}_{k_{1}}(t)$ is the extrapolated value of the previous input at the current communication time $t$, and $y_{k_{1}}(t)$ is the computed output. Subsystem 2, in turn, reports that the power it received from Subsystem 1 via the power port $k_{2}$ is

$$ P_{k_{2}}(t) = \tilde{u}_{k_{2}}(t) y_{k_{2}}(t)\text{,} $$

(17)

where $\tilde{u}_{k_{2}}(t)$ and $y_{k_{2}}(t)$ are the respective inputs and outputs.

Should the extrapolation within the subsystem be error-free, no energy would be lost or gained at the interface, and balance $-(P_{k_{1}} + P_{k_{2}}) = 0$ would hold. However, as the extrapolation error cannot be, in general, avoided, this balance will be violated, and therefore

$$ -(P_{k_{1}} + P_{k_{2}}) \neq 0 . $$

(18)

The above violation in the generalized power can now be exploited to monitor the co-simulation accuracy:

$$ \delta P_{k} = -(P_{k_{1}} + P_{k_{2}}), $$

(19)

where $\delta P_{k}$ is a generalized power residual of the power bond $k$, that is, the spuriously introduced power due to extrapolation errors. For multiple power bonds, this figure can be computed as

$$ \delta P = -\tilde{\mathbf{u}}{^{\mathrm {T}}}\mathbf{y}, $$

(20)

where $\tilde{\mathbf{u}}$ and $\mathbf{y}$ are the vectors of extrapolated inputs and computed outputs, respectively. Should the number of inputs and outputs not be equal, as it will be the case with the used examples, only a subset of inputs or outputs is used to define a power bond. Following the notation introduced in Fig. 4, now we can write Eq. (20) as

$$ \delta P_{k}(t_{n+1}) = -\tilde{ \mathbf{u}}_{k}(t_{n+1}){^{\mathrm {T}}}\mathbf{y}_{k} (t_{n+1}), $$

(21)

which can be implemented in the co-simulation manager in a straightforward manner. The generalized residual energy is, in turn, defined piecewise between each communication instance with integral

$$ \delta E_{k} (t_{n+1}) = \int _{t_{n}}^{t_{n+1}} \delta P_{k} (t)\, dt, $$

(22)

which now exactly defines the incorrectly added or dissipated energy that results from the extrapolation errors during the communication step $t_{n} \rightarrow t_{n+1}$. A numerical approximation can be used to compute Eq. (22) in the practical implementations. For $m = 0$, a rectangular rule can be used, for example,

$$ \delta E_{k} (t_{n+1}) \approx \delta P_{k}(t_{n+1}) H (t_{n+1}) \text{,} $$

(23)

where $H (t_{n+1})$ is the communication step size. For $m > 0$, higher-order quadrature rule should be used. Trapezoidal rule, which reads as

$$ \delta E_{k} (t_{n+1}) \approx -\frac{H}{2} \big(\delta P_{k}(t_{n}) + \delta P_{k} (t_{n+1})\big), $$

(24)

is a straightforward solution with linear extrapolation. Although higher-order formulas, such as Simpson’s rule, could provide a more accurate approximation, the outputs $\mathbf{y}_{k}$ are not known between communication instants, and therefore in this work, for extrapolation orders $m > 1$, we use the trapezoidal rule.

It should be noted that the residual energy in Eq. (22) defines the spurious energy added to the system or removed from it due to the decoupling of the original problem, regardless of the internal energy conservation in the subsystems. This means that the power sources or sinks in the subsystems are not visible in the residual, and therefore in this paper, two definitions exist for simulation accuracy. Firstly, there is the coupling accuracy defined in terms of Eq. (22) and its approximations. Secondly, there is system-level accuracy, which is defined in this paper by the energy balance of the entire system. While the former has an effect on the latter, since the decoupling error is added to the simulation error, the latter may not be visible in the former. These two definitions are jointly used to provide a comprehensive view to the accuracy of the simulation.

3 Examples

We implemented two nonlinear problems to evaluate co-simulation performance under different coupling approaches. Firstly, we used a planar model of a crane with a single actuator to examine the co-simulation with two interacting subsystems only. Secondly, we obtained a system with multiple actuators by attaching a second hydraulic cylinder to the crane. We modeled and integrated the dynamics of the two hydraulic actuators separately. In practice, hydraulic actuators are coupled through the dynamics of the mechanical system. In a Jacobi co-simulation setup, the integration of these two hydraulic subsystems proceeds separately between communication points, without information about the evolution of the other actuator. For this reason, keeping the co-simulation of the two-actuator model stable is more challenging than in the case of its single-actuator counterpart.

3.1 Multibody models

3.1.1 Single-actuator model

The first example is a two-link hydraulically actuated crane, shown in Fig. 5. A similar model was used as a test problem in [24] and [25].

Fig. 5

A single-actuator crane

The system has two degrees of freedom, one of which is controlled by a hydraulic actuator between the center point $\mathsf {P}$ of the first link and point $\mathsf {B}$ on the ground. The first link in the system is a rigid rod of length $L$ with a homogeneously distributed mass $m$. Link 2 is a massless rod of length $L_{\mathrm{h}}$. Point masses $m_{\mathrm{p}}$ and $m_{\mathrm{h}}$ are located at the tips of the first and second links, points $\mathsf {Q}$ and $\mathsf {R}$, respectively. The set of generalized coordinates chosen to describe the example motion is $\mathbf {q}= \left [\begin{array}{ ccccc } x_{\mathsf {P} } & y_{\mathsf {P} } &\theta _{1} &x_{\mathsf {R} } & y_{\mathsf {R} } \end{array}\right ] {^{\mathrm {T}}}$; three kinematic constraints $\mathbf{\Phi }=\mathbf{0}$ are used to enforce that the distance between points $\mathsf {Q}$ and $\mathsf {R}$ remains constant and that angle $\theta _{1}$ is proportional to the $x$ and $y$ coordinates of point $\mathsf {P}$. The interface velocity $w_{\mathrm{i}}$ in this example corresponds to the actuator rate $\dot{s}_{1}$ and can be related to the generalized velocities as follows:

$$ w_{\mathrm{i}} = \dot{s}_{1} = \left [ \textstyle\begin{array}{ ccccc } \dfrac{x_{\mathsf {P} }-x_{\mathsf {B} }}{s_{1}} & \dfrac{y_{\mathsf {P} }-y_{\mathsf {B} }}{s_{1}} &0 &0 &0 \end{array}\displaystyle \right ] \dot {\mathbf {q}}= \mathbf{A}_{\mathrm{i}} \dot {\mathbf {q}}. $$

(25)

The force transmitted at the interface $\lambda _{\mathrm{i}}$ corresponds to the force exerted by the cylinder.

The physical properties of the single-actuator crane are summarized in Table 1.

Table 1 Mechanical parameters of the single-actuated model

3.1.2 Double-actuator model

The second example is a crane with two hydraulic actuators obtained after adding a second cylinder, identical to the previously existing one, between point $\mathsf {P}$ and the center point of link 2, as shown in Fig. 6. Such a system was used in [24] to assess the behavior of interface models in the co-simulation of multiply actuated systems.

Fig. 6

A double actuator crane

The set of variables used to describe the system motion remains unchanged. However, a new component of the interface velocity is required to describe the interaction between the second actuator and the mechanical system. The rate $\dot{s}_{2}$ of the second cylinder is selected to this end:

$$\begin{aligned} \mathbf{w}_{\mathrm{i}} = \left[\textstyle\begin{array}{c} \dot{s}_{1} \\ \dot{s}_{2}\end{array}\displaystyle \right] = \left[\textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} \dfrac{x_{\mathsf {P} }-x_{\mathsf {B} }}{s_{1}} & \dfrac{y_{\mathsf {P} }-y_{\mathsf {B} }}{s_{1}} &0 &0 &0 \\ 0 &0 &0 &\dfrac{x_{\mathsf {R} }}{4~s_{2}} &\dfrac{y_{\mathsf {R} }}{4~s_{2}}\end{array}\displaystyle \right] \dot {\mathbf {q}}. \end{aligned}$$

(26)

The mechanical properties of the system remain the same as those in the single-actuator model.

3.2 Hydraulic circuit

Figure 7 depicts the hydraulic actuator circuit used in both examples in Sect. 3. In the latter, this circuit is duplicated with slightly different parameters to add the second actuator to the system. The circuit includes a pump and a tank, which are assumed to have constant pressures $p_{P}$ and $p_{T}$, respectively. Oil volumes flow from or to these sources via a $4/3$ directional valve, which is controlled with a control signal $U$ produced with a joystick, the positive direction of the signal being depicted in the figure. The circuit is split into three volumes, one between the directional and a throttle valve and one at each side of the piston head. Additionally, the figure depicts the total cylinder length $l$ and lengths at rod and piston sides, $l_{2}$ and $l_{1}$, respectively. A full list of the used parameters is presented in Table 2.

Fig. 7

A hydraulic circuit

Table 2 Parameters of the example circuit

The hydraulic circuit has three volumes, and therefore applying Eq. (8) yields three first-order differential equations, which describe the system dynamics:

$$\begin{aligned} \dot{p}_{1} = \frac{B_{e1}}{V_{1}}(Q_{31} - A_{1} \dot{s}), \end{aligned}$$

(27)

$$\begin{aligned} \dot{p}_{2} = \frac{B_{e2}}{V_{2}}(A_{2} \dot{s} - Q_{2V}), \end{aligned}$$

(28)

$$\begin{aligned} \dot{p}_{3} = \frac{B_{e3}}{V_{3}}(Q_{V3} - Q_{31}), \end{aligned}$$

(29)

where $A_{1}$ and $A_{2}$ are piston and rod side areas in the cylinder, and $\dot{s}$ is the actuator rate. The volume flows over the valves, $Q_{2V}$, $Q_{V3}$, and $Q_{31}$, are computed with the semiempirical model [38], and their positive directions are depicted in Fig. 7. The expression for the volume flow over the throttle valve can be written as

$$ Q_{31} = \textstyle\begin{cases} C_{t} \sqrt{ \vert \Delta p \vert }, & \Delta p > 0, \\ 0 , & \Delta p = 0, \\ -C_{t} \sqrt{ \vert \Delta p \vert }, & \Delta p < 0, \\ \end{cases} $$

(30)

where $C_{t}$ is the semiempirical coefficient, and $\Delta p$ is the pressure difference over the valve. In the case of a simple throttle valve the semiempirical coefficient is computed as $C_{t} = C_{d} A_{t} \sqrt{\frac{2}{\rho }}$, where $C_{d}$ is the discharge coefficient, $A_{t}$ is the throttle area, and $\rho$ is the oil density. The volume flows $Q_{2V}$ and $Q_{V3}$ over the directional valve take the form

$$ Q = \textstyle\begin{cases} C_{v} U \sqrt{ \vert \Delta p \vert }, & \Delta p > 0, \\ 0 , & \Delta p = 0, \\ -C_{v} U \sqrt{ \vert \Delta p \vert }, & \Delta p < 0, \end{cases} $$

(31)

where $C_{v}$ is the semiempirical coefficient computed based on one operational point of the valve, and $U$ is the control signal computed as

$$ \dot{U} = \frac{U_{\mathrm{ref}} - U}{\tau }, $$

(32)

where $U_{\mathrm{ref}}$ the reference, in this case from the joystick, $U$ is the achieved position, and $\tau = \left (2\pi f_{-45}\right )^{-1}$ is the time constant computed based on the −45°phase shift frequency.

The force produced by the cylinder can be computed as follows:

$$ F_{s} = p_{1}A_{1} - p_{2} A_{2} - c \dot{s} - F_{d} \text{,} $$

(33)

where $c$ is a viscous coefficient of $1\times 10^{5}$ Ns/m, and $F_{d}$ is the end damper force computed in this case as

$$ F_{d} = \textstyle\begin{cases} k_{d}(l_{d} - l_{1}) - c_{d} \dot{s}, & l_{1} \leq l_{d}\text{,} \\ -k_{d}(l_{d} - l_{2}) - c_{d} \dot{s}, & l_{2} \leq l_{d}\text{,} \\ 0 & \text{otherwise,} \end{cases} $$

(34)

where $k_{d}$ and $c_{d}$ are, respectively, the spring and damping coefficients of the dampers, and $l_{d}$ is the damper length. Note that the cylinder lengths $l_{1}$ and $l_{2}$, as depicted in Fig. 7, include the dampers.

3.3 Co-simulation setup

The described systems are coupled by the co-simulation manager via a noniterative Jacobi scheme. The mechanical system sends at each communication point the actuator length $s$, velocity $\dot{s}$, and acceleration $\ddot{s}$ to the manager. The hydraulic systems, in turn, send either the cylinder force $F_{s}$ or the pressures $\mathbf{p} =[p_{1}\ p_{2} \ p_{3}]$ to the interface. In the latter option the cylinder force is computed within the multibody system; accordingly, information about the properties of the hydraulic system needs to be available to the multibody subsystem, so that it can evaluate Eq. (33) to obtain the cylinder force. These two options are hereafter referred to as velocity–force and velocity–pressure couplings, and both can be used with and without the reduced interface model. Polynomial extrapolation, based on the Lagrange polynomials, is used on the inputs of the hydraulic system when the IM is not used.

Regarding time integration, the mechanical system is solved with the index-1 and index-3 augmented Lagrangian algorithms described in Sect. 2.1 and shares the constant communication step-size of the co-simulation manager, that is, $h_{m} = H$. The ALI1 formulation is integrated with a semi-implicit forward Euler formula, and the ALI3-P with the trapezoidal rule integrator. The hydraulic system, in turn, has a smaller internal step-size, that is, $h_{h} < H$, which forms a matching grid with the communication instances. To evaluate the effect of the internal integration scheme of the hydraulic subsystem, three different integrators are implemented for experimenting, namely the forward Euler formula (labeled Euler) as a single-step explicit integrator, the fourth-order Runge–Kutta method (RK4) as a multistep explicit integrator, and the trapezoidal rule (TR) as an implicit integrator. Fixed step-sizes are used during the integration process. These configurations are summarized in Table 3.

Table 3 Configuration variables

4 Results

The given examples are used in the benchmark to gain insight into the noniterative Jacobi-scheme co-simulation of multibody and hydraulic dynamics. The focus of the benchmark is on the influence of the configuration on the achievable communication step-size, since it is known to have a significant effect on co-simulation accuracy and efficiency [12, 18]. To this end, the internal step-size $h_{h}$ of hydraulics is locked to a value, in this case 0.1 ms, small enough to obtain a stable solution, whereas $H$ is varied under each tested co-simulation configuration. The generalized residual energy at the end of the simulation is taken as an indicator of the coupling accuracy, whereas the overall solution accuracy is measured by the energy balance of the entire system. In addition to the study of the achievable communication step-size with a predetermined value of $h_{h}$, here the effect of the internal integrator of the hydraulics subsystem on the achievable internal step-size is also briefly touched upon. A monolithic solution, which follows the integration scheme shown in [26], is used to produce a set of reference results with the ALI3-P formulation. As the case examples are implemented in a Matlab environment and the used implementation does not produce comparable computation times, here CPU-times are omitted.

4.1 Single actuator model

For the first studied case, that is, the single actuator crane in Sect. 3.1.2, a reference signal $U$ is produced for the directional valve with a joystick. The resultant actuator length that corresponds to this input is depicted in Fig. 8. During the work cycle, the actuator hits both end dampers, the models of which are similar to contact models, which can be assumed to be challenging for the integrators. Figure 8 shows the reference result obtained with a monolithic simulation at 1 ms step size. This reference solution is compared to others obtained with Jacobi-scheme co-simulation using different coupling strategies.

Fig. 8

Simulated work cycle

The effect of the communication step-size $H$, used integrators $\int _{m}$ and $\int _{h}$, and extrapolation methods, that is, $m = 2$ and IM, on the coupling accuracy are shown in Figs. 9a–9f, where the generalized energy has been normalized so that the configuration of ALI1 with second-order extrapolation and forward Euler integrator for hydraulics at $h_{h} = 1~\mbox{ms}$ internal step size equals zero. It must be noted that the generalized energy residuals are directly comparable when the coupling variables are the same, that is, within left- and right-hand side columns of figures, and in both cases, the lower energy residual indicates a better solution. The achievable step-size is taken to be the segment of the curve where the accumulated energy residual is relatively low, almost linearly dependent on the communication step size $H$. When the curve “takes off”, accuracy of the coupling and therefore the simulation can be assumed to be lost.

Fig. 9

Effect of the communication step-size $H$ on coupling accuracy with different hydraulic integrators and coupling variables; single-actuator model

As Fig. 9 shows, the index-1 formulation can achieve significantly larger communication step-sizes with a similar level of coupling error introduced at the co-simulation interface, when compared to the index-3 formulation. The lower energy residual likely results from the explicit integrator used with the ALI1 formulation, which eliminates the need for signal extrapolation, and thus is indifferent to the co-simulation scheme. However, as Fig. 10 shows, in turn, the index-1 solution introduces losses in the total system energy as the simulation progresses, whereas the index-3 solution maintains system-level accuracy until the stability of the solution is lost. This difference may stem from the dissipative terms in the ALI1 formulation and, in addition, from the used explicit integrator. This demonstrates the fact that the coupling accuracy alone cannot be used to fully determine if the solution for the system is accurate, even if it can be indicative of it, as is the case with the ALI1 solutions. However, on the other hand, the coupling accuracy can be computed even if the system level energy balance is not known and can be used, as Fig. 9 shows, as an indirect indicator of solution stability.

Fig. 10

Energy balance with Euler integrator at $h_{\text{hyd}} = 0.1\text{ ms}$, polynomial extrapolation, and velocity-force coupling; single-actuator model

Regarding the use of the different coupling variables, the use of velocity pressure coupling seems to yield more accurate coupling and stable results, as Fig. 9 shows. This is true for both implemented formulations for the mechanical system, although in this case the index-1 formulation seems to benefit more from the use of pressures in the coupling, especially when the reduced interface model is used. This demonstrated that difference originates from the hydraulic force description in Eq. (33), i.e., $F_{s} = p_{1} A_{1} - p_{2} A_{2} - c \dot{s} - F_{d}$, which couples the mechanical and hydraulic states. If velocity-force coupling is used, then this equation is evaluated within the hydraulic subsystem, and for this reason, the actuator rate $\dot{s}$ is that extrapolated from the previous time step. This causes a direct feed-through, that is, the output of the hydraulic subsystem directly depends on its inputs. It has been shown in [8, 39] that algebraic loops caused by direct feed-through are detrimental for the stability of the cosimulated solution. Conversely, if a velocity–pressure coupling is used, then Eq. (33) is evaluated within the multibody system with the integrated value of $\dot{s}$ and pressures $\mathbf{p}$ obtained from the previous time step, which removes the direct feed-through and therefore, as demonstrated here, yields a more stable simulation.

As demonstrated in the literature, the used extrapolation method has a significant effect on co-simulation accuracy [12]. For the purposes of this study, polynomial extrapolation is used to compare to the interface model [24], which allows us to extrapolate the data from the mechanical system based on the system dynamics. The order of polynomial extrapolation was decided based on experiments with a configuration of ALI1, Euler, and $h_{h} = 0.1~\mbox{ms}$. As the results show, with the single-actuator model, polynomial extrapolation is performed in most cases better than the reduced interface model, both in the sense of solution and coupling accuracy and achievable communication step-size. Figure 11 exemplifies the difference in results between the polynomial extrapolation and reduced interface model. As we can see, use of the reduced interface model may lead to degraded velocity level solution, whereas smaller deviations can be seen at the configuration level. This behavior also explains the increase in the residual energies seen in Fig. 9a.

Fig. 11

Polynomial extrapolation vs. IM; single-actuator model

Figures 12 and 13 summarize the effect of the internal configuration of the fast subsystem on co-simulation accuracy. As already demonstrated in Fig. 9, the internal integrator of the fast subsystem, hydraulics in this case, in turn, seems to have only little impact on the achievable communication step-size $H$, should the internal step $h_{h}$ be sufficiently small to ensure a stable integration. For a clearer demonstration, in Figs. 12a–12b the normalized energy residuals with different integrators for the hydraulic subsystem are presented when the rest of the system configuration is kept constant. This suggests that changing the hydraulic integrator may not allow larger communication step-sizes to be taken, should the internal step-size be small enough for the used integrator to produce an accurate solution for the system. Regarding the possibility to increase the internal step-size $h_{h}$, in turn, as demonstrated in Figs. 13a–13b, it seems to be possible only when $H$ is sufficiently small. Notice how the solution with the trapezoidal rule integrator at $H = 4~\mbox{ms}$ becomes unstable, whereas with $H = 1~\mbox{ms}$ a stable solution can be found. The solution accuracy behaves in a similar manner. At $H = 1~\mbox{ms}$ the solutions with $h_{h} = 0.1~\mbox{ms}$ and $h_{h} =0.5~\mbox{ms}$ are closer to the monolithic solution than with $H = 4~\mbox{ms}$. This shows that in the studied case, there is a coupling between achievable $H$ and the used $h_{h}$. This further indicates that depending on the location of the computational bottleneck of the co-simulation setup, it may be possible to optimize the overall efficiency by tuning step-sizes in fast and slow subsystems together.

Fig. 12

Effect of the hydraulics integrator on the coupling accuracy; single-actuator model

Fig. 13

Effect of the step size to the solution accuracy; single-actuator model

4.2 Double actuator model

Similarly to the first case, for the double-actuator model, reference signals for the directional valves of the hydraulic circuits are produced with a joystick. Figure 14 shows the resultant actuator lengths obtained with a monolithic simulation of 1 ms step-size. Although the system still seems fairly simple, the addition of the second hydraulic circuit as a new subsystem introduces an indirect coupling between the cylinders via the mechanical system, which from the co-simulation point of view, is significantly more challenging than the single-actuator case. The problem is solved again with the Jacobi-scheme co-simulation with different coupling strategies, as described before. To limit the number of considered variables, the same settings are used for both hydraulic subsystems.

Fig. 14

Simulated work cycle with the double-actuator model

Figure 15 shows the resultant generalized energy residuals, normalized in a similar manner to the previous case, that is, the ALI1 formulation with forward Euler integrator for the hydraulics, and first-order extrapolation at $H = 1~\mbox{ms}$ is set to one in the chart. As done with the single-actuator model, the internal step-size for both hydraulic subsystems has been kept at a constant of $h_{h} = 0.1~\mbox{ms}$ in all cases, whereas other configuration parameters are varied.

Fig. 15

Effect of the communication step-size $H$ on the coupling accuracy with different hydraulic integrators and coupling variables; double-actuator model

As can be seen from Fig. 15, the achieved communication step-size $H$ is, in general, smaller than those shown in Fig. 9 for the single-actuator model. This result can be expected, as the previously explained indirect coupling of the hydraulic subsystems makes this case more challenging for the used noniterative Jacobi scheme. An exception to the rule is the reduced interface model combined with the ALI1 formulation, which seems to be able to achieve communication step-sizes close to those of the single-actuator model without loss of stability. This result, in turn, agrees to the original proposal of the method [24] and, in addition, shows that the method seems to be viable with stiffer hydraulic models. The explanation behind the large communication step-sizes achieved with IM lies in its capability to maintain the coupling of the hydraulic subsystems between the communication instants.

Otherwise, the results follow the trends set by the single-actuator case. A clear difference can be seen in the results between the velocity–force and velocity–pressure couplings, the latter being able to achieve larger communication step-sizes. In this case, however, the benefit of the velocity–pressure coupling is more clear than in the single-actuator case, and it suggests that in more complex cases the elimination of the direct feed-through may be advisable in the used co-simulation scheme. The effect of the internal solver of the fast subsystem in this case, as Fig. 15 shows, also has little effect on the achievable communication step-size. This is again demonstrated in Figs. 16 and 17, where the generalized energy residuals are compared between the different subsystem integrators. The relation between the achievable communication step-size $H$ and the internal step-size $h_{h}$ is also exemplified. As we can see, should the communication step-size be small enough, the trapezoidal rule integrator can still achieve 1 ms step-sizes without losing stability, although with the increase in the communication step-size this property is quickly lost.

Fig. 16

Effect of the hydraulics integrator on the coupling accuracy; double-actuator model

Fig. 17

Effect of the step size to the solution accuracy; double-actuator model

5 Conclusions

In this work, we presented a benchmark study to gain insight into noniterative Jacobi-scheme co-simulation of multibody and hydraulic dynamics. Two nonlinear examples are used for this purpose, with different approaches to the integration procedures and to the coupling. Their effects on the coupling and system accuracy have been investigated, and the key conclusions can be summarized as follows:

• The selection of variables for the information exchange can have a significant effect on the achievable communication step-size, since they may introduce or remove direct feed-through in the coupling scheme. In the light of the presented results, it seems that the use of the state variables as subsystem outputs in the coupling should be recommended.

• The used multibody algorithm, and especially the used integrator, can have a significant effect on the achievable communication step-size and also on the accuracy of both the solution and the coupling itself. However, an accurate system level solution may not be the most stable one and vice versa, as the results show.

• The interface model (IM) can allow a significant increase in the communication step size, but with a cost of degraded velocity level solution in some cases. The acceptability of this trade-off depends on the application, and therefore no definitive guidelines can be given regarding its use.

• If the step-size of the hydraulics integrator is small enough, here 0.1 ms, then the selection of the internal integrator seems to have no practical effect on the achievable communication step-size. In the examples reported in this research the multibody algorithm, the choice between polynomial extrapolation or use of an IM, and the selection of the coupling variables had a more noticeable effect on the achievable $H$ when $h_{h}$ was kept as the above-mentioned constant.

• Should the internal step-size of the hydraulic subsystem be increased, in turn, the results suggest that a relation exists between the largest achievable communication step-size and the internal step-size used in the fast subsystem. Thus, depending on the location of the computational bottleneck of the co-simulation setup, it may be possible to optimize the overall efficiency by tuning step-sizes in fast and slow subsystems together.