A generalized constraint reduction method for reduced order MBS models

In this paper we deal with the problem of ill-conditioned reduced order models in the context of redundant formulated nonlinear multibody system dynamics. Proper Orthogonal Decomposition is applied to reduce the physical coordinates, resulting in an overdetermined system. As the original set of algebraic constraint equations becomes, at least partially, redundant, we propose a generalized constraint reduction method, based on the ideas of Principal Component Analysis, to identify a unique and well-conditioned set of reduced constraint equations. Finally, a combination of reduced physical coordinates and reduced constraint coordinates are applied to one purely rigid and one partly flexible largescale model, pointing out method strengths but also applicability limitations.


Introduction
Multibody system simulation (MBS) nowadays is an important part in various industrial fields like the automotive industry, aerospace, training simulators, and several more. Within the automotive industry, the main areas of application are load data acquisition for further investigations like finite element computations, component design or the identification of production-related uncertainties like inertia parameters, the body mass, the exact location of the center of mass, etc. As many of those applications are of a repetitive nature, it is desirable to optimize such recurring simulation tasks.
In this sense, model order reduction of the underlying system of differential algebraic equations (DAEs), created by the MBS pre-processor is one potential technique. We point B D. Stadlmayr daniel.stadlmayr@fh-wels.at 1 out that the present work covers multibody system models, created by automated modeling strategies using a redundant set of coordinates. Further, such models are processed by solvers using the Hilber-Hughes-Taylor [9] (HHT) algorithm, as it is the case for commercial software like MSC Adams [18] or the freeware tool FreeDyn [7]. Suitable model order reduction techniques, like Proper Orthogonal Decomposition (POD) [2,14,28] or Smooth Orthogonal Decomposition [4] cannot be applied a priori, and need at least one full system simulation. Hence, model order reduction is only meaningful in terms of repetitive simulations. In the past, both methods were applied to dynamical systems by several authors, and the interested reader is referred to [3,5,6,10,11,13,17] and references therein.
In [17], POD model order reduction was recently applied to a suspension system. Therein, the kinematics and dynamics are treated with special attention and hence, the governed equations of motion are very distinct from those treated in this work. We recently compared POD and SOD [26] and further proposed a physical and constraint coordinates reduction method [27] for MBS models. Therein, we focused on rather good-natured models, subject to idealistic boundary conditions. The present contribution demonstrates the applicability and difficulties concerning redundant constraint equations, when dealing with practical, rigid, and flexible models, subject to superimposed vibrations or resonance effects. To overcome the issue of redundant constraint equations in terms of reduced order models, we propose a generalization of the constraint reduction method introduced in [27], based on a singular value decomposition of the reduced constraint Jacobian of the system under consideration. The use of the singular value decomposition method in terms of the elimination of Lagrange multipliers or to handle redundant constraint equations was discussed in the past by various authors; see, for instance, [8,12,15,22,30,31].
In the first part of this work, the system of DAEs under consideration is presented. Following [27], the second part of this contribution introduces the model order reduction technique, and presents a generalization of the constraint reduction method, which improves the conditioning of the reduced order model. In the third part, the reduction process is first applied to an oscillatory supported, rigid crank drive, which introduces superimposed vibrations into the model. The second example shows a flexible offroad vehicle suspension, which is loaded by a frequency sweep at its wheel hub. Finally, the model order reduction method is discussed critically, pointing out the strengths, but also known limitations to the method's applicability.

Redundant multibody systems
Available MBS software, e.g., MSC Adams [18], FreeDyn [7], etc., bear much of the burden of modeling multibody systems. Unfortunately, such universal modeling strategies are not able to identify a minimum set of coordinates representation and, hence, are based on a redundant set of coordinates description. In our case, the six degrees of freedom (DOFs), of a rigid body motion in 3D space, are described by three translational and four rotational coordinates, utilizing Euler Parameters [21]. Body interactions (boundary conditions, joints, etc.) are represented by constraint equations which, in combination with the redundant set of coordinates, lead to the following redundant, nonlinear, index-three differential algebraic system of second order: Let q(t) ∈ R n be the vector of (in the rigid case) n = 7 · k generalized coordinates with k > 0 bodies and its time derivativesq(t) andq(t). Furthermore, let M(q(t)) ∈ R n×n be the global mass matrix, and Q(q,q, t) ∈ R n be the vector of generalized and gyroscopic forces. C(q) ∈ R m represents the vector of constraint equations, C q ∈ R m×n is the constraint Jacobian, and λ(t) ∈ R m is the vector of Lagrange multipliers. For better readability, time dependencies will be omitted throughout the paper. Note that due to the time and state dependency of the system matrices, they must be updated repetitively throughout the simulation.
The MBS model is generated using our in-house MBS software package FreeDyn [7] and both time integration and model order reduction are carried out in Scilab 5.5.1 [23], utilizing a Hilber-Hughes-Taylor (HHT) solver algorithm, which is able to directly handle the present system of DAEs. For further insight into the topic of multibody system modeling, the interested reader is referred to [24].

The reduced order model
A redundant MBS model, as introduced in Sect. 2, often consists of too many coordinates and/or constraint equations. A simple example would be a 2D mechanism modeled in 3D space. Hence, the MBS solver has to deal with a larger set of equations than actually necessary. To reduce this set of equations, a POD model order reduction method, in combination with a modified Galerkin projection and a constraint reduction method following [27], is applied.

A brief review on the Singular Value Decomposition
The model order reduction methods applied to both the physical coordinate and the constraint coordinate reduction are based on the Singular Value Decomposition (SVD) of data, from which we derive basis vectors characterizing the reduction subspaces to project onto. Hence, we briefly summarize the SVDs main features as they are needed for the model order reduction method.
First, recall that the rank of any matrix Y ∈ R m×n is calculated by the SVD where the rank of Y is given by the number of nonzero entries in the diagonal matrix B ∈ R m×n , which is the row/column size of the quadratic matrix D ∈ R l×l . The matrix of left singular vectors L ∈ R m×m spans the column-space of Y and can be split into independent (superscript i ) and dependent (superscript d ) left singular vectors. The matrix of right singular vectors R T ∈ R n×n spans the row-space and can be split accordingly. With no loss of precision, see, e.g., [1,28], Eq. (2) can be rewritten in terms of the independent components as where the independent column space L i ∈ R m×l spans a subspace solely consisting of the unique information from Y.
Beside their usefulness to calculate the matrix rank, the left singular vectors, together with the related singular values, give information about the unique content of the entry matrix. As the singular values are collected in descending order, the largest (first) singular value relates to the first left singular vector, which gives the highest representation of the entry data information. The second singular value relates to the second left singular vector, which gives the highest representation of the set of remaining left singular vectors, and so on. These features can be used to generate a low-rank approximation of the information stored in the entry matrix by a set of left singular vectors.

Physical coordinate reduction based on POD
The physical coordinate reduction method was presented in [27], and it is used in this work with no changes.
The Proper Orthogonal Decomposition can be formulated in terms of the SVD of a given snapshot matrix X. In that sense, and referring to the preceding subsection, POD investigates a high-dimensional data time-series in order to identify a low-dimensional description, which covers the essential characteristics of the original data. We sample the data-time se-riesẊ ∈ R n×h , with n = n trans + n rot + n flex , with the number of time points h > n, on the velocity level. Note that this time-series must be collected in an initial simulation of the original system in Eq. (1). Afterwards, we distinguish between translational, rotational and, if present, flexible coordinate dataẊ trans ∈ R ntrans×h ,Ẋ rot ,Ẋ flex which are processed separately by the SVDẊ In terms of POD, the left singular vectors Φ trans = [ϕ 1 , . . . , ϕ ntrans ], with Φ trans ∈ R ntrans×ntrans are called Proper Orthogonal Modes (POMs) and are uniquely related to the singular values Σ trans ∈ R ntrans×h or Proper Orthogonal Values (POVs). POVs are sorted in descending order in the diagonal matrix Σ trans = diag[σ 1 , . . . , σ h ], Σ rot , Σ flex with their magnitude representing the importance of the corresponding POM in terms of the data signal energy. In terms of the POD model order reduction applied in this work, the right singular vectors W are not used.
The subspace of the reduced order model is found by investigating the POVs decay. Unfortunately, there is no general algorithm or criteria available, which allows to automatically identify the essential number of POVs/POMs. Actually, one has to get a feeling of POVs and their descending behavior. We experienced good results by first normalizing the POVs vector and then searching for drops in POV magnitude by several powers of ten. However, we can give no clear statement on the exact drop-size, which characterizes a sufficient drop. In other cases, if there is no clear drop detectable, it can be helpful to define a maximum error in signal energy and sum up the POVs "from back to front" until the maximum error is reached. The remaining first POVs then characterize the reduced order model's dimension. Again, the maximum error in signal energy is no hard fact but model-dependent. It might also be the case that one has to try out various sets of POMs to achieve the best reduced order model. Further, as the POVs are not only model dependent but depend also on the simulation inputs, the set of POVs and POMs is not necessarily unique. This fact must be kept in mind, if changes to the model or the inputs are made without evaluating the reduction subspace. Therefore, we will deal with the topic of subspace-robustness in the context of parameter identification in the future.
Once a sufficient number of POVs is found (e.g., the first five POVs, r trans = 5), the POM matrix, e.g., Φ trans ∈ R n×n is truncated after its fifth column and the reduced subspace Φ r,trans ∈ R n×rtrans , to project the system onto, is defined by these first five POMs. The model order reduction approach is further applied in analogous manner to the rotational and flexible data, finally resulting in the reduced subspaces Φ r,trans , Φ r,rot , Φ r,flex . Note that the dimension of each subspace is unique and not necessarily equal to that of the others. Finally, the reduction subspaces are then combined in the global projection matrix U ∈ R n×r with r = r trans + r rot + r flex To achieve a reduced order model of Eq. (1), the global projection matrix is used in the modified linear Galerkin projection to project the original coordinates q ∈ R n to reduced order coordinates q r ∈ R r . Note that due to truncating the POM matrix, the projection onto the subspace and back is not exact but only an approximation with an error related to the magnitude of truncated signal energy in the POVs. For this reason, by projecting the reduced coordinates back into the full space can only deliver an approximationq ∈ R n of the original, full coordinates q.
In order to ensure that initial conditions are fulfilled, the linear Galerkin projection is extended by the residual term R q 0 = q 0 − UU T · q 0 , with R q 0 ∈ R n and the original vector of initial conditions q 0 ∈ R n . If needed, the projection can be applied to velocity coordinates accordingly, using the related residual term Rq 0 =q 0 − UU T ·q 0 and hencė q ≈q = U ·q r + Rq 0 .
As the state dependent system matrices are defined in terms of the original system, their evaluation throughout the simulation is a bottleneck of the present method. To update the system matrices, the reduced coordinates q r are projected back into the approximated full space. The system matrices are then updated, using these approximated coordinates. To highlight the fact of an approximation, the system matrices are therefore also denoted with a tilde (M,C, etc.). Finally, the approximated system matrices are projected into the reduced space, to enter the solver algorithm. To summarize, by applying the Galerkin projection and pre-multiplying the ODE part of Eq. (1) with U T , the projected MBS model At this point, we briefly summarize our reasons, originally presented in [27], to apply POD for each coordinate-type separately and on the velocity level: Zero-valued POVs describe possible simplifications in the orientation of a body-fixed coordinate frame, as motion along the corresponding POM is zero. In other words, the body coordinate frame orientation is marginal in terms of the acting motion if the corresponding POV is zero. Hence, a more significant frame orientation can be indicated by POVs, as, for instance, a 3D motion along an arbitrary vector in space may be reduced to a onedimensional motion along a new directional vector, which is described by the corresponding POM. For this reason, we decide to handle translational, rotational, and flexible DOFs separately.
Further, while POD is frequently applied to position data, we use velocity data as input to the POD process. Velocity coordinate data contains more "dynamic information" of the system under consideration. Especially when investigating models which combine small (fast) fluctuations with large (slow) body motions, position data may overvalue large movements, which was, at least to us, unascertainable using velocity data. For a more detailed discussion on our choice of the POD input-data, and the reduction method in general, the interested reader is referred to [27].

Constraint reduction
The MBS model described in this work may introduce unstressed constraint equations, as in the case of a planar problem modeled in 3D space. Furthermore, equal coordinate motion could be forced by appropriate constraint equations leading to redundancy in the time-series of the related coordinates. Such equal coordinate time-series are detectable by the POD process, resulting in the global projection matrix U. Such constraint equations are negligible, as they are fulfilled by the Galerkin projection, which itself acts as a kind of constraint, see [27]. Most importantly, if the vector of generalized coordinates q ∈ R n is reduced to q r ∈ R r , without reducing the vector of constraint equations C ∈ R m , the resulting DAE system might become overdetermined and ill-conditioned, if m > r. This is due to the solver algorithm, which basically has to handle a system Jacobian (see [20]) of the form which can become singular if the physical coordinate reduction results in m > r.
To overcome this issue, we proposed a constraint reduction procedure in [27], which uses the projection of the constraint Jacobian onto the subspace spanned by U. Therein, if U T · C T q equals zero for all instances of time, a smaller set of constraint equations C l (q) ∈ R l with l < m can be defined. However, this method tends to fail if the constraint equations are defined in any direction other than the global axis directions. So, in terms of a more complex model, which consists of boundary conditions in various directions in space, hardly any constraint equations might be identified as reducible by the algorithm proposed in [27].
In order to overcome this issue, the constraint reduction method is generalized in the following section, ensuring that the reduced MBS system is determined and well-conditioned.

A generalized constraint reduction method
The generalization of the constraint reduction method presented in this section, is closely related to Principal Component Analysis [16,25,29]. It is based on the singular value decomposition of the projected constraint Jacobian C q,r = C q U and uses the notation of Sect. 3.1. Recalling the property of the singular value decomposition, to precisely describe the information content of an entry matrix by the independent part of the output matrices and further pre-multiplying Y with (L i ) T , the original redundant data is re-expressed as a linear combination of its independent basis vectors L i ∈ R m×l as with Y unique ∈ R l×m solely containing the essential (unique) information from Y.
Switching to the present problem of constraint redundancy, the method is applied to the constraint Jacobian C q,r ∈ R m×r in an offline phase of the simulation procedure. Therein, m is again the number of original constraint equations and r is the number of reduced coordinates. By setting Y = C q,r and proceeding as proposed, the therewith related independent right singular vectors L i ∈ R m×l and singular values D ∈ R l×l are generated. The constraint Jacobian is then re-expressed as a linear combination of its independent basis vectors C lq,r = L i T C q,r .
As a result, C lq,r ∈ R l×r , l ≤ r, solely consists of the non-redundant constraint information. Note that, in the case of l = r, the reduced MBS model is fully constrained, as the number of remaining constraints equals the number of reduced coordinates. In such a case, the singular values in D should be investigated in a similar way to the physical coordinates reduction in Sect. 3.2, to check whether a drop in the singular value decay indicates a smaller set of constraint equations. Furthermore, note that any constraint equation evaluating to zero in the sense of the projection U T · C T q in [27] is equivalent to a redundant column in C T q , and hence, results in a zero singular value in B, see Eq. (2). Therefore, the constraint reduction method proposed in [27] is a special case of the present method, which possibly allows a subset of the redundant constraints to be identified.
To ensure consistency, the proposed projection is applied to the vector of constraint equations C as It should also be noted that if the algebraic constraint equations (C = 0) hold true for the original system, any linear combination of these terms in the sense of Eq. (13) meets the zero-condition C l = 0. Finally, the physical and constraint coordinate reduced model is defined as Note that in Eq. (14), one could set U T C T lq,r λ l = U T C T q,r L i λ l and argue that the last part is similar to a Galerkin Projection of the vector of Lagrangian multipliers λ ≈λ = L i λ l . In fact, this term arises due to the consequent application of the constraint reduction method and is not intended as a Galerkin projection. Nevertheless, the vector of Lagrangian multipliers has to change its dimension according to the remaining number of constraint equations C l .

Comments with respect to the constraint forces in the reduced model
The model order reduction to physical and constraint coordinates, as presented in this section, allows us to find a reduced order model, which reproduces the original systems coordinate movement with high consistency. Nevertheless, the advantage of finding a smaller set of coordinates by a linear combination of the original coordinates, goes hand in hand with a loss of system information, as we truncate the POM matrix.
To ensure correct constraint forces is still an open field of research and several authors (see [8,12,30,31]) deal with the problem of redundant constraint equations, whether or not they result from singular configuration, overconstrained modeling (as in the case of a Fig. 1 One mass oscillator door being modeled using both hinges), etc. Also in the present case of the reduced order multibody system, the constraint forces turn out to differ from the constraint forces of the original system. This is a major limitation, as the reduced order model cannot be used to generate, for instance, load data in joints, as they might be needed for special FEM simulations. Nevertheless, we can show that the change in the constraint forces is already caused by the physical coordinate reduction (the reduction of q) and not introduced due to the constraint reduction (the reduction of C, etc.).
As a matter of fact, if any original coordinate is weighted with factor zero in the reduced order model, which is a zero entry in all related entries of the projection matrix U, forces acting on this coordinate no longer enter the reduced order model. Therefore, the corresponding constraint force UC T q λ equals to zero. This issue can be demonstrated on the following example of a single mass oscillator ( Fig. 1) with mass m, which is coupled to ground by a linear spring with spring stiffness k. It is actuated by an external force F in the (x, y) direction, acting in the mass's center of gravity. Let the original model consist of two degrees of freedom, which are the (x, y) movements of the cart. The model is subject to one constraint equation, which is C 1 = y = 0. Hence, the mass is forced to move in the x direction only. The constraint Jacobian of this example is of the form The velocity result fileẊ ∈ R 2×h , with the number of time points h, of this example is of the formẊ 16) and hence, the SVD of this data matrix delivers one nonzero and one zero-valued POV. The projection matrix, out of the POM corresponding to the nonzero POV is of the form and transforms the system to a one-DOF reduced order model, as all zero-valued POVs are omitted. Due to this choice, which is also the minimal set of coordinates representation, the constraint Jacobian is projected onto the reduction subspace as As the reduced order model consists of the acting x DOF only, the constraint force U T C T q λ, which is purely depending on the y DOF equals to zero for all instances of time. Hence, Fig. 2 Numerical example-V8 crank drive when transforming the reduced system solution back into the full space, the movement of the single mass oscillator is equal to that of the original system, but the constraint forces equal to zero.
With these limitations in mind, we point out that this model reduction method is not usable for simulations concerning constraint forces as important results.
Beside this limitation, the presented example further points out the necessity to adapt the constraint equations according to the method proposed in Sect. 3.3.1. Without adapting the vector of constraint equations and the constraint Jacobian of the present example, the system Jacobian in Eq. (9) of the reduced system is of the form of and hence singular. Applying the constraint reduction method from Sect. 3.3.1, the SVD of the reduced constraint Jacobian delivers zero singular values only. Therefore, the constraints and the constraint Jacobian are completely eliminated from the solver algorithm and the reduced model is again well-defined and well-conditioned as the system Jacobian reduces to

Numerical examples
In this section, the physical coordinate reduction, together with the newly proposed generalized constraint reduction method, is applied to both a vibrating rigid V8 crank drive (Fig. 2) and a flexible model, based on the front axis suspension of an offroad vehicle (Fig. 3). Both examples are set up in our in-house MBS software tool FreeDyn [7] and processed in the open-source software Scilab 5.5.1 [23].

Example: rigid V8 crank drive under superimposed vibrations
The purely rigid example is based on the V8 crank drive model introduced in [27]. In the present work, the crank drive, Fig. 2, is extended by an engine housing and no longer fixed to the ground, but supported by engine and gearbox mounts, as well as one torque roll restrictor. The engine supports are modeled by virtual spring and damper elements, indicated by three coordinate systems in Fig. 2. Spring and damping parameters of the supports are chosen such that the static sinking is less than five millimeters and the rotation in the cranks longitudinal axis is less than five degrees due to maximum rotation speed. Due to these limitations, Due to superimposed vibrations, the constraint reduction method in [27] fails to identify any reducible constraint equations. This is related to the fact that the constraints are no longer acting along the ideal (global) axis system, but tumble in space. In contrast, the herein introduced generalized constraint reduction method identifies a total of 72 constraint equations as redundant or insignificant. According to Fig. 4(c), the reduced order model consists of 46 constraint equations. Although the reduced order model with 26 translational, 25 rotational, and 46 constraint coordinates is significantly smaller than the original model, it only saves about 5 % of the overall simulation time.
The choice of reduced translational/rotational coordinates was made after testing various reduced coordinate combinations, as the presence of several notable drops ruled out the possibility of defining a set of reduced coordinates in a straightforward manner. Further reduced order models, consisting of (a) 18 translational, 22 rotational, and 35 constraint coordinates or (b) 14 translational, 25 rotational, and 34 constraint coordinates, are compared in Fig. 5. The relative error of various coordinates throughout the model does not exceed 2 % and stays below 1 % bound for most of the simulation time. The overall simulation time is reduced by almost 28 %, using the reduced order model consisting of 14 translational, 25 rotational, and 34 constraint coordinates. Compared to the results in [27], where the fixed-to-ground model was reducible to four translational and six rotational coordinates, it becomes obvious that superimposed vibrations lower the model reduction capability. Nevertheless, an appreciable time saving is still possible.

Example: flexible front axis suspension
The flexible suspension model, see Fig. 3(a), consists of the flexible lower wishbone, the flexible pushrod and the rigid components tie-rod, top wishbone, steering knuckle, and suspension linkage. The flexible lower wishbone consists of four interaction nodes, representing the connection to the frame, to the push-rod and to the steering knuckle. The flexible pushrod consists of two interaction nodes, representing the connection to the linkage and the lower wishbone. The modal basis of each flex body is derived by MSC. Nastran [19] and consists of 30 modes for the wishbone and 26 modes for the pushrod. Note that flexible MBS simulation in FreeDyn needs the so-called zero-inertia bodies. These bodies, with mass and inertia equal to zero, are needed to establish constraints between a flex bodies' interface node and any other body or ground. Hence, any boundary condition between a flex body and a rigid body brings up a set of seven coordinates and at least seven constraint equations, including the inner Euler-constraint [21]. Summing up, the model consists of four rigid, two flexible, and six zero-inertia bodies. It should be noted that flexible bodies, modeled in FreeDyn, consist not only of the flexible coordinates but also of seven physical coordinates which describe the rigid body motion. Hence, together with the 56 flexible coordinates, a total of 140 coordinates arise. As the vehicle's frame is not modeled, the suspension is fixed to the ground, bringing up a total of 84 constraint equations. The high number of constraint equations originates from the zero-inertia bodies, which are coupled to the flexible bodies using fix joints, with six constraint equations each. As the flexible lower wishbone has four coupling points, and the flexible pushrod two, 36 constraint equations arise between these points and the zero-inertia bodies. Further, due to the use of Euler parameters, each body brings up the need of one internal constraint equation, which causes another six constraint equations. The remaining 42 constraint equations result from the actual constraints between the single physical/zero-inertia bodies and the environment.
The suspension strut is replaced by its characteristic spring and damping force curves, and the frictionless model is actuated at the wheel hub by a frequency sweep from 1 to 15 Hz. Due to the chosen mass and damping parameters, the simulation passes through a resonant frequency of the model, see Fig. 3(b).
Investigating the POV plots in Figs. 6(a)-6(c) high valued drops occur after 22 translational and 44 rotational POVs. The POVs plot of the flexible coordinates shows several drops starting with six flexible coordinates, seven and eight flexible POVs. Again, the choice on the final set of reduced coordinates is not clear but requires trying out several combinations of reduced coordinates. In fact, the use of 22 translational, 44 rotational POVs, and six flexible POVs shows good convergence to the original solution. Again, no general statement can be given on how to identify the "right" set of reduced coordinates, as in this case the use of the flexible POMs associated with POVs 1 to 6 results in a reduced order model, which is as acceptable as if using flexible POMs 1 to 8. The generalized constraint reduction method identifies a total of 68 necessary and unique constraint equations, see Fig. 6(d). The reduced order model saves almost 70 % of the overall simulation time. Similar to the previous example, the POVs are further investigated to see whether any other POM constellation may lead to better or faster results. In Fig. 7, two possible reduced order models are compared to the obvious reduced dimension. Although the errors for characteristic model coordinates are still below 1 % limit for all reduced order models found, the overall simulation time does not improve.
This paper extends the model order reduction method, presented in [27], to models in arbitrary position in space. As the method originally proposed focuses mainly on rigid models, with constraints modeled in global axis directions, its application to models other than these is not as straightforward as expected. Especially when dealing with vibrations or resonance phenomena, the resulting tumbling behavior critically influences the constraint reduction potential. Therefore, we herein propose a generalization to the constraint reduction method in [27], which is based on the ideas of principal component analysis, showing that the original constraint reduction idea is a special case of the present one. The proposed constraint reduction method ensures that the resulting reduced order models are determined and wellconditioned. One no longer has to distinguish between internal and external constraint equations, as the constraint reduction method also allows the Euler-constraints to be reduced.
However, the reduced order models are no longer close to the minimal set of coordinates representation. Finally, it must be stated that, although the reduced order models show high consistency in the system states, the constraint forces of the reduced order model are no longer comparable to the original system. Hence, this method is not suitable to deliver the same joint reaction forces, which are then used for upcoming considerations like FEM simulations, as the original model.
Further, the accuracy of the reduced order models is closely related to the processed snapshots. In cases with critically changing excitation directions or model parameters, omitted constraints may become necessary. In such a case, it might be necessary to collect a new set of snapshots in order to renew the reduction subspace and the set of unique constraints.
The numerical examples, with which we demonstrated the method, are modeled arbitrarily in space. The proper orthogonal value plots suggest that for practical examples, the reduced order models are not necessarily indicated by a set of high-valued POVs or huge drops in magnitude, but could also be indicated by, e.g., flat spots. We point out that in such a case it may be necessary to identify the reduced order model by repetitive simulation runs, using a varying number of reduced coordinates. The numerical examples show that although vibrating systems and resonance phenomena complicate the identification of a reduced order model, the reduction method still performs well.