## 4 Conclusions  
In some cases, the solution to EHL is required to properly assess the functionality and performance of the design of the mechanism, and examples are not hard to find. An internal combustion engine has several components, such as a piston group, cam followers, gears, etc., where EHL plays a critical role. To analyze these components, general-purpose MBD software is often used. However, not every general-purpose MBD software has a prepared option to solve the advanced lubricated contact pair. Therefore, it must be augmented with these physical effects. To keep the simulation within one software, monolithic modelling is advantageous. In this case, the EHL is implemented in the MBD solver using the user-written subroutine. This brings challenges that may be solved as follows.  
During the solution to the MBD simulation, the overlapping of structural and hydrodynamic grids occurs. Thus, data from one grid must be mapped onto another grid. To transfer the deformation and true shape data, a cubic Hermite spline is a precise option.  
A cubic spline is not sufficient to map the force effects and leads to solution instabilities. On the contrary, mapping by analytical integration respects the equivalence between the source and target mapped field and is suited for mapping force effects. This type of mapping is recommended when different structural and hydrodynamic grid sizes are used. This is advantageous as the size requirements of these grids differ heavily. The deformation is well represented by the coarse mesh. However, the hydrodynamic solution requires a fine mesh. In this case, where the same fine mesh is used for the structural grid, the computational time becomes excessive.  
The introduction of EHL in the MBD model causes an increase in the number of degrees of freedom, thereby causing an increase in the computational time. If the sensitivity analysis is solved serially by the finite forward difference method, the computational time may be unacceptable. If the MBD solver allows for the use of a more advanced method, such as the adjoint variable method, the direct differentiation method or the automatic  
Figure 11: Piston thrust side minimum oil film thickness dependent on the crankshaft angular position.  
Figure 10: Piston lateral displacement dependent on the crankshaft angular position. The firing occurs near a crank angle of 360\({}^{\circ}\).  
differentiation method, it is highly recommended to do so. The _Parallel_ method proposed and verified in this paper reduces the computational time by an order of magnitude.  
The parallelization of the design sensitivity analysis requires minor changes to the subroutine's code. The accuracy of this method is determined as per the perturbation step size used. In general, the optimal perturbation step size changes during the MBD simulation. Nevertheless, it does not change significantly for the measurements related to the deformation. Therefore, in this case, the perturbation step sizes estimated at the beginning of the simulation can be used for the rest of the simulation without significantly reducing the accuracy of the Jacobian.  
Conceptualization, J.D.; methodology, J.D.; software, J.D.; validation, P.N.; writing--original draft preparation, J.D.; writing--review and editing, P.N. All authors have read and agreed to the published version of the manuscript.  
The research leading to these results has received support from Specific research of Brno University of Technology, project No. FSI-S-20-6267.  
Not applicable.  
Data sharing is not applicable to this article.  
The authors declare no conflict of interest.