Benchmarking of two flexible multibody dynamic simulation software in engine simulations

In this paper, two different commercial multibody dynamic (MBD) simulation software cases are studied. Due to the restrictions determined in the conditions of contract, the names of the software are not revealed, instead being called Software S and Software E. The central purpose of this research was to investigate the abilities of Software S in the simulation of a large engine, as a part of the strength analysis process. The abilities were studied by comparing the program with another, here called Software E, which is designed primarily for engine simulations. The capabilities of Software E have been proven after years of usage at Wärtsilä, resulting in its essential role in the strength analysis process today. The aim was to find the shortcomings and restrictions of Software S but also advantages it could bring to the strength analysis process for Wärtsilä. Similar simulation models were also built using both programs during this research. A 16-cylinder V-engine was selected as the subject because of its size in order to obtain further information about the behavior of the program when working with extensive model files. The components of the engine were flexible and were reduced FE models, also called super elements. The forces and contact situations that occur inside the engine were modeled using elements provided by the MBD programs. Different levels of detail of the modeling elements were used to obtain information about the flexibility of the program. The results obtained from time integrations were compared to ensure the similarity of both modeling elements used. Also, this paper reports the calculation times. In addition, a small-scale study was performed for Software S to clarify the effect of the modes used in time integrations towards results accuracy and calculation times. Simulation models were built successfully in both programs, and the results obtained correlated with each other on an adequate level. Significant differences appeared in the features and usability of the programs in general. The GUI of Software S is advanced and user-friendly, whereas Software E is not focused on these features. On the other hand, the modeling element library of Software E covers all of the required features related to large engine simulations, some of which Software S is lacking. This work can be used in assistance when considering buying new software for a company as well as when investigating new development areas that could be improved with new software.


Introduction
Benchmarking is a practice widely used in business and can be said to be a key component for organizations in looking for improvement in quality, competitive position, or market share.Benchmarking is a continuous improvement process called "striving to be the best of the best" [1].In this article, the term "benchmarking" is used for the qualitative comparison between two simulation programs.Due to the scope of the comparison, methods of qualitative research have been used.Qualitative research allows for an open and flexible design and, in doing so, stands at odds with the rigorous guidelines of quantitative research [2].
In many mechanical systems, such as the reciprocating engine, the components undergo large rotational and translational displacements and produce kinematic constraints for each other with different types of joints.The definition of a multibody dynamic system is a collection of subsystems called bodies and components, each subsystem possibly being subject to substantial rotational and translational displacements.[3] Computer-based simulations of reciprocating engines have developed along with computers, software, and calculation techniques.[4] At Wärtsilä, engine development is based on highly detailed computer-based simulations of a running engine.The results obtained from simulations correspond very well with results measured from a running engine, as can be seen in multiple articles [5][6][7][8].However, there is still a need to speed up the development process at Wärtsilä [9,10], and as an additional challenge, [11] to show the steadily increasing number of customer-specific variations of product configurations.Further, the need to capture different type of physical phenomena [12,13] sets more and more requirements for the software.
In this work, the suitability of the MBD simulation program Software S is investigated as a part of the engine development process at Wärtsilä.In practice, it is examined whether Software S could replace the MBD simulation program, Software E, which has been used at Wärtsilä for years.The aim is to find the advantages and disadvantages that Software S could bring to the engine development process.When comparing the capabilities of the programs, the focus was on usability in general, and the performance, available modeling elements, and accuracy of the results in particular.
This paper shows MBD simulation models of the Wärtsilä engine modeled using both programs and the results of their comparison.To obtain comparable results, the Structural Analyst built the models to be as similar as possible.The Abaqus analysis produces the reduced FE models, which describe most of the engine components.A reduced FE model contains information about the shape functions of the body as well as the reduced mass and stiffness matrices.Different levels of detail of the modeling elements provided by the MBD software determined the interactions between bodies.
Because the verification of the simulation models was not the main focus of this study, and based on the previous verifications for Software E, the simulation model done by Software S was verified by the results obtained from Software E. Again, the parameters needed for the simulation model in Software E are based on the simulation models verified earlier by measurements [5][6][7][8].

Generic functionality of the MBD software
Software S, as well as Software E, are widely used within a wide range of industry areas such as automotive, power transmission, and engines.They are capable of handling investigations and full non-linear analyses of static, quasi-static, and dynamic systems over a complete frequency range.However, Software S has focused on more general multibody simulations, whereas Software E has been more tailored to combustion engines, transmissions, and powertrains.
The main workflow is similar for both programs and is outlined in figure 1.At the first stage, the input files have to be generated to define flexible body definitions.A common reducing technique for the coordinate reduction, the Component Mode Synthesis technique (CMS), is used.The CMS method assumes that a specific set of modes can describe small elastic deformations in the body.The reduction process results a reduced Finite Element model (FE-model) called a "substructure".The Structural Analyst used Abaqus 2016 to obtain the substructures used in this work.Many different methods have been developed based on the CMS technique, e.g., the fixed interface method (Craig-Bampton [14]) and the free interface method (Craig-Chang) [15].The method chosen is dependent on the boundary conditions used and are explained in more detail [16,17].
To account for the flexibility of the bodies in the multibody dynamic system, the same principle is applied for both software, i.e. the "floating frame of reference formulation".Two types of the equation of motion can be distinguished: an equation of motion for the large motion of the reference body and vibration equations for the small elastic vibrations/deformations. [18] After this, the Structural Analyst used the GUI of the multibody simulation program to build the actual simulation model.Also, the Analyst defined the bodies, their connections to each other, contacts and external excitations of the model.The third step is setting the time integration parameters such as integration settings and tolerances.At this stage, the Structural Analyst also defined the desired results output of the result file.The final stage, resulting from post-processing and exporting, can be performed in a separate post-processor or by exporting the results to a third-party program in an appropriate file format.

Simulation model
Table 1 shows the components used in the simulations.Flexible bodies are defined in the substructures, and rigid bodies were modeled using a few input parameters to obtain their mass and inertia properties.Eigenmodes selected to be included in the substructures are defined depending on the scope of the analysis.In this case, the engine block includes eigenmodes up to 200 Hz other components -up to 2000 Hz.Due to the strong dynamic behavior of the reciprocating engine, the substructures have to include both the static and the dynamic properties of the full FE-model.Therefore, the substructure generation is a two-stage process.The first step is a natural frequency extraction, and the second is the actual substructure generation analysis.Natural frequency extraction is performed in order to include dynamic behavior, whereas static behavior is included in the substructure generation analysis by defining the degrees of freedom for the retained nodal points.In this simulation model, degrees of freedom of the substructure retained are not constrained in the frequency extraction step, therefore, the Craig-Chang method was utilized [15].In table 2, the modeling elements which were used to describe the contacts and excitations in the virtual engine are presented.The aim was to use similar modeling elements in both software.However, it can be seen in the table that differences exist in piston liner contact and the bearing descriptions.Loading conditions of the engine were a constant torque added to the end of the crankshaft, corresponding the torque generated by the cylinder pressure.This approach does not fully describe the real situation of the loading of the engine fluctuating highly due to the variation of the electric grid conditions.The effect of the loading variation on the dynamic engine behavior is a complicated process, which is discussed more in the article.[19]  In Software S, a similar piston liner contact was not available as in Software E, which will probably have a small effect on results.Thus, the Analyst can use constraints in Software S, meaning the elimination of the x-direction movement of the big end of the conrod relative to the cylinder liner (Figure 2).For the constraint points, the script performed linear interpolation of markers from the four nodes on the cylinder surface.By doing so, it accounts for the effects of the side forces due to the cylinder pressure;

Force pair between engine block and conrod Camshaft excitations
Force and moment into the camshaft Force and moment into the camshaft however, not for the deformation of the cross-section of the cylinder because of the linear interpolation of the constraint markers from the four nodes located on opposite sides [20].
In Software E, the side forces are added only to the contact area of the cylinder wall.
In Wärtsilä engines, intermediate gears transfer the torque from the crankshaft to the camshafts.Due to the torsional vibrations and fluctuating loads, gear forces are usually non-constant and include high impacts.These impacts will produce high-frequency vibrations, which again produce structure-borne noise and may potentially introduce conditions for the fretting phenomenon [5].An accurate simulation of the gear pairs is needed to obtain stable dynamic behavior of the reciprocating engine.[21] The Structural Analyst modeled two different detail levels of the gear pairs in both programs.The simpler modeling approach is similar between programs and does not consider tooth deformations.The more advanced gear wheel model also takes tooth deformations into account.Figure 3 shows the retained nodes of the advanced gear wheel in Software S. Teeth forces are calculated based on the stiffness of these nodes with the help of analytical calculation of the compliance of the teeth.In Software E, the advanced gear wheel is described purely analytically, and the nodes on the teeth roots are not needed [18].
Real Wärtsilä engines are equipped with an injection pump, which are used to pump fuel into the cylinder.Each cylinder has its own injection pumps, driven by a camshaft such as exhaust and intake valves.Due to these, they produce forces and moments on the camshaft, which are needed to take account of and to capture their effect on the results.In order to keep the complexity of the model as low as possible, these kinds of forces and moments are modeled with the help of excitations.Also, the cylinder pressure, which is the principal source for the variable mechanic load, has to be modeled.[22] To define cylinder pressure, there are ready add-on modules in both programs.In Software E, the workflow for defining other excitations according to the crank angle is very straightforward.However, in Software S, the definition of the excitations according to the crank angle requires more effort.Figure 4 shows the work stages needed in defining the excitations.To connect to the instantaneous crank angle method, an expression element (2.) is used here to measure the crank angle during the solver run.This is used again as an input in a new expression (3.) with the input function (1.), which contains the excitation force as a function of time.The output signal from this expression element (3.) is input into a force element assigned to the node on the rotating body.

Elastohydrodynamic bearing
The description of the hydrodynamic bearing phenomena plays a significant role in MBD simulations at Wärtsilä.There is a demand for higher load, higher reliability, and lower losses for engine bearings.[23] Bearing phenomena are also a central interest when analyzing crankshaft dynamics.Furthermore, bearing results are boundary conditions in non-linear finite element analysis, [7,8] which is the reason for the extensive study of the EHD modules in this paper.A hydrodynamic bearing is the entire or partial separation of the sliding surfaces by a hydrodynamic lubricant, where the hydrodynamic pressure in the lubricant keeps the sliding surfaces separated from each other [24].Based on the theory of the hydrodynamic and Reynolds equations, this hydrodynamic pressure results in low friction and theoretically zero wear.[25] In both software, the Reynolds equation determines the hydrodynamic behavior of the oil film on the bearing surfaces.Usually, the most important results over time are oil film thickness and pressure distribution in the oil between sliding surfaces.When applying the EHD module, solving the equations of motion takes into account the local elastic deformations of the bearing surfaces.In Software S, the accuracy of the solution depends strongly on the flexible body modes used in time integration.[20] The main difference of the hydrodynamic bearing modules between programs was the distribution of the bearing forces.To be able to use a similar bearing model, the Structural Analyst needed to use a different type of node connectivity.In Software E, the resultant load of a single bearing module was distributed between five center nodes on the crankshaft, while in Software S, the Analyst could only use one center node.Also, again in Software E, the resultant load to the engine block distributed to the surface nodes of the bearing surfaces, while in Software S, the distribution was to the one center node.Figure 5 presents the description of the main bearing in Software E and Software S. Red discs in the two first pictures presents the center nodes on the crankshaft (in Software E, five center nodes are used for each main bearing instead of one in Software S).Again, below in the Figure 5, red nodes present the connectivity nodes on the engine block.In Software E, 200 nodes on the bearing surfaces were used for the bearing load distribution, whereas in Software S, one center node is defined, which is coupled to 200 nodes located on the bearing surface.The EHD results produced by Software E have been noted to correspond well with the real situation, and the Structural Analyst will also obtain the local elastic deformations of the bearing surfaces in a reasonable calculation time.In contrast, there is a need of at least two hundred flexible body modes in a substructure for each EHD module to describe the local deformations in Software S correctly.This number of modes will increase the calculation times dramatically, especially for large engines that can contain dozens of sliding bearings.

Time integration
Because of the complexity and non-linearity of the dynamic multibody system, time integration has to be performed in the time domain.In both software cases, the implicit integration method was used to solve the equation of motion, the Backward Differentiation Formulation (BDF).In Software S, BDF is the default integrator, whereas in Software E, it is the only option available.BDF interpolates the solution of k of previous time steps t i+1−j with j = 1, ..., k and, at the next point of time, discretization t i+1 .For this purpose, the polynomial degree k needs to be found, which also denotes the BDF order.[18] In Software E, the possibilities to modify the solver settings are shown in Figure 6, as well as parameters used for those.In Software S, the solver settings are customizable and much more widespread, such as the order of the BDF, and a formalism between explicit and residual can be changed.In this work, default solver settings for Software S were used, except the BDF order, which was changed to two from five, as in Software E.

Results
This section presents some results from time integrations of both programs.The simulation models used were as similar as possible, except for the differences presented earlier, and the same modes were used (Mode set 1 and Mode set 2, in Table 3) to describe the deformations.
Table 3 presents four different simulation cases.Mode set 1 presents Software E, and the other three represent Software S. Taking into account the camshaft excitations, the engine block was simulated as rigid (without any modes).Without the camshaft excitations, the Structural Analyst used a flexible engine block model.This paper shows this separation because of integrator stability problems present in Software E. Software S simulation times are strongly dependent on the amount of modes considered, and therefore, the Structural Analyst performed the simulations using different combinations of modes (both eigenmodes and IRM).The simulation times, integration steps and iterations can be also seen in Table 3.The number of iterations in Software S is higher, which is a consequence of a higher number of integration steps.In both cases, the integration step size is adjusted automatically by the solver in order to obtain a solution within given tolerances.By adjusting tolerances in Software S, the number of integration steps might decrease but could not be included due to the time constraints of this study.The time integrations were run for six full engine cycles, which means twelve full crankshaft rotations.In Software S, a special force element monitors and stops time integration, when the crankshaft reaches six full rotations.
Figures 7 and 8 present the results from the hydrodynamic bearing module located in the main-bearing five.Both figures show the correlation between programs.The behavior of the curves is cyclic, they are quite similar, and the magnitudes generally correspond to each other reasonably well.However, both figures show small difference.One reason for the differences undoubtedly originates from the force distribution of the hydrodynamic bearing module, as described in the previous section.Because of this, the base of the bearing modules in Software E is five times stiffer than in Software S, and the forces are distributed differently between the other main bearings.Figure 8 shows that highfrequency vibrations occur more often in Software E results, which may also indicate a stiffer bearing module.
Figure 9 presents the maximum force in the axial bearing.The Structural Analyst applied a surface-to-surface method and used EHD modules in both programs.The bearing forces calculated by Software S are approximately two times higher than in Software E. Going by previous experience, Software E results seem to be more correct.After careful study, the Structural Analyst could not find an apparent reason for this difference, but it is sensible to simulate the axial bearing module with a long crankshaft, and even a tiny modification of some parameter can change the behavior of the axial bearing.With this in mind, there are many parameters related to the axial bearing, defined differently between the programs.In order to obtain more similar results, the parameters should be adjusted further and the effect of those to the results checked.
Figure 10 shows the total moment produced by the gear wheel on the camshaft.High impacts and fluctuation of the tooth forces due to the hammering of the gear wheels within the clearance can be seen in the figure.The modeling of the gear pair is sophisticated and requires a lot of different parameters.Some differences appeared between them, like in Software S, the "Addendum Modification Factor" is called the "Profile Shifting Factor" in Software E [26].For Software E, these parameters have been verified earlier [21] and utilized in this simulation model.When using flexible bodies, retained nodal displacements are the primary results of the simulation and are used directly in post-processing such as stress analysis.Thus, the relative displacements between the nodes enable the study of the deformation of the body.With relative node displacements, the web opening can also be measured, providing a good description of the stresses in the fillets of the crank pin.The main reason for the web opening is the cylinder pressure force acting on the crank pin, see Figure 11, defined by measuring the relative displacement of two adjacent main bearings in the axial direction.[27] Figure 11.Simplified web opening on the crankshaft due to the cylinder pressure Figures 12 and 13 show the web opening for two webs (between main bearings 1-2 and 8-9).The figures show that the noticeable differences between the programs appear for the first web opening, whereas the difference is getting smaller between main bearings four and five and even smaller between eight and nine.Taking into account the different results of the axial bearing between the programs, it becomes clear that the axial bearing is the reason for the differences in the web openings.The author also evaluated the user experience of the programs.The Structural Analyst used both programs via their graphical user interface (GUI), and therefore, its functionality and user-friendliness were also significant focuses of the research.
There are some advantages and disadvantages in both programs.In Software S, the big size of the simulation model lead to the slowing issues with model handling in the GUI.Thus, the stress recovery is problematic for big components.Regardless of the issues, the GUI of Software S is modern and sophisticated, and model creation is easy even for an inexperienced user.On the contrary, the GUI of Software E contains many features that are unique to it, and the user would require training or usage of a significant period of time to use the software without issue.
In general, the creation of a single modeling element is effortless in both cases.However, in a big engine, there are multiple modeling elements of the same type, acting in different nodes, for example, the bearings of the connecting rods (sixteen instances) or the bolted oil sump.In Software E, this kind of repetition is partially automated (the automatic searching of nodes for the conrods), and it is possible to assign one modeling element between multiple nodes (for example, one force element to a hundred bolt connections).In Software S, the search of node numbers is not automated, and for every single connection, the Structural Analyst had to create a separate modeling element.The solution could be automated scripting, but the model created manually requires more work in Software S compared to Software E, regardless of the superior user-friendliness of Software S. On the other hand, one advantage of software S is the rigid joint that is easy to define and usually needed when modeling, for example, stiff bolted joints.In Software E, the Analyst needed to use stiff springs for this purpose, which often leads to convergence issues.
An advantage of the MBD approach is that the calculated results are valuable boundary conditions for more detailed FEM analyses.Hydrodynamic pressures, acceleration fields, nodal displacements and forces can be transferred automatically to FEM models using Software E. This capability is missing from Software S, although it is possible to create scripts as a workaround.On the other hand, software S is more a general-purpose tool and contains some features that could be utilized directly in the MBD model instead of creating a separate FEM model.These include, for example, dry frictional contacts, which could be useful for bolt joints.
In Table 4, the main advantages and disadvantaged for both software are summarized.

Conclusion
The central purpose of this research was to investigate the abilities of one MBD simulation program, Software S, in the simulation of a large engine, as a part of the strength analysis process.The ability was studied by comparing the software with another MBD simulation program, called Software E, designed primarily for engine simulations.The capabilities of Software E have been proven as a result of years of user experience at Wärtsilä, today an essential part in the strength analysis process.The aim was to find the shortcomings and restrictions of Software S, but also advantages it could bring into the strength analysis process for Wärtsilä.Similar simulation models were built using both programs, Software S and Software E, during this research.A 16-cylinder V-engine was selected as the subject of this work because of its size and also to obtain information about the behavior of the program using extensive model files.The components of the engine were flexible, reduced FE models, also called "superelements".The forces and contact situations that occur inside the engine were modeled using elements provided by the MBD programs.Different levels of detail of the modeling elements were used to obtain information about the flexibility of the program.The results obtained from the time integrations were compared to ensure similarity of both modeling elements used.Furthermore, this paper reports the calculation times.In addition, a small-scale study was performed for Software S to clarify the effect of the modes used in time integrations toward result accuracy and calculation times.The simulation models for the big engine were built successfully using both programs and the results obtained correlated with each other on an adequate level.
In general, significant differences appeared in the features and usability of the programs.The GUI of Software S is more advanced and user-friendly, whereas Software E is not focused on these features.However, it still has some shortcomings compared to Software E, in relation to MBD simulations of the big engine.Perhaps the influential restrictive factors of Software S worth mentioning are the problems with hydrodynamic as well as the elastohydrodynamic bearing modules.Also, multi-assigning modeling elements is one of the most desired features, which the software does not have.Software S lacks an interface for FEM solvers to transfer loadings, such as hydrodynamic pressures, to make further, more detailed FEM analyses, while Software E has an excellent interface with different FEM solvers.To summarize, Software S is a more general-purpose tool, thus, it could possibly solve dry frictional contacts on the MBD side.Software S also has a rigid connection, which would be necessary for stiff bolt connections.Possibly future research could be performed to investigate the effect of different mode combinations in time integration regarding calculation times and result accuracy further.Also, the possibility to determine a unique damping value for a single frequency would require more research.Some potential features might also be co-simulation with Software S and Abaqus and different contact situations based on surface geometries, such that the Structural Analyst could perform a more realistic simulation and more accurate analysis of a big reciprocating engine.

Figure 2 .
Figure 2. Interpolated cylinder marker for constraint forces

Figure 3 .
Figure 3. Teeth nodes on the root of the intermediate gear

Figure 5 .
Figure 5. Main bearing approach in software E (left) and in software S (right)

Figure 6 .
Figure 6.Integration settings in software E

Figure 10 .
Figure 10.Simple gear wheel, moment on the camshaft

Table 2 .
Used modeling elements

Table 3 .
Integration statistics and mode sets used in time integration

Table 4 .
Comparison between programs