Comparing double-step and penalty-based semirecursive formulations for hydraulically actuated multibody systems in a monolithic approach

The simulation of mechanical systems often requires modeling of systems of other physical nature, such as hydraulics. In such systems, the numerical stiffness introduced by the hydraulics can become a significant aspect to consider in the modeling, as it can negatively effect to the computational efficiency. The hydraulic system can be described by using the lumped fluid theory. In this approach, a pressure can be integrated from a differential equation in which effective bulk modulus is divided by a volume size. This representation can lead to numerical stiffness as a consequence of which time integration of a hydraulically driven system becomes cumbersome. In this regard, the used multibody formulation plays an important role, as there are many different procedures for the constraint enforcement and different sets of coordinates to choose from. This paper introduces the double-step semirecursive approach and compares it with a penalty-based semirecursive approach in case of coupled multibody and hydraulic dynamics within the monolithic framework. To this end, hydraulically actuated four-bar and quick-return mechanisms are analyzed as case studies. The two approaches are compared in terms of the work cycle, energy balance, constraint violation, and numerical efficiency of the mechanisms. It is concluded that the penalty-based semirecursive approach has a number of advantages compared with the double-step semirecursive approach, which is in accordance with the literature.


Introduction
The computer simulation of dynamic systems has proven its value in product development, from early prototyping phase to user training, and emerging digital twins and artificial intelligence applications. In practice, the modeling of mechanisms can be carried out by using multibody system dynamics in which equations of motion describe a force equilibrium for the mechanical system under investigation. The use of multibody system dynamics also allows describing actuator models such as hydraulics or electric drives.
The multibody methods can be, in general, categorized to two main groups according to the selected coordinates [25]. Firstly, there is the group of formulations based on global coordinates, where the absolute positions, velocities, and accelerations of each body are used. In the second case, relative coordinates are employed. Should high computational efficiency be required, as is the case in real-time applications, such as those presented in [22,23], the use of the relative coordinates is often considered to be an appropriate approach. The selection of an approach is case-dependent, as demonstrated in [17]. The computational efficiency can be significantly effected by implementation details [18], and, among others, the use of the automated differentiation tools [6], and sparse and parallelization techniques [20].
Within the family of methods based on the relative coordinates, the semirecursive approach is often used. In the semirecursive approach, closed loops need to be handled by employing constraint equations. Constraint equations, in turn, can be accounted in the semirecursive approach in many different ways, such as by using the Lagrange multiplier method, the penalty-based approach as proposed by Cuadrado et al. [10] or the double-step approach, which is using a coordinate partitioning as proposed by Rodríguez et al. [33]. The latter two approaches are originated from Featherstone's articulated inertia method [13]. The penaltybased approach utilizes the index-3 augmented Lagrangian formulation with projections [3,9] to enforce the constraints. After their original introductions, both methods have also been successfully used in practical applications, such as in real-time vehicle simulations [20,29]. In this study, the double-step approach based on coordinate partitioning method [24,33] is referred as the double-step semirecursive formulation, and the penalty-based augmented Lagrangian approach [9,10] is referred as the penalty-based semirecursive formulation.
Regarding the solution of the multiphysics problem that system-level simulations require, two main approaches exist in the literature. A straightforward one is the monolithic approach, or sometimes referred as unified scheme, where a single set of differential equations is formed and integrated forward in time as a whole. In an alternative approach, namely cosimulation approach, in turn, the system is split into two or more subsystems that are each integrated separately. In this approach, the required variables, such as the state variables, are communicated in predetermined time intervals. Multiple instances of both approaches can be found from the literature. Cosimulation, due to the discrete time information exchange and the resulting coupling error, has especially been under keen research interest in recent years [16,36]. The studies include important aspects, such as the multirate cosimulation [18], cosimulation configuration [4,30], and energy-based coupling error minimization [4,34,35]. Monolithic schemes, in turn, have been under less active development since the simple coupling requires less research effort. Nevertheless, in [12] a multiphysics model was derived for a semiactive car suspension. Naya et al. [28], in turn, presented a monolithic method for a coupled multibody and hydraulic simulation based on the index-3 augmented Lagrangian [3], followed later in [32] and [31], where the index-3 semirecursive formulation [10] was used.
In the case of a coupled multibody and hydraulic dynamics, which is a common case in mobile working machinery, one of the key aspects to consider is the numerical stiffness introduced by the hydraulics. This problem can be alleviated by proper selection of a multibody approach. While these aspects have been discussed in some sources, such as in [32], few comparisons between the available methods in this context exist in the literature. In the context of the absolute nodal coordinate formulation, in turn, a study in this direction has been done by Matikainen et al. [27], where coordinate partitioning, Lagrange multiplier method with Baumgarte's stabilization, and penalty formulation were used for constraint enforcement. The results indicated best performance for the Lagrange multiplier method, closely followed by the penalty-based approach, while less efficient solution was sought with the coordinate partitioning method.
The objective of this paper is to introduce the double-step semirecursive formulation [24] and compare it with the penalty-based semirecursive formulation [10,32] in the framework of hydraulically driven systems. A monolithic scheme for the coupled simulation of the double-step semirecursive formulation and hydraulic systems is introduced in this study. As explained earlier, the modeling of hydraulic actuators often leads to numerically stiff systems. In this study, the hydraulic system will be described by using the lumped fluid theory. While variable-step integrators often provide more efficient solutions, especially with stiff systems when compared with fixed step-size solutions, since the author's interests lie within the field of real-time simulation, only fixed step-sizes are considered in this work. The coupled systems are referred after the name of the multibody formulations used, such as the double-step semirecursive approach and the penalty-based semirecursive approach. As case studies, hydraulically actuated four-bar and quick-return mechanisms are illustrated. Using the numerical examples, the two approaches are compared based on the work cycle, energy balance, constraint violation, and numerical efficiency of the mechanisms.

Semirecursive multibody formulations
The dynamics of a constrained mechanical system can be described by using a multibody system dynamics approach. In the semirecursive formulations, the dynamics of the openloop multibody systems are formulated in relative joint coordinates, which are independent. In the case of closed-loop multibody systems, in turn, the relative joint coordinates are not independent and cut-joint method is often used to open the loop. In this study, the cut-joint constraints are incorporated by using the coordinate partitioning method [24,33], referred as the double-step semirecursive formulation, and by using the penalty-based augmented Lagrangian method [2,9], referred as the penalty-based semirecursive formulation. Since the system is hydraulically actuated, the internal dynamics of the hydraulics are computed and the resultant force, as well as the stroke and stroke velocity, is used to combine the hydraulics to the multibody equations of motion. Therefore, in this study, the constraints are assumed scleronomic. Should kinematic constraints be needed, both methods can easily be extended to rheonomous systems, as shown in the literature [10,24].

Semirecursive formulation for an open-loop system
In the semirecursive formulation, a system is considered as an open-loop multibody system with N b bodies. The following Cartesian velocity Z j and accelerationŻ j are used to describe the absolute velocity and acceleration of each body [24] whereṡ j ands j are, respectively, the velocity and acceleration of the point of the body j , which at that particular time is coincident with the origin of the inertial reference frame, and ω j andω j are the angular velocity and angular acceleration, respectively, of the body j . In this approach, the kinematics of the open-loop system is calculated in a recursive form, either from the base to the leaves, as performed in this study, or from leaves to base, by applying the classical kinematic relations as in [11]. Figure 1 shows an example of an openloop system. In general case, the position of the system can be described by using the relative joint coordinates z = z 1 , z 2 , . . . , z N b T .
The absolute velocity Z j and accelerationŻ j for body j ∈ [1, N b ] can be recursively expressed in terms of the previous bodies as [24] where the scalarsż j andz j are the first and second time derivatives, respectively, of the relative joint coordinate z j , and the vectors b j and d j depend on the type of joint [11] that connects the bodies j − 1 and j . Note that the indexes j − 1 and j may not be successive, as the system may branch. The absolute velocities Z ∈ R 6N b and accelerationsŻ ∈ R 6N b of the open-loop system can respectively be expressed in the matrix forms as A velocity transformation matrix R ∈ R 6N b ×N b that maps the absolute velocities into a set of relative joint velocities can be written as [11,24] whereż ∈ R N b is the vector of relative joint velocities, T ∈ R 6N b ×6N b is the path matrix that represents the topology of the open-loop system, and R d ∈ R 6N b ×N b is a diagonal matrix whose elements are the vectors b j arranged in an ascending order. The path matrix T is a lower block-triangular matrix [11] whose elements from diagonal to left are (6 × 6) identity matrices I 6 , if the corresponding body is in between the considered body and the root of the system. Note that the termṘż in Eq. (5) can be expressed in terms of the vectors d j by using Eq. (3) [11]. The absolute velocity Y j and accelerationẎ j of the center of mass of body j can be written as where g j is the position vector of the center of mass of body j with respect to the inertial reference frame, I 3 is a (3 × 3) identity matrix, andġ j ,g j , and g j , respectively, are the first and second time derivative, and the skew-symmetric matrix of the position vector g j . By using Eqs. (6) and (7), the virtual power of the inertia and external forces acting on the open-loop system can be written as [24] where the virtual velocities are denoted with an asterisk (*), and the matrices M j ,M j , and Q j can be written as andQ where m j is the mass of body j and J j is the inertia tensor of body j . The inertia tensor J j can be written as J j = A T jJ j A j , whereJ j is the constant inertia tensor in the body reference frame of body j and J j is expressed in the inertial reference frame. Note that both J j and J j are defined with respect to the center of mass of body j . In Eq. (10), f j is the vector of external forces applied on body j and τ j is the vector of external moments with respect to the center of mass of body j . By substituting Eqs. (4) and (5) in Eq. (8), a set of ordinary differential equations that describes the motion of the open-loop system can be written as [24] whereM ∈ R 6N b ×6N b is a diagonal matrix that consists of the mass matricesM j , andQ ∈ R 6N b is a column vector that consists of the force vectorsQ j . Equation (11) can be rewritten where

Double-step semirecursive formulation
In the double-step semirecursive formulation, the constraints are introduced by using the coordinate partitioning method [24]. This method can be used in frameworks of the explicit and implicit integrators. A set of m constraint equations = 0 are used for the closure of an open-loop. For the sake of simplicity, the constraint equations are assumed holonomic and scleronomic. To account for rheonomic constraints, the reader is referred to [24]. The constraint equations = 0 can be expressed in terms of the relative joint coordinates as (z) = 0. Figure 2 shows an example of a closed-loop system. By using the coordinate partitioning method, the dependent velocities can be written in terms of the system's degrees of freedom f as [25] where d z ∈ R m×m and i z ∈ R m×f are, respectively, the dependent and independent columns of Jacobian matrix z , andż d ∈ R m andż i ∈ R f are, respectively, the dependent and independent relative joint velocities. It is assumed that neither redundant constraints nor singular configurations exist, which guarantees that the inverse of d z can be found. A velocity transformation matrix R z ∈ R N b ×f is introduced to transform the dependent relative joint velocities into independent ones as [25] Similarly, accelerations can be written as wherez i ∈ R f are the independent relative joint accelerations. In this study, the independent relative joint coordinates are identified by using the Gaussian elimination with full pivoting to the Jacobian matrix z . Note that Gaussian elimination with full pivoting is a lower-upper (LU) matrix factorization technique, where the rows and columns of a matrix are interchanged to use the largest element (in absolute value) in the matrix as the pivot [21]. Accordingly, the (N b − m) columns of the Jacobian matrix z , where the m pivots do not appear, determines the independent relative joint coordinates [15,25]. Fisette and Vaneghem utilized the same technique in the identification of dependent and independent coordinates [14]. This can be considered as a relative drawback of the double-step semirecursive formulation because the penalty-based semirecursive formulation utilizes the full set of coordinates. By substituting Eqs. (14) and (15) into Eq. (11), the dynamic equations of the closed-loop system can be written as [7] where + TṘ dż are the absolute accelerations, when vectorz in Eq. (5) is set to zero. Equation (16) can be expressed in a simplified form as where Note that the relation between Q and Q can be written as

Penalty-based semirecursive formulation
In the penalty-based semirecursive formulation, the constraints are introduced by using the penalty-based augmented Lagrangian method [10]. In this formulation, the time integration scheme is carried out by using the trapezoidal rule. The loop-closure constraints, a set of m constraint equations = 0, are introduced in Eq. (12) with a penalty method similar to the index-3 augmented Lagrangian with projections to satisfy the constraints on velocity and acceleration levels [3,9]. The equations of motion for the closed-loop system can be written as where z is the Jacobian matrix of (z) = 0, α is the penalty factor that can be set the same for all constraints, and λ is the vector of iterated Lagrange multipliers. In this method, these multipliers are obtained at each time-step k as where h is the iteration step. The value of λ (0) k is the final value of λ k−1 , calculated in the previous time-step [10]. As mentioned earlier, the system is integrated by using an implicit single-step trapezoidal scheme [10]. In this approach, the relative joint velocitiesż and accelerationsz are corrected by using the mass-damping-stiffness-orthogonal projections as [3,9] whereż andz are, respectively, the relative joint velocities and accelerations obtained from the Newton-Raphson iteration, and W = M + t 2 C + t 2 4 K, where C and K represent the damping and stiffness contributions in the system. Note that in Eqs. (20) and (21), timedependent constraint terms are not incorporated because the constraints are assumed scleronomic. To account for rheonomic constraints, the reader is referred to [10,32].

Hydraulic actuators
In this study, the hydraulic pressures within a hydraulic circuit is computed by using the lumped fluid theory [37]. By using this approach, a hydraulic circuit is divided into discrete volumes where pressures are assumed to be equally distributed. The effect of acoustic waves is thus assumed to be insignificant. In a hydraulic volume V k , the pressure p k can be computed asṗ where Q ks is the sum of incoming and outgoing flows associated with the volume V k , n f is the total number of volume flows going in or out of the volume V k , and B e k is the effective bulk modulus associated to the volume V k . The effective bulk modulus can be written as where B oil is the bulk modulus of oil, n c is the total number of subvolumes V s that form the volume V k , and B s is the corresponding bulk modulus of the volume V s .

Valves
In this study, the valves are described by using a semiempirical modeling method [19]. By using the semiempirical modeling approach, the volume flow rate Q t through a simple throttle valve can be written as where p is the pressure difference over the valve, sgn( p) is the sign function that determines the sign of p, and C vt is the semiempirical flow rate coefficient of the throttle valve that can be calculated as where C d is the flow discharge coefficient, A t is the area of the throttle valve, and ρ is the density of the oil. Similarly, the volume flow rate Q d through a directional control valve can be written as where C v d is the semiempirical flow rate constant of the valve procured from the manufacturer catalogues, and U is the relative poppet/spool position. If the pressure difference is less than 2 bar, the volume flow is assumed to be laminar, and Eqs. (24) and (26) are modified so that the volume flow and the pressure difference follows a linear relation. Equation (26) is complemented by the following first order differential equation that describes a spool positionU where U ref is the reference voltage signal for the reference spool position, and τ is the time constant, which can be obtained from the Bode-diagram of the valve that describes the dynamics of valve spool.

Cylinders
The volume flow produced due to the motion of a hydraulic cylinder (shown in Fig. 3) can be written as where Q in and Q out are, respectively, the volume flow rate going inside and coming out of the cylinder,ẋ is the piston velocity, and A 1 and A 2 are, respectively, the areas on the piston and piston-rod side of the cylinder. The force F s produced by the cylinder can be written in terms of its dimensions and chambers pressure as where p 1 and p 2 are, respectively, the pressure on the piston and piston-rod side that can be calculated by using Eq. (22), and F μ is the total friction force caused by sealing.

Coupling of multibody formulations and hydraulic actuators
In this section, the multibody formulations described in Sect. 2 are extended to incorporate the dynamics of the hydraulic actuators described in Sect. 3 in a monolithic approach. The coupling of the double-step semirecursive formulation with the lumped fluid theory is inspired from [28] and [32]. Whereas the coupling of the penalty-based semirecursive formulation with the lumped fluid theory was already carried out in [32]. The force vectorQ in Eqs. (17) and (18) is incremented with the pressure variation equations, leading to the combined system of equations as follows where p is a vector of the pressures in the hydraulic subsystem and h (p, z,ż) are the pressure variation equations. It is assumed that the dependency of both Q and function h with respect to z,ż, and p are known.
In this study, both the approaches are integrated by using an implicit single-step trapezoidal rule [26], which is second order and A-stable method. While the trapezoidal rule was often used in structural dynamics, it was, however, seldom used in multibody system dynamics until the study by Bayo et al. [1]. In the study [1], it was agreed that the trapezoidal rule will lead to poor convergence characteristics when applied to multibody system dynamics in a similar way as other multistep integrators, that is, by considering the accelerations as primary variables. However, Bayo et al. [1] and Cuadrado et al. [8][9][10] demonstrated that the trapezoidal rule performs very satisfactorily when it is combined directly with the equations of motion by taking the positions as the primary variables, as shown below. Similarly, for the hydraulic subsystem, pressures are taken as the primary variables, as shown in [28,32]. Note that this study is more inclined to use the above approaches for real-time applications, such as [23], in future studies. Therefore, a single-step integration method is preferred that can use the same computational cost in each integration step [1]. Furthermore, even though the explicit, multistep integrators can be inexpensive and accurate, however, they do not demonstrate good stability, which is a limiting factor for real-time integration [1], especially for stiff systems. Thus, an implicit method is used. Moreover, A-stability is crucial for a numerically stiff system [1], such as presented in this study.
In the double-step semirecursive approach, the trapezoidal rule can be written as where t is the time-step, z i k are the independent relative joint coordinates,ż i k are the independent relative joint velocities,z i k are the independent relative joint accelerations, and p k andṗ k+1 are, respectively, the pressures and pressure derivatives. Equation (32) can be rewritten by considering z i k+1 and p k+1 as the primary variables and, respectively, solving forż i k+1 ,z i k+1 , andṗ k+1 at time-step (k + 1) aṡ Note that in the double-step semirecursive approach, the rules z i k+1 = z i k +ż i k t + 1 2z i k t 2 and p k+1 = p k +ṗ k t are applied to z i k+1 and p k+1 , respectively. In the double-step semirecursive approach, given z i k+1 , the dependent relative joint coordinates z d k+1 are obtained by iteratively solving the loop-closure constraint equations (z) = 0, which are highly nonlinear [11]. In this study, the Newton-Raphson method is used to solve the loop-closure position problem and convergence is achieved by providing a reliable estimate for z d k+1 by using information from the previous time-step as [15] The dependent relative joint velocitiesż d k+1 and accelerationsz d k+1 are computed from Eqs. (14) and (15), respectively. In the penalty-based semirecursive approach, the full set of relative joint coordinates z k+1 , velocitiesż k+1 , and accelerationsz k+1 are used in the above discussion of Eqs. (32), (33), and (34), instead of z i k+1 ,ż i k+1 , andz i k+1 . The dynamic equilibrium, established at time-step (k + 1), for both the approaches can be written as Equations (35) and (36) are a nonlinear system of equations that can be denoted as for the double-step semirecursive approach and x = z T , p T T for the penalty-based semirecursive approach. Such a nonlinear system of equations can be iteratively solved by employing the Newton-Raphson method as The residual vector [f (x)] (h) k+1 can be written as where λ is obtained as shown in Eq. (19). The approximated tangent matrix ∂f(x) ∂x (h) k+1 can be obtained numerically by using the forward differentiation rule where is the differentiation increment. To avoid ill-conditioning of the tangent matrix, is computed as in [32], namely where 1 × 10 −2 limits the minimum value for the differentiation increment to 1 × 10 −10 . Equation (41) is a modification of a method presented in [5]. In the penalty-based semirecursive approach,ż andz are corrected by using the mass-damping-stiffness-orthogonal projections [3,9] as shown in Eqs. (20) and (21).

The case studies of hydraulically actuated four-bar and quick-return mechanisms
In this study, a hydraulically actuated four-bar mechanism, as shown in Fig. 4, and a hydraulically actuated quick-return mechanism, as shown in Fig. 5, are used for a comparative study between the two multibody formulations in a numerically stiff coupled environment. The numerical stiffness in the coupled environment is introduced by the hydraulic subsystem. The mechanisms are modeled, first by using the double-step semirecursive formulation, and later by using the penalty-based semirecursive formulation, as explained in Sect. 2. For the planar system in Fig. 4, three joint coordinates are used in the modeling of the structure and two loop-closure constraints are used for a cut-joint (revolute joint) at point E. Whereas for the planar system in Fig. 5, five joint coordinates are used in the modeling of the structure and four loop-closure constraints are used for two cut-joints (translational joints) at points J and M. Both the mechanisms have one degree of freedom. In Fig. 4, bodies 1, 2, and 3 are assumed as rectangular beams, whose lengths and masses are    2 12 and mL 2 6 , respectively, where m is the mass and L is the length. The gravity is assumed to act in the negative Y direction, whose value is g = 9.81 m/s 2 . Both mechanisms are actuated by using the hydraulic actuators, as shown in Figs. 4 and 5. For simplicity, identical hydraulic actuators are used in both mechanisms. In this study, a simple hydraulic circuit is accounted, which consists of a pump, with a constant pressure source p P , a directional control valve, with a control signal U , a throttle valve, a double-acting hydraulic cylinder, connecting hoses, and a tank, with a constant pressure p T . The control volumes, V 1 , V 2 , and V 3 , used in the modeling of the hydraulic circuit are marked in Fig. 4. The pressure in the respective control volumes are p 1 , p 2 , and p 3 ; and their respective effective bulk modulus are B e1 , B e2 , and B e3 , that are calculated by using Eq. (23). For simplicity, the hydraulic circuit is assumed ideal, that is, the leakage is neglected.
In the hydraulic subsystem, the control volumes V 1 , V 2 , and V 3 are calculated as where V h 1 , V h 2 , and V h 3 are the volumes of the respective hoses; A 2 and A 3 are, respectively, the areas of the piston side and the piston-rod side within a cylinder; and l 2 and l 3 are, respectively, the lengths of the chambers of the cylinder, piston and piston-rod side. The length of the hydraulic cylinder is l such that l 2 + l 3 = l; and the variable chamber lengths, l 2 and l 3 , are calculated as where | s | is the actuator length of the hydraulic cylinder (see Figs. 4 and 5); s 0 is the actuator length at t = 0; and l 2 0 = s 0 − l and l 3 0 = l − l 2 0 are, respectively, the length of the piston and piston-rod side of the cylinder at t = 0. The differential equations of the pressures p 1 , p 2 , and p 3 are computed based on Eq. (22) aṡ where the volume flow rates Q d1 and Q 3d are calculated from Eq. (26), the volume flow rate Q t is calculated from Eq. (24), andṡ is the actuator velocity. The actuator length | s | and actuator velocityṡ of the hydraulic cylinder are computed as a function of the relative joint coordinates. For example, in Fig. 4, the upper end of the hydraulic cylinder is attached to body-3 at point F , while, the lower end is attached to ground at point G, such that r G = [0, −2, 0] T m. Therefore, s andṡ for the four-bar mechanism can be computed as where the position r F and velocityṙ F are calculated by applying the classical kinematic relations as in [11,22]. For simplicity, the force F s produced by the hydraulic cylinder in Eq. (29) is expressed in the form of Eq. (10) as where s X , s Y , and s Z are the components of vector s along the axes of the inertial reference frame. The initial value of the force F s 0 produced by the hydraulic cylinder is calculated from the static configurations, as shown in Figs. 4 and 5. For example, in case of four-bar mechanism, F s 0 = √ 2g (3m 1 + 2m 2 + m 3 ). In the static configuration, the initial value of pressure p 1 is equal to the initial value of pressure p 2 , which can be calculated based on Eq. (29) as p 2 = F s 0 + p 3 A 3 /A 2 . Note that the friction is neglected in static configuration and the initial value of pressure p 3 is assumed 3.5 MPa. The directional control valve, parameter U in Eq. (27), controls the movement of the cylinder through volume flows and is actuated for 10 s by using the following reference voltage signal U ref as t<1 s, 2.5 s ≤ t < 5 s, t ≤ 10, 10, 1 s ≤ t < 2.5 s, −10, 5 s ≤ t < 8 s, t<1 s, 4 s ≤ t < 6.5 s, t ≤ 10, −10, 1 s ≤ t < 4 s, 10, 6.5 s ≤ t < 8 s, where t is the simulation run time. The parameters of the hydraulic circuit are shown in Table 1. The set of variables used to solve the combined system of equations with the proposed integration scheme (Sect. 4) are Note that, in the double-step semirecursive approach (Eq. (48)), z i is identified by using the Gaussian elimination with full pivoting to the Jacobian matrix z . In this study, the error tolerance is considered as 1 × 10 −7 rad for positions and 1 × 10 −2 Pa for pressures. The voltage that corresponds to the spool position is integrated by using the trapezoidal rule and its error tolerance is considered as 1 × 10 −7 V. Furthermore, the penalty factor α in Eq. (18) is considered as 1 × 10 11 . Note that in the penalty-based semirecursive approach, the penalty term is analogous to a spring constant by considering that there is a spring attached to the cutjoint location to fulfill the constraints. Due to the high numerical stiffness introduced by the hydraulics, a high penalty term (spring constant) is used. In this study, both the approaches are implemented in the Matlab environment.

Results and discussion
This section presents the simulation results of the hydraulically actuated four-bar and quickreturn mechanisms presented in the previous section. Figures 6a and 6b show the simulation frames of the four-bar and quick-return mechanisms, respectively, presenting the position of the bodies at different instants of time. Here, the two approaches, namely, the double-step semirecursive approach and the penalty-based semirecursive approach, are compared based on the simulation work cycle, energy balance, constraint violation, and numerical efficiency.

The work cycle
In the four-bar mechanism, the hydraulic cylinder lifts the structure between 1-2.5 s and lowers it down between 5-8 s. Whereas in the quick-return mechanism, the hydraulic cylinder pulls the structure between 1-4 s and pushes it between 6.5-8 s. In the subsequent Fig. 7 The relative joint coordinates and pressures for the four-bar mechanism Fig. 8 The relative joint coordinates and pressures for the quick-return mechanism plots, the regions between the opening and closing of the directional control valve are highlighted in purple for the four-bar mechanism and in cyan for the quick-return mechanism. For both the mechanisms, the relative joint coordinates showed a good agreement in the two approaches, as shown in Figs. 7a and 8a. Thus, the solutions of both the approaches are accurate with respect to each other. The pressures in the hydraulic control volumes are shown in Figs. 7b and 8b, respectively, for the four-bar and quick-return mechanisms. Note that the pressures are identical in both the approaches and the choice of the multibody formulation does not affect the results of the hydraulics.
In the double-step approach, the independent joint coordinate is identified by using the Gaussian elimination with full pivoting on the Jacobian matrix of the constraints. For the presented work cycle, the independent coordinates are identified as z 3 for the four-bar mechanism and as z 4 for the quick-return mechanism, throughout the simulation. In the doublestep approach, if the independent and dependent coordinates are not identified adequately, then it can lead to the numerical problems during the integration because of the poor conditioned matrices. Therefore, the adequate identification of the independent and dependent coordinates is considered a relative drawback of the double-step approach compared with the penalty-based approach because the latter uses the full set of coordinates.

Energy balance
The energy balances in both the approaches are compared by analyzing the kinetic energy, potential energy, and work done by the actuator. The energy comparison of the four-bar mechanism is shown in Fig. 9 and of the quick-return mechanism is shown in Fig. 10. The energy balance in the double-step and penalty-based approaches showed a good agreement with each other for both the mechanisms. The peak energy drift in case of the four-bar mechanism is 0.84 J, which is 0.09% of the maximum actuator work, 982.95 J, and it occurred at the closing of the valve around 8 s. Whereas in case of the quick-return mechanism, the peak energy drift is 0.23 J, that is, 0.06% of the maximum actuator work, 357.47 J, and it also occurred at the closing of the valve around 8 s. An analogy can be drawn between the hydraulic system and a stiff spring, which supports the structure of the mechanisms.

Constraint violation
The basic difference between the double-step and penalty-based approaches is the way they handle constraints in their multibody formulation. Thus, it is worth showing a comparison of the constraint violations in both the approaches for the two mechanisms as shown in Figs. 11 and 12. Note that in the quick-return mechanism, two cut-joints are used in the modeling, but for demonstration purposes, the results for only one of the cut-joint are presented. In both the approaches, the constraints are fulfilled with good accuracy for both mechanisms. Thus, the robustness of the multibody formulations, as explained in the literature [10,24], is preserved in the application of a monolithic simulation of coupled multibody and hydraulic systems. Furthermore, the double-step approach fulfills constraints to the level of machine precision, or the precision specified, and it maintains a good error control. Therefore, this can be considered as the relative advantage of the double-step approach compared with the penalty-based approach because the latter approach is relatively more relaxed on the fulfillment of constraints.
In the double-step approach, for a possible mapping between independent and dependent relative joint coordinates, there is an assumption that redundant constraints and singular configurations do not exist. This assumption is a relative disadvantage of the double-step approach compared with the penalty-based approach because the latter approach can handle redundant constraints and can deliver accurate solutions in the vicinity of singular configurations [15].

Numerical efficiency
The numerical efficiencies of both the approaches are compared in Figs. 13 and 14, for the four-bar and quick-return mechanisms, respectively. In the four-bar mechanism, the average and maximum iterations, and the total integration time for the double-step approach are 1.39, 3, and 27.45 s and for the penalty-based approach are 1.56, 4, 21.54 s. Whereas in the quickreturn mechanism, the average and maximum iterations, and the total integration time for the double-step approach are 1.49, 3, and 29.03 s and for the penalty-based approach are 1.57, 6,   Even though the number of iterations are lower in the doublestep approach, its integration time is greater compared with the penalty-based approach, which is in accordance with the literature [15,27]. The poor numerical efficiency of the double-step approach is attributed to the iterative solution of the dependent joint coordinates by using the Newton-Raphson method.
Note that the examples are implemented in the Matlab environment and for this reason, the result about the integration time may be unreliable because of the implementation details. However, both approaches are carefully implemented such that they have the same number of function calls in each iteration. This was achieved by symbolic precomputation and function generation of the matrix product, as in [32]. Furthermore, both approaches can be potential candidates for real-time applications by implementing them in a lower level language such as C++ or Fortran and by limiting the number of iterations to a lower value.

Conclusion
This study introduced the double-step semirecursive formulation and compared it with the penalty-based semirecursive formulation, in a numerically stiff coupled environment. The numerical stiffness was introduced by a hydraulic system. A monolithic scheme for the coupled simulation of the double-step semirecursive formulation and hydraulic systems was introduced in this study. To this end, the hydraulic system was described by using the lumped fluid theory. As case studies, hydraulically actuated four-bar and quick return mechanisms were modeled to compared the double-step semirecursive approach with the penalty-based semirecursive approach. The two approaches were compared based on the work cycle, energy balance, constraint violation, and numerical efficiency of the mechanisms, and similar results were obtained in both the mechanisms.
The relative joint coordinates and energy balances showed a good agreement in both the approaches. In the double-step approach, the adequate identification of the independent and dependent joint coordinates was carried out by using the Gaussian elimination with full pivoting on the Jacobian matrix of the constraints. Its identification is considered a relative drawback compared with the penalty-based approach where the full set of coordinates were used. In both the approaches, the constraints were fulfilled with good accuracy. However, the double-step approach had an advantage of fulfilling the constraints to the level of machine precision, or the precision specified compared with the penalty-based approach. Nevertheless, the double-step approach had an assumption that the redundant constraints and singularity configuration do not exit. This is a relative disadvantage of the doublestep approach because the penalty-based approach can handle redundant constraints and can provide accurate solutions in the vicinity of singular configurations. Furthermore, the double-step approach suffered from poor numerical efficiency, which was attributed to the iterative solution of the dependent joint coordinates by using the Newton-Raphson method. To improve the numerical efficiency of the double-step approach, it would be of utmost importance to enhance the iterative solution of dependent coordinates and this is left as a topic for future studies. In conclusion, the penalty-based semirecursive approach has a number of advantages over the double-step semirecursive approach, which is in accordance with the literature [15,27].
For future studies, alternate multibody formulations can be coupled with hydraulic systems to study on the optimal formulation for simulating coupled multibody and hydraulic systems. Note that the selection of the integrator type is bound to have an effect on the results. In this study, an implicit integrator is used, whereas the double-step semirecursive formulation is usually integrated by using the fourth order Runge-Kutta method, that is, an explicit integrator. Therefore, different integrators can be implemented to study on the optimal integrator choice for such coupled systems. This study utilized the Matlab environment, therefore, to make a firm conclusion on the numerical efficiency of the two approaches, a large-scale, three-dimensional example needs to be investigated in programming languages, such as C++ or Fortran.