Abstract
This paper addresses the development of an elastodynamic model of a motorcycle engine cranktrain aimed at accurately evaluating the interactions between the crankshaft and the engine block, thus allowing an improved structural design. A rigid multibody model is first implemented and simulated; only kinematic joints are involved at this stage, leading to a statically determinate assembly of the mechanism. Such a modelling approach prevents the loads at certain interface locations to be evaluated; furthermore, high-frequency dynamic effects cannot be predicted. These drawbacks can be removed by introducing bushing-like elements and/or modelling component flexibility. In this paper, this latter aspect is the objective of the investigation; in particular, a finite element model of the crankshaft is implemented as a replacement for the corresponding rigid member. The well-established Craig-Bampton model reduction technique is used to represent the elastodynamic behaviour of the component with a limited number of coordinates. The mode selection procedure is emphasized here: a measure of modal dynamic importance, namely the effective interface mass fraction, is used to rank fixed-interface normal modes based upon their contribution to loads at the substructure interface; choosing the modal base according to such ranking leads to a minimal yet accurate representation.
1. Introduction
Modern powertrain design is facing increasingly strict requirements in terms of emissions, fuel consumption, noise, and vibration levels. In recent years, this trend is extending towards the motorcycle industry, in which competitive design focused on achieving a high power-to-weight ratio calls for optimized engine components. This in turn requires the adoption of a multidisciplinary approach early in the conception phase, and the use of advanced simulation tools which might help the analyst in gaining a deeper insight into the physical phenomena associated with the engine operation. Concerning powertrain dynamics, modern analysis techniques involve the use of multibody simulation tools, which allow an accurate prediction of the operational loads acting on the engine components, leading to an optimal structural design.
Several approaches are described in literature dealing with multibody modelling of internal combustion (i.c.) engine powertrains. Some papers deal with the construction of fully coupled cranktrain models through the use of commercial multibody dynamics codes, which provide a general modelling platform for mechanical systems, see, for example, [1, 2]; the equations of motion are in this case implicitly synthesized by the software kernel, and solved by means of some standard integration scheme. As an alternative, some studies describe the development of specialized modelling codes, see [3, 4]; the system equations of motion are retrieved analytically, converted to computational code, and solved numerically.
In this work, the former approach has been followed, the model being developed by means of ADAMS multibody simulation software; this package offers standard performance in terms of results’ accuracy versus simulation time, also allowing for integration with other simulation tools, for example, finite element codes, and a straightforward procedure in defining customized subroutines. These aspects play an important role in the perspective of developing an integrated, multidisciplinary i.c. engine simulation platform, which is the ultimate goal of the research framework this work is part of.
In the context of multibody modelling, the definition of a system made up of rigid links connected to each other via kinematic joints typically represents the first step in the process; in fact, commercial multibody software platforms offer both CAD interfaces and joint libraries, permitting the analyst to set up a basic dynamic model with little time and effort. Clearly, this modelling approach is affected by some important limitations: first of all, the adoption of rigid bodies in combination with kinematic joints prevents some interface loads to be evaluated whenever the mechanism under study exhibits some static indeterminacy; furthermore, any structural dynamic effect, which might contribute significantly to the actual loads, is evidently lost. These shortcomings can be eliminated by embracing a refined modelling methodology, based upon the introduction of nonideal joints at interface locations, for example, bushings, and the inclusion of component flexibility. The latter, in particular, is a necessary task if one wants to capture dynamic effects which arise in high-speed applications, and is the main concern of this paper.
In flexible multibody dynamics, the floating frame of reference formulation is the most widely used method to account for component flexibility [5]; according to such approach, the motion of each flexible body is subdivided into a reference motion, which can be described according to rigid multibody formalism, and a deformation. The finite element method is used to describe such a deformation; since small displacements and rotations are assumed to occur about the floating frame, linear methods borrowed from structural dynamics can be exploited, such as modal superposition techniques. These techniques are used to reduce the number of coordinates required to describe the component deformation; this is accomplished by means of some coordinate transformation, which normally involves the definition of proper component mode shapes. Several component mode synthesis methods have been proposed in the past [6–9], the Craig-Bampton approach being the most popular amongst the multibody dynamics community. The coordinate transformation matrix is built up in this case starting from the computation of constraint modes and fixed-interface normal modes; proper selection of the latter is a nontrivial task, and several papers exist on the subject, see, for example, [10]. In this work, the selection procedure is carried out in accordance with a modal ordering scheme based on the effective interface mass measure of dynamic importance [11, 12], which ranks modes based upon their contribution to loads at component interface.
This paper discusses the first steps in the development of a multibody model of a Ducati L-twin, four-stroke engine cranktrain. The ultimate purpose of such model is the accurate prediction of the loads acting on the main components in the system, with special regard to the engine block; this will eventually enable a refined structural design, resulting in a weight reduction combined with improved overall performances. After a brief description of the system under study (Section 2), the modelling process will be reviewed; in Section 3, a rigid multibody model is presented, and inherent limitations are highlighted; Section 4 depicts the flexible multibody implementation, mainly focusing on the adopted methodology. Some concluding remarks and future research directions are eventually given in Section 5.
2. System Description
As mentioned, this study deals with the cranktrain of a motorcycle i.c. engine. Figure 1 shows a schematic of the mechanism. The main component in the system is a single throw crankshaft supported by four main journals, three of which equipped with ball bearings (the two central ones and the external main bearing 1), the other one being hydrodynamically lubricated (external main bearing 2). The shaft carries a flywheel, and the two pinions 1 and 2 that transmit power to the valvetrain and the clutch, respectively. Two connecting rods, arranged in a 90° V configuration, act on the crankpin through hydrodynamically lubricated journal bearings; wrist pins connect the two pistons (not visible in Figure 1) to the respective conrods. The overall mass of the crankshaft equipped with the flywheel and the two pinions is 7.56 kg whereas the conrods and the pistons weight, respectively, 0.32 kg and 0.50 kg each. The crankshaft length is about 300 mm whereas other dimensions are manufacturer’s confidential data and cannot be reported.
(a)
(b)
3. Rigid Multibody Modelling
The first step in the modelling process was the definition of a multirigid-body model of the system. Starting from 3D CAD models of the components, these were imported into the multibody environment and assembled by means of kinematic joints (Figure 2). Great care was taken here to build a single degree-of-freedom (DOF) system: joints and joint primitives were arranged with such purpose, preventing the software to automatically detect, and remove, any redundant constraints.
Model inputs are the combustion forces acting on both pistons, resulting from an experimental test campaign, and the external torques acting on pinions 1 and 2; the former was estimated by the manufacturer by means of some analytic model of the valvetrain while the latter was obtained by simply enforcing periodic speed oscillations over the engine cycle.
Estimates of the loads acting on both the connecting rod and the crankshaft were retrieved as a result of this first dynamic model; comparable results might be obtained by means of analytic models (as those currently used in the design phase by the manufacturer) based on classic formulations, see, for example, [13].
The described model shows some evident limitations. Firstly, loads at two out of four main bearing locations cannot be evaluated, since static indeterminacy of the assembly cannot be treated with a purely kinematic approach: only the central main bearings were thus introduced in the model, neglecting the constraints associated with the external bearings. Secondly, the elastodynamic effects, which might have a major impact on the estimated loads, provided the high rotational speed of the engine (10000 rpm), are completely neglected.
These drawbacks can be removed by introducing nonideal joints at the four main bearing locations, and by modelling component flexibility. Both procedures concur towards the definition of a more realistic model, which should be capable of predicting the interactions between the crankshaft and the engine block with improved accuracy. Only the latter has been considered thus far, and will be presented in the next section.
4. Flexible Multibody Modelling
The development of a flexible multibody model usually requires a preprocessing stage in which a finite element (FE) model of the designated component is produced. Whenever the FE discretization results in a large model order (as in the case of the crankshaft FE model represented in Figure 3), which is often the case when dealing with solid meshes, some reduction procedure is needed to reduce the number of coordinates required to describe component flexibility; this can then be accounted for in nonlinear dynamic simulations at acceptable computational cost. The assumption here is that the deformation of the flexible body keeps small with respect to a body local reference frame, which in turn undergoes large, nonlinear motion relative to an inertial global reference frame; linear elastic theory is used to describe such deformation, which is approximated as a linear combination of a number of shape vectors: where is a vector of physical coordinates, that is, translations (and rotations, if any) of the FE nodes, is a matrix holding the shape vectors as columns, and is a vector of generalized coordinates. Usually, the dynamic behaviour of the FE component in a certain frequency range of interest can be captured using a much smaller number of generalized coordinates compared with the original number of physical coordinates: there lies the model reduction.
Several model reduction methods were proposed in the past [6–8] which rely on a coordinate transformation as the one defined by (1); they basically differ from each other in the selection of the mode set used to build the transformation matrix , which can include normal modes, constraint modes, and attachment modes. A comprehensive overview of those techniques, which are referred to as component mode synthesis (CMS) or substructuring techniques, can be found in [9].
4.1. Craig-Bampton Model Reduction
In this study, the Craig-Bampton (CB) approach has been used and will be briefly reviewed here. The equilibrium equation for the free-free, undamped structure holds: where and are the FE mass and stiffness matrices, respectively, and is the external force vector. The physical DOFs are partitioned into two complementary sets, the interface DOFs (subscript , “active”) and the interior DOFs (subscript , “omitted”): where it has been assumed that external loads, that is, external forces/torques, are applied only at the interface DOFs; clearly, this requires a proper -set definition. The following coordinate transformation is then introduced: in which is a matrix of shapes obtained considering the lower partition of the static portion of (3) and solving for the -set displacements whereas is a matrix of shape vectors satisfying the -set eigenproblem The first columns of the transformation matrix in (4) represent the static deformation shapes of the component when subjected to unit displacements at each of the -set DOFs, the other being restrained; these are termed constraint modes in the literature. The last columns are fixed-interface normal modes, that is, eigenvectors, representing the dynamics of the substructure interior relative to the interface; the corresponding eigenvalues are collected in the diagonal matrix . As one might notice, the generalized coordinate vector comprises both physical displacements, , and modal displacements, ; the fact that interface DOFs are retained in the reduced representation greatly facilitates component coupling, this being probably a major reason for the success of the CB method.
Introducing the transformation in (4) into (3), and premultiplying by the transpose of the transformation matrix, the CB substructure representation is obtained: where the fixed-interface modes have been normalized with respect to the -set mass matrix. In (7), the following positions hold: and being the statically reduced mass and stiffness matrices, respectively. The matrix is a modal participation factor matrix, containing the multiplication factors for the acceleration inputs at the interface DOFs governing the response of the fixed interface modal coordinates:
A modal selection problem arises at this point. Clearly, the two most important aspects of such a problem are model order and model accuracy: an optimal reduction would result in the minimal set of component modes which ensures acceptable accuracy in the simulation results. Concerning the CB mode set, the definition of the constraint modes directly comes from the coordinate partitioning process, which in turn is dependent upon the choice of interface DOFs; usually, these are selected based upon the knowledge of constraint and applied load locations. On the other hand, the choice of fixed-interface normal modes to retain in the reduced representation is somewhat arbitrary. In this work, the normal mode selection is carried out in accordance with a modal ordering scheme based on the effective interface mass (EIM) measure of dynamic importance, which constitutes a generalization of the effective mass concept that has been traditionally used by structural dynamicists to assess the completeness of a reduced mode set. The EIM approach was introduced by Kammer and Triller [11, 12] already in the mid nineties, with main goals being pure structural dynamics and control-structure interaction applications; here, the authors propose a possible extension towards multibody dynamics, which is the main novelty character of this paper.
4.2. Normal Modes Selection Procedure
The EIM measure was developed by considering (10), from which the time-domain response of the th fixed-interface mode to the -set inputs can be computed as: being the mode eigenfrequency. Denoting the convolution integral as , the previous can be expressed in a more compact form: from which the corresponding modal acceleration results: Once again considering (7), its upper partition holds: in which the product has been designated as , representing the portion of the load at the interface due to the response of the fixed interface modes. Using (13), this can be expressed as The norm of the matrix thus gives a relative measure of the contribution of the th fixed interface mode to the loads at the interface; summing all the contributions produces the so-called reduced interior mass (RIM) matrix: Substituting (9) into (16) and bearing in mind that, as mentioned, normal modes are mass normalized, give after some simple manipulation: It is worth noting that such expression can be computed based solely upon the partitioned FE mass and stiffness matrices, and is thus totally independent of any eigenvalue solution. By using some appropriate matrix norm, for example, the trace norm, it can be used as an absolute reference with respect to which the dynamic importance of each mode shape can be computed. The EIM value of the th fixed interface mode is then introduced: By using an analogous expression, a measure of dynamic completeness of a reduced representation in which normal modes are retained is given by the EIM cumulative sum: where is obtained by simply extending the summation in (16) to the -set modes.
A threshold value can be specified for the expression in (19) to set a required level of dynamic completeness for the reduced model. This implicitly defines how many normal modes, and possibly which ones, should be retained in order to achieve the target value, thus completing the CB mode set selection procedure.
It should be noted that though represents a cumulative sum of the -set EIM values, it can be computed independently from the specific criterion used for modal selection, and can thus be used as a general modal completeness indicator whether the EIM approach is adopted or not.
4.3. Orthonormalization of the CB Mode Set
It is worth recalling that the CB mode set obtained by the raw union of the constraint modes and the selected normal modes is not suitable for direct use in a multibody system simulation [14]. For the implementation in a multibody code, an orthonormalization of the CB mode set is generally applied (by solving a further eigenvector problem) thus leading to an orthogonal basis whose modal coordinates are no more directly split into the physical displacements related to the constraint modes and the modal displacements associated with the normal modes. The resulting orthogonal modes (hereinafter, called orthogonal reduced solution) are not eigenvectors of the original system, but they are the eigenvectors of the Craig-Bampton representation of the system. Their physical description is difficult, but, according to the authors’ experience, it can be observed that, in general, the first elastic modes and their corresponding natural frequencies are basically identical to the mode shapes and natural frequencies of the original system analyzed in free-free boundary conditions.
4.4. Application to Crankshaft Modelling
An FE model of the crankshaft was produced first by means of the commercial software PATRAN. The CAD geometry was simplified by removing all those features that do not contribute significantly to the dynamic behaviour of the component, in order to simplify the meshing phase and to limit the total number of elements (Figure 3). Solid elements TET10 were used for the shaft whereas the two pinions and the flywheel were modelled as point elements (MASS2, with properly defined inertial properties retrieved from an accurate CAD model) connected to the mesh nodes by means of rigid spiders (RBE2, not represented in Figure 3).
The discussed model reduction and modal selection methodologies were implemented as following:(i)interface nodes with 6 DOFs each were defined where the loads transmitted to the crankshaft are applied, that is, at the four main bearings and the two crankpin journals, as well as at clamping locations of the two pinions, for a total number of interface DOFs;(ii)a CB solution including 80 normal modes was computed trough NASTRAN FE solver;(iii)the RIM matrix was computed starting from component FE matrices, by implementing an ad hoc NASTRAN DMAP alter;(iv)the EIM values were computed in MATLAB starting from CB and RIM matrices exported from NASTRAN in the previous steps. Slightly different expressions with respect to those in (18) and (19) were used in order to properly weight translational and rotational DOFs; for further details, see [11].
Figure 4 shows EIM values for the computed 80 normal modes (those having an EIM value above 0.001 are highlighted): as one might have expected, in this case the low frequency modes give a major contribution to forces at the substructure interface, exhibiting high values of . However, clearly there can be some high-frequency mode contributing as much as some lower-frequency one (see, e.g., modes 20 and 54, showing similar EIM values).
Figure 5 shows the EIM cumulative sum as a function of the number of retained normal modes, , for modal ordering schemes based upon increasing eigenfrequency and decreasing EIM value, respectively. Curves are very close in this case, both approaching quite quickly the maximum value of 1 (corresponding to the unreduced case ). Of course, for any given number of retained modes, the line related to the modal selection scheme based upon eigenfrequency lies below, or equals, the one related to the EIM-based approach: it is evident that the mode set chosen based upon EIM values constitutes the minimal set providing any specified dynamic completeness expressed in terms of EIM cumulative sum.
In this study, the target normal mode set was defined as the one holding at least 90% of dynamic completeness (), corresponding to 18 normal modes; in this particular case, the selected normal mode set results are independent of the chosen modal ordering scheme, since the 18 lowest frequency eigenmodes are also those exhibiting the highest EIM values, see Figure 4. The mode selection procedure could then be easily implemented by using NASTRAN standard functions to produce the CB solution.
4.5. Nonlinear Dynamic Simulations
At this point, the rigid crankshaft within the multibody model was replaced by its corresponding flexible component. Some multibody models were developed following both the proposed EIM-based procedure and the standard practice consisting of truncating the modal basis of a refined model at a certain cutoff frequency, defined as a multiple of the maximum frequency of interest. In the latter case, the EIM cumulative sum cannot be computed since truncation acts directly onto the orthogonal mode set. Five models are presented in this study:(i)Model A: 80 normal modes are retained in the crankshaft CB solution corresponding to an EIM cumulative sum .(ii)Model B: following the discussed mode selection procedure, only 18 crankshaft normal modes are retained to provide an EIM cumulative sum .(iii)Model C: obtained truncating the orthogonal reduced solution of the crankshaft in Model A at the frequency of 10 kHz.(iv)Model D: the same as Model C but with a cut frequency chosen at 20 kHz.(v)Model E: the same as Model C but with a cut frequency chosen at 30 kHz.
Model A was considered as reference in this study, due to the high dynamic completeness. The reaction forces acting on the four main bearings and those acting on the conrods were calculated by means of dynamic simulations of the described models. An assessment concerning model accuracy was then performed, taking Model A as a reference. Relative deviations were, therefore, evaluated for Models B÷E with respect to Model A over an averaged engine cycle (i.e., two crankshaft full revolutions). In each case, these errors were normalized based upon the peak load value as estimated for Model A, in order to assign a lower importance to deviations occurring at low levels. As an example, Figure 6(a) shows the force magnitude on the central main bearing 2 as evaluated for the aforementioned models; the load value is normalized for data confidentiality. Figure 6(b) shows the relative percentage deviations exhibited from Models B, C, D, and E, respectively. Concerning this specific result, it is evident that Model B performs much better than Models C, D, and E. In particular, due to its capability of capturing higher-frequency dynamic effects, it is able to replicate with extreme accuracy the load trend estimated by the reference model whereas Models C, D, and E exhibit significant deviations. Moreover, the peak load value is perfectly captured by Model B (error of 0.03%) and significantly underestimated by the other three models (errors of 15.8%, 8.6%, and 7.6% for Models C, D, and E, respectively). Similar results hold for the loads acting at other interface locations as well. In particular, for the external main bearing 2 (Figure 7), Model D and E globally overestimate the load values while Model C provides results that can be hardly related to the reference ones; on the other hand, Model B proves to be very accurate. It is worth mentioning that both the peak load values and the load trends (strictly related to the frequency content) are essential input data for a proper design of components subjected to dynamic loads (due to fatigue and durability issues).
(a)
(b)
(a)
(b)
The RMS value of the relative percentage deviations was then taken as a simple and effective metrics to measure the result accuracy of the simplified Models B÷E with respect to the reference Model A. Table 1 reports the values of these parameters as computed for Models B÷E concerning the reaction forces at the crankshaft main bearings and those transmitted by the conrods (for instance, values on the second row are related to Figure 6(b)). It is evident that Model B outperforms Models C, D, and E, resulting in all cases extremely close to Model A. It should be noted that, for all the models, the largest errors are those related to loads acting on the external main bearings, which likely present a greater sensitivity to high-frequency dynamic effects than the other interface locations.
It is worth mentioning that the computational cost of Model A is much higher with respect to the other models simulated here, the reason for that being the large modal basis associated with the crankshaft it incorporates. Such a model, which acts as a reference here due to its nature, might then be too expensive for a practical use in the design process, where the requirement of optimized component design paralleled with the need of meeting tight deadline calls for efficient simulation models. Within the context of structural dynamics and multibody system simulation, it is a common practice to reach this degree of efficiency by truncating the modal basis of the components of interest according to an arbitrary cutoff frequency; this approach was used to obtain the crankshaft models included in Models C, D, and E, which differ from each other solely by the definition of such arbitrary parameter. The aim in putting those three models into operation was here to stress the fact that the choice of this parameter heavily affects the accuracy of the results, which represents another major requirement a simulation model should satisfy. Therefore, a good dose of experience is required to choose a proper cutoff frequency, allowing the analyst to be confident in the results of his simulations; this might not always be the case, especially for components and/or systems the analyst might not be familiar with. Furthermore, no direct indication of effectiveness is provided by such an approach, so the analyst has to rely entirely on his engineering judgment skills. Alternatively, one could develop a number of models having increasing cutoff frequency, in order to determine an “asymptotic trend” of the results so as to define the most appropriate model (Models C, D, and E were developed in this perspective).
The EIM-based approach, adopted here to define the crankshaft reduced solution in Model B, allows overcoming these drawbacks: relying on a quantitative assessment of the dynamic completeness of a reduced representation of the flexible body at issue, it directly provides an indication of effectiveness, referring in particular to the capability in predicting accurate interface loads. Furthermore, the adoption of such a technique complies with the requirement of model efficiency: computational costs are reasonable when compared to those associated with the reference model, and the possible need for the analyst in defining multiple models for accuracy assessment (something which might remind our “refinement” from Model C through D through E) is removed.
5. Conclusions
The development of a flexible multibody model of a motorcycle engine cranktrain has been presented. A rigid-body model was firstly implemented, providing quite similar results as those obtained using some analytical approach by the manufacturer. The main limitations of such model have been discussed, which require the adoption of a refined modelling approach. In particular, the inclusion of a flexible model of the crankshaft has been reviewed in this paper, with an emphasis on the model reduction methodology. A procedure based on EIM values has been employed to select the normal mode set used to construct a CB reduced representation of the component: such a systematic approach, which relies upon the definition of a measure of dynamic completeness of the reduced model, leads to improved load prediction and optimal simulation times, as demonstrated by the comparison with a refined model taken as reference. In fact, in this case some simplified models, employed to assess the degree of accuracy of different model reduction approaches, were compared with a model having a large modal basis, whose results were considered as the reliable reference. It should be pointed out, however, that some validation would also be needed for such reference model, this representing a best practice for every modelling activity; this is outside the scope of the present paper and will be the objective of future investigations.
Further developments involve the integration of bearing models into the multibody system, and the implementation of flexible models for the connecting rods and, mainly, for the engine block; in this case, the described model reduction and modal selection procedures might lead to significant benefits in terms of result accuracy versus computational cost. Nonlinear dynamic simulations of the resulting model will ultimately allow a more accurate constraint reaction prediction, which will constitute the basis for the subsequent stress evaluation and design optimization process.
Acknowledgments
The research has been performed in the framework of the FIRB project no. RBIP068WAA titled “Definition of an Integrated Platform for the Design of Engine Components of Motorvehicles Characterized by a Low Weight/Power Ratio and Reduced Environmental Impact, by Means of Novel Modelling Methods and by Carrying Out Research on Innovative Materials and Process Technologies, Also Transferable to Other Vehicle Components.” The support of the Italian Ministry of Education, University, and Research is gratefully acknowledged.