Abstract
This paper introduces a novel method for ranking and selecting the interior modes to be retained in the Craig-Bampton model reduction, in the case of linear vibrating systems under periodic excitation. The aim of the method is to provide an effective ranking of such modes and hence an optimal sequence according to which the interior modes should be progressively included to achieve a desired accuracy of the reduced-order model at the frequencies of interest, while keeping model dimensions to a minimum. An energy-based ranking (EBR) method is proposed, which exploits analytical coefficients to evaluate the contribution of each interior mode to the forced response of the full-order system. The application of the method to two representative systems is discussed: an ultrasonic horn and a vibratory feeder. The results show that the EBR method provides a very effective ranking of the most important interior modes and that it outperforms other state-of-the-art benchmark techniques.
1. Introduction
The synthesis of accurate finite element (FE) models of complex vibrating mechanisms and structures usually imposes the use of fine meshes, leading to large dimensional models. Unfortunately, large dimensional models are often undesirable since they are difficult to handle and often numerically ill-conditioned because of the presence of large condition numbers (i.e., the ratio between the maximum and the minimum singular values of a matrix). Therefore, they can be too cumbersome for use in simulations, synthesis of controllers [1] (mainly of real-time model-based controllers), model-based design [2], and optimization techniques [3–6]. In order to cope with such an issue and to optimally trade off model accuracy and dimension, several model reduction techniques have been developed in the last decades [7]. In the field of structural dynamics and multibody system dynamics one of the most implemented reduction techniques is the Craig-Bampton (CB) method [8] because of its effectiveness and easiness of implementation.
Two peculiar reasons make the use of the CB method particularly advantageous. First of all, it allows retaining physical displacements in the reduced-order coordinate vector. Physical coordinates may be of interest, for instance, for coupling a system with other systems designed either in the mechanical domain or in other domains (consider not only the CMS [8, 9] but also the coupling between electromechanical systems in multiphysics simulations [10]). Physical coordinates are also of interest when modifications of physical parameters should be computed through inverse structural modification techniques [5]. Secondly, the CB method makes direct use of second-order model formulations (i.e., models with mass and stiffness matrices) rather than first-order ones: second-order models are usually preferred in structural dynamics [7]. These features are not met by many other model reduction techniques, such as, for instance, Krylov subspace and the balanced truncation [9].
Basically, the CB method is a combination of Guyan’s static condensation and modal truncation. Indeed, it uses the static deformation shapes of some degrees of freedom (dofs) of the system (the so-called master dofs), as in Guyan’s condensation, and then enriches this space with a reduced set of coordinates referred to as the interior modes of the system. The practical implementation of the method imposes, as a first step, partitioning the displacement vector into a subset of master dofs , named the external dofs, and a subset of slave dofs , the interior dofs, with , and then transforming it into a set of hybrid coordinates through the CB basis (the transformation matrix):The definition of the transformed hybrid coordinates introduces vector , which includes the aforementioned interior modal coordinates (often referred to as the fixed interface modal coordinates) which replace the slave dof coordinates . In contrast, the master dof coordinates are entirely retained in and are, therefore, usually chosen as those lying at the interfaces that are to be coupled to other systems and/or as those where external forces are applied. In (2), is the nonsingular CB transformation matrix, is a Guyan’s reduction basis, and is the eigenvector matrix obtained by solving the eigenvalue problem of the system with all the external dofs constrained (which will henceforth be referred to as the constrained subsystem). The columns of are the interior normal mode shapes, henceforth referred to as , associated with the interior modal coordinates collected in vector . Finally, and represent, respectively, the identity and the null matrices.
In order to reduce model dimensions, the interior modal coordinate vector is truncated to a smaller vector , , by leading to the reduced coordinate vector , which is expected to provide an effective approximation of :In (3) is the rectangular CB reduction matrix, obtained from (2) by removing the columns of associated with the interior neglected modes.
The crux in the practical implementation of the CB method is performing an optimal selection of the interior modes to be retained in reduced models. Typically, as a widespread rule of thumb, model reduction relies on retaining the interior modes whose eigenfrequencies are not greater than about twice the highest operating frequency [11]. Such a sorting rule based on the eigenfrequencies of the interior modes (henceforth referred to as SBE, for brevity), however, is not based on rigorous principles and neglects the characteristic of the external excitation influencing the system response. As a consequence, it may lead to rough approximations of the full-order model by discarding high-frequency interior modes whose participation in the system dynamics is important, or, in contrast, may lead to large dimensional reduced models including low-frequency modes providing negligible contributions.
In order to improve the effectiveness of the reduction, some methods have been proposed in literature to rank and select the interior modes. These techniques can be seen as “auxiliary” methods for CB reduction, since they operate between the two steps of the CB method. This idea is schematically depicted in Figure 1. Among the ranking methods proposed in literature one should at least recall the “Component Mode Synthesis χ” (CMS χ) [12], the “Effective Interface Mass” (EIM) [13, 14], and the “Optimal Modal Reduction” (OMR) [15], which have been proved to be effective techniques that allow reducing model dimensions while preserving accuracy. Their effectiveness is certainly higher than that of the traditional SBE approach. Basically, these approaches rank the interior modes on the basis of some terms representing the contributions of each mode to the dynamic loads at the interface. They are general purpose methods that can be applied to several applications but can lead to approximate and less effective results in some particular cases. Indeed, they cannot handle requirements on the frequencies at which the reduced model should be accurate or on the external forces in terms of both spatial distribution and frequencies. Such requirements are often known, mainly in those systems which are designed to operate excited at a specific frequency, for instance, in the neighborhood of a resonance frequency, with a prescribed shape and in the presence of a known external force. Typical examples are open-loop resonators, such as ultrasonic horns or vibratory feeders, or systems in the presence of a single-point harmonic force excitation. For these systems, the availability of reduced models is an essential need for performing numerical simulations, such as multiphysics simulation coupling mechanical and electromechanical models or when simulating complex plants (e.g., manufacturing plants), and for developing model-based design and optimization techniques (mainly when inverse structural modification is employed, see, e.g., [4, 5]). Indeed, although the accuracy of a model might decrease when large modifications of some critical physical parameters are made, it has been widely demonstrated that using simplified, though approximate, models increases the numerical reliability and solvability of the design techniques and also allows applying some design methods that do not work, or may fail, for large scale models. For instance, using smaller problems improves the conditioning number of matrices and allows solving least-squares problems, which are often ill-posed by nature [5]. As for the computational effort, reducing models reduces more than linearly the analysis time. Therefore, the computational cost for reducing models is balanced out by the improvement in the application of the techniques, whenever the reduction is capable of preserving model accuracy.
In order to address such an issue, this paper introduces a novel energy-based ranking (EBR) method for the selection of the CB interior modes to be retained in the reduction of dynamic models of vibrating systems under periodic excitation (with either single- or multiharmonic components), whose harmonic content and spatial shape are assumed to be known and constant, at least approximately. The goal is finding an optimal sequence according to which the interior modes should be progressively included to achieve a desired accuracy of the forced response of the reduced-order model at the frequencies of interest, while keeping model dimensions to a minimum. The underlying idea is that the most important interior modes are those providing the largest energy contributions to the system forced response. The interior modes are therefore ranked by using scalar coefficients representing the contribution of each interior mode to the mean kinetic and potential elastic energy stored by the complete system in a period of the external force. These coefficients provide a measure of the importance of each interior mode to the forced response of the full-order system and account explicitly for both the frequencies and the spatial distribution of the force.
The interior modes with the highest values of the coefficient are those to be retained in the reduced models. The analytical formulation of the coefficients ensures the straightforward practical applicability of the method.
Clearly, in the presence of large modifications of either the spatial distribution or the shape of the external force vector, the reduced-order model should be updated to account for the changed contributions of the interior modes. Nonetheless, the proposed method solves a relevant issue, and several vibrating systems meet the conditions it requires.
The paper is organized as follows: the theoretical development of the EBR method is proposed in Section 2, while the EBR assessment is discussed in Section 3. Two different test cases of industrial interest are investigated: an ultrasonic horn, of the kind of those usually employed in ultrasonic welding and modeled through solid FEs and a vibratory linear feeder as those usually adopted in manufacturing or packaging plants for conveying small products or parts. A comparison with all the aforementioned state-of-the-art ranking methods is also proposed to prove the benefits of the theory developed. Concluding remarks are given in Section 4.
2. The EBR Method
2.1. System Energy in Physical Coordinates
Let us consider -dimensional linear time-invariant and undamped vibrating system, represented through its stiffness and mass matrices and the physical coordinate vector . The total system mechanical energy is given by the sum of the elastic and kinetic energy contributions:The system is supposed to be excited on the master dofs by a set of periodic external nodal forces , which are represented as the sum of finite numbers of harmonic components :In (5) is the th harmonic component of the periodic force, and , , are, respectively, the amplitude, the angular frequency, and the relative phase of the th harmonic component exciting the th dof.
By applying the superposition principle, the system steady-state response to such a force is represented as the sum of the responses to each single-harmonic component :where , are, respectively, the amplitude and the relative phase of the response of the th dof to the th harmonic component. On the basis of (6), the system total energy defined in (4) can be rewritten highlighting the contributions at each frequency:The mean value of the time varying energy in a period of excitation τ is, therefore,In the development of the ranking method, the mean energy will be accounted for, since it is a time-independent scalar value providing a meaningful and concise measure of the elastic and kinetic energy stored in the system in a period. Equation (8) can be rewritten in the following more compact form:where is the amplitude vector of the response to the th harmonic component. The notations and denote diagonal matrices having, respectively, the sine and cosine of on the main diagonal.
2.2. Energy Contributions of the Interior Modes
In order to evaluate the contribution of each interior mode to the system energy, and in particular to its mean value , the set of physical coordinates in (9) is transformed in the CB basis through (2):where is the vector of the amplitude of the response to the th harmonic component (in hybrid coordinates), and are diagonal matrices defined as follows:(the notation represents a square diagonal matrix with the elements of an arbitrary vector v on the main diagonal). The terms in (11) represent the relative phase of the response of the th interior modal coordinate to the th harmonic component. Equation (10) also introduces matrices and , which are the stiffness and the mass matrices in the CB basis, computed in accordance with the usual CB theory. If, without any lack of generality, matrices and are partitioned in accordance with the definition of the external () and interior () coordinates,and the modal matrix is normalized with respect to the mass matrix of the interior dof subsystem , then the following expressions are obtained for and [8]:The diagonal matrix collects the squared angular eigenfrequencies of the constrained subsystem. In order to provide a clearer expression of the equations, the submatrices of the matrices in (13) will be hereafter referred to through the following compact notation, with the obvious meaning of the symbols:By making explicit the entries of , , and in (10), and by employing the notation introduced in (14), the contributions to the mean energy provided by the master dofs and by the interior modal coordinates can be split as follows:In (15) three distinct terms can be recognized within brackets. The first term (first line) represents the contribution of the -dimensional external dof subsystem (represented through and ) regardless of the motion of the interior modes. Therefore, it can be discarded in the evaluation of the contribution of each interior mode to . Conversely, all the other terms within brackets explicitly depend on the interior modes and therefore must be accounted for. In particular, the second term (second line) represents the inertial coupling between the master dofs and the interior modal coordinates, while the third term (third line) only depends on the interior modes. The second and third terms within brackets of all the th harmonic components will be henceforth collected in the scalar variable named , which can be also expressed as the summation of the contributions of each interior mode (indexed through ), over the harmonic components of the force (indexed through ):The amplitudes and in hybrid coordinates of the forced responses to each harmonic component in (16) can be rewritten as a function of the external force, by means of receptance matrices:where is the amplitude vector of the th harmonic component of the force vector (see (5)) and and are receptance matrices. By substituting (17) in (16), can be explicitly written as a function of the external force: Such an expression clearly shows that the contribution of the th interior mode to the system mean energy can be evaluated through the scalar coefficients :In the equation above, is the th row of the receptance matrix , relating the response of the th interior mode to the th harmonic component of the force, .
Clearly, the larger the value of is, the more the th interior mode contributes to the system response in the presence of the periodic force defined through (5). It is here therefore proposed to rank the interior modes in descending order based on the values of . Then, they can be progressively included in the reduced model until a desired model accuracy is achieved.
3. EBR Method Application and Assessment
This section proposes two different test cases for demonstrating the EBR method effectiveness: an ultrasonic sonotrode (Section 3.2) and a vibratory feeder (Section 3.3). Both systems are often employed in industry or research laboratories and are designed to generate suitable vibrations excited by periodic forces. Such forces usually have a few harmonic components and their spatial distribution is known. They are, therefore, well-suited for the application of the proposed method. On the other hand, the availability of reduced-order models is essential for model-based design or optimization of these devices, both whenever direct approaches are adopted (i.e., through extensive simulations) and when inverse dynamic structural modification techniques are employed [5].
3.1. Criteria for Result Evaluation
A reduced-order model should provide an accurate description of the system forced response to the periodic force, in terms of both spatial distribution and amplitude. Two parameters are adopted to evaluate the correctness of the approximation provided by the reduced models and hence the method capability to select the most important modes to be retained. The first one is the modal assurance criterion (MAC) between the vectors of the forced responses in all the FE model nodes (at an arbitrary time ) computed through the full-order model () and through the reduced-order ones (): In practice, the MAC is the squared cosine between the two mentioned vectors and should hence approach 1 to ensure identical spatial distributions of the forced response (i.e., parallel vectors). In the case of the reduced models, whose dimension is smaller than , the -dimensional vector is computed by mapping the reduced set of hybrid coordinates into an approximated set of physical coordinates by means of (3), where vector includes the interior modal coordinates retained in each reduced model evaluated.
Since the MAC does not provide any information on the amplitude of the system response, and therefore of the receptances, the relative gain error is introduced and defined as follows:Such a parameter is the relative percentage error between the forced response computed through the full-order model and through the reduced-order ones , evaluated at an arbitrary time for one representative, or sample, master dof (denoted through index ). Clearly, since the same force is considered for both models, can be also seen as the relative percentage error of the receptances.
Both MAC and obtained by the proposed method have been compared with those obtained by adopting the other methods available in literature and quoted in the Introduction. Although these methods are not specifically developed for these kinds of applications, they are the most important techniques for ranking the interior modes available to date, and their effectiveness has been extensively proved.
3.2. First Test Case: Ultrasonic Horn
3.2.1. Device Description and Experimental Model
The theory proposed is firstly applied to the device shown in Figure 2(a), which is an ultrasonic horn (or sonotrode). Sonotrodes are, for example, employed in welding of plastics and nonferrous metals, cleaning, and cutting. A sonotrode is one of the components of the so-called ultrasonic stack, consisting of a piezoelectric transducer, a booster, and a sonotrode. Piezoelectric transducers are the actuators transforming electrical energy into high-frequency mechanical vibrations, boosters amplify the amplitude of vibrations, and finally, sonotrodes, brought into contact with the workpieces, transfer mechanical vibrational energy from piezoelectric transducers to workpieces. All the components of the stack are tuned to resonate at the same ultrasonic frequency. For example, in plastic welding, application for which the experimental test-bed studied is designed, the ultrasonic stack is tuned to the eigenfrequency of sonotrode first longitudinal mode [16]. Therefore, transducer must be driven at a known frequency which should match or be as close as possible to such an eigenfrequency.
(a)
(b)
In order to correctly model the sonotrode studied, some experimental measurements have been carried out, which have led to the identification of the correct properties of the sonotrode material (titanium alloy Ta6V): mass density 4340 kg/m3, Young’s modulus 115 GPa, and Poisson’s ratio 0.32. Measurements have also corroborated the typical assumption that these systems are almost undamped. In particular, this assumption is confirmed by the experimental frequency response functions (FRFs) measured in the neighborhood of the sonotrode first longitudinal eigenfrequency (19886.0 Hz). An example of FRF is shown in Figure 3. This FRF represents the ratio between the velocity of one of the sample nodes at the sonotrode tip and the voltage exciting the piezoelectric transducer. A modal damping ratio has been identified through the half-power method. Similar results (not shown here for brevity) have been obtained for the other nodes along the tip. The experimental setup adopted to estimate the experimental FRFs is shown in Figure 4 and is composed by the following:(i)an Agilent 33220A waveform generator, used to drive the transducer through sine sweep voltage signals in the frequency range 19 kHz to 21 kHz;(ii)a Polytec laser Doppler vibrometer (LDV), used to measure the vibration velocities of the nodes at the sonotrode tip;(iii)a LMS SCADAS SCR02 acquisition system, interfaced to a PC running the proprietary software LMS Test and lab for performing the experimental modal analysis.
A model of the sonotrode has been obtained through FEs by employing solid tetrahedral elements with eight nodes and three dofs (i.e., Cartesian coordinates) per node. In order to achieve an adequate accuracy in the representation of the system dynamics and to match the experimental response, it has been necessary to mesh the model very finely (see Figure 2(b)). The resulting FE model has 8685 dofs, and the predicted eigenfrequency of the first longitudinal mode is 19887.71 Hz. Clearly, handling such a large dimensional model with mass and stiffness matrices of dimensions equal to the number of dofs is cumbersome, not only for their size but also for the large condition number, which causes large ill conditioning whenever numerical algorithms for model-based design are to be applied (see, e.g., [5]).
Following the CB approach, a set of 21 Cartesian coordinates has been chosen as the master dofs, while the remaining 8664 coordinates are the slave dofs. The set of master dofs includes the displacements of the nodes lying at the physical interface between the sonotrode and the booster, through which the mechanical vibrations are transmitted from the booster to the sonotrode. In particular, the force applied to such master dofs is a single-harmonic force whose spatial distribution along the longitudinal direction is uniform in all the nodes, while it is zero in the other directions. The force frequency is set almost equal to the natural eigenfrequency of the sonotrode first longitudinal mode ( is not set equal to it, to avoid trivial numerical problems arising from exciting an undamped system at a natural frequency). The sonotrode is designed to respond to that force with uniform displacements along the tip to ensure regular welding.
3.2.2. Method Application and Results
The application of the EBR method leads to the results summarized in the logarithmic plot shown in Figure 5, where the two evaluation parameters defined in Section 3.1 are plotted as functions of the dimensions of the reduced-order models. In the same figure the evaluation parameters obtained by adopting the other ranking methods available in literature are plotted too. These models with increasing dimensions have been obtained by adding the interior modes progressively, following the ranking order provided by the methods investigated.
(a)
(b)
The capability of the EBR method to ensure accuracy through a minimum set of properly selected interior modes is clearly proved by the results obtained. Indeed, the convergence of the EBR method to the ideal results (i.e., and ) outperforms the ones of the benchmark methods, which, in turn, are more effective than the empirical sorting based on eigenfrequency.
Figure 5 shows that error levels are not monotonically decreasing with respect to the number of interior normal modes. Although the common sense would suggest that the model accuracy improves each time an additional interior mode is included in the model, other papers on model reduction have already shown results with nonmonotonic behaviors (see e.g., [17]), especially when the size of the reduced model is much smaller than the one of the full model. Clearly, these nonmonotonic behaviors should be regarded only as “local” effects. In contrast, the overall trend confirms the expected accuracy improvement when dimension increases.
Table 1 collects some more results. For instance it proves that, by adopting reasonable accuracy thresholds for the MAC and of, respectively, 0.999 and 1%, a 775-dimensional model (with 21 master dofs and just 754 out of the 8664 interior modes) turns out to be adequate if the retained modes are selected through the EBR method. This leads to a percentage order reduction ratio (i.e., the ratio between the number of dofs neglected in the reduced model and the number of dofs of the full-order model) of 91.1%. Conversely, significantly higher model dimensions are needed to achieve the same accuracy through all the other benchmark methods. Table 1 also states the percentage relative errors in the estimation of the eigenfrequency of the first longitudinal mode, with respect to the frequency provided by the full-order FE model (i.e., 19887.71 Hz), denoted as . This information is provided just for such a mode since it is, in practice, the sole mode of interest for the forced response. The eigenfrequencies have been computed for the proposed method and for the benchmark ones, by assuming the model dimensions that allow meeting the accuracy thresholds previously stated (as reported in the same table). This result proves that an almost exact representation of the eigenfrequency of the mode of interest is provided by the proposed ranking techniques with the reduced set of interior modes assumed.
A further proof of the effectiveness of reducing the order of a dynamic model through the CB reduction technique together with the EBR method comes from a comparison with a FE model directly synthesized with a reduced set of dofs. Indeed, small size models might be also obtained by coarsening the mesh at the FE modeling stage. However, this approach usually results in a significant reduction of the model accuracy. For example, if the investigated ultrasonic sonotrode is modeled with a coarse but rather uniform mesh resulting in a 1383 dof model (notice that these dofs are about twice the ones of the reduced model previously obtained by the EBR method), the longitudinal natural eigenfrequency computed through the mass and stiffness matrices of such a FE model is 20541 Hz. This result is affected by a considerable error of 3.3% in the estimation of the natural frequency of the vibrational mode of interest, whose actual value identified experimentally is 19886 Hz. In contrast, the eigenfrequency of the 8685-dimensional model with fine mesh is 19887.71 Hz, while the 775-dimensional model obtained through the proposed EBR method is 19887.74 Hz. The latter value leads to a negligible 0.01% error in the estimation of the actual experimental frequency.
3.3. Second Test Case: Vibratory Feeder
3.3.1. System Description
The second vibrating system to which the EBR method is applied is the one sketched in Figure 6. It represents a linear vibratory feeder of the type usually employed in packaging or manufacturing plants for conveying small components or products. The conveyed products move along the upper beam (the so-called tray), forced by three concentrated electromagnetic exciters, modeled as three independent lumped masses connected to the beam through three linear springs. Six linear springs also connect the tray to a lower beam, which is a support beam, connected to the rigid frame by means of two elastic supports modeled as linear springs. Both of the beams are modeled through a suitable number of four-dof Euler-Bernoulli beam finite elements (see Figure 6). The resulting model has 39 dofs. The system is supposed to be excited by three in-phase forces, having the same amplitude for all the six actuated dofs, which are chosen as the master dofs. Opposite directions are instead assumed for the forces exciting the three master dofs in the upper beam and for those exciting the three lumped masses. Indeed, the electromagnetic exciters are usually driven by identical and in-phase periodic currents, so as to generate the proper motion of the tray and of the conveyed products. Mode ranking and model reduction have been carried out focusing on the two harmonic components of the excitation force at 50 Hz and 100 Hz, which are the first two harmonic components due to the electromagnetic actuators forcing the system and, hence, the two harmonic components with respect to which maximum system response accuracy is required.
3.3.2. Method Application and Results
The results of the investigation are summarized in Figure 7, where the two evaluation parameters are shown as functions of the reduced-order model dimensions.
(a)
(b)
As in the first test case, the EBR method outperforms the other benchmark techniques: a significantly smaller number of interior modes are needed to get a very accurate representation of the system forced response. Once again, Table 2 summarizes the number of interior modes required by the different ranking strategies to meet the accuracy requirement of 0.999 for the MAC and of 1% for . This leads to a model order reduction ratio of 69.2%. The most effective benchmark method (i.e., the CMS χ) requires a number of modes that is 171% higher than the number of modes required by the EBR method.
4. Conclusions
This paper has introduced a new energy-based ranking method, named EBR, developed for performing effective model reduction through the Craig-Bampton method. The EBR method identifies the most important interior modes to be retained in a reduced-order model and is suitable for synthesizing reduced-order models providing a reliable description of the forced response of a linear vibrating system in the presence of a periodic force whose frequency content and spatial distribution are known, as it is often the case in vibrating systems of industrial interest.
The EBR method ranks the interior modes on the basis of scalar analytical coefficients representing the contribution of each interior mode to the mean mechanical energy stored by the full-order system in a period of the external force. Therefore, the EBR method provides the optimal sequence according to which the interior modes should be progressively included in a reduced-order model to achieve a desired level of accuracy at the frequencies of interest, while keeping model dimensions to a minimum.
The theory presented here is general and has been applied to two different devices where model reduction is often of interest for model-based design and simulation: an ultrasonic horn and a vibratory feeder. It has been shown that the EBR method ensures very effective model reduction. The method effectiveness has been also corroborated through the comparison with other “general purpose” ranking techniques available in literature. It has been proved that, in both cases investigated, the EBR method achieves the desired levels of accuracy with a minimal set of interior modes, always considerably smaller than the ones of the other benchmark methods.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors acknowledge the financial support from Fondazione Studi Universitari di Vicenza (FSU) and Fondazione Cariverona, within the frame of the project “Tre Poli.”