## 1 Introduction
[7] combined the MBD solution and smoothed particle hydrodynamics solution to calculate the fluid-structure interaction. Moreover, another way to generate an MBD model with EHL is through monolithic modelling. Monolithic modelling, sometimes called uniform, is based on the representation of all subsystems of a multiphysical problem using the same time integrator [8]. [PERSON] et al. [9] investigated two coupling techniques: those in which only one tool performs integration (function evaluation), and those where each tool uses its own integrator (co-simulation). Function evaluation means that one software is performing numerical integration, while the other software returns values on request, as per the states passed to the integrator tool. This type of coupled simulation cannot be classified as co-simulation since only a single integrator performs the integration. This technique is popular in expanding the original simulation code with additional phenomena, which is originally not included. The MBD software with a prepared EHL solution is commercially available. Often, the software is designed to be used in a specific industry, e.g., internal combustion engines [1; 10; 11; 12; 13]. These types of software offer fast computational time, though they are unable to solve problems of dynamics in other industries. Commercial general-purpose MBD software often do not include a solution to EHL. Therefore, it is necessary to create a coupled simulation in one of the above-mentioned ways. Another large group of MBD software with EHL are scientific software that have been developed at universities or research organizations. Generally, this type of software is intended for a specific group of users, as often there is no user-friendly interface or user manual [14]. The aforementioned MBD programs extended by EHL differ in their approaches to solving this phenomenon. The simplest approach is to use empirical approximation relationships that are designed for a given application. [PERSON] et al. [15] applied such a simplified solution to the elastohydrodynamic contact of the teeth in a geared multibody system. A more advanced lubrication solution considers the premise of an infinitely short or infinitely long bearing. As per this assumption, it is possible to derive an analytic expression for the hydrodynamic pressure of a bearing [16]. Nevertheless, in the real world, the bearings have finite dimensions. Therefore, there is a more complex numerical solution to EHL that requires that the Reynolds equation to be solved either directly or iteratively. The Reynolds equation is numerically discretized, such as by the finite difference method [17], the finite volume method [18] or the finite element method [19]. Available MBD models with EHL differ according to the method of calculation of elastic deformation. [PERSON] et al. [15] calculated the tooth deflection of the geared multibody system as the superposition of bending deflection (nonuniform cantilever beam), the deflection of the fillet and foundation, and the local contact and compression deflection (linearized Hertzian contact model with constant stiffness). At this point, it is worth noting that the Hertz contact theory [20] is adequate for the concentrated contact between two ellipsoids. For concentrated line contacts, as between two cylinders, the Hertz theory does not afford such a relationship, and the formulae for ellipsoids are invalid [21]. [PERSON] et al. [17] calculated the bearing deformation using a simplified Winkler model [22], where the deformable body is modelled as a set of spring elements neglecting the shearing action. The actual deformation at a point is only dependent on the contact pressure at that point, the material and the geometric parameters of contact bodies. The component mode synthesis is a universal and efficient method for calculating deformation in the MBD system.