Dynamic substructuring methods encompass a range of techniques, which allow for the decomposing of large structural systems into multiple coupled subsystems. This decomposition of structures into smaller domains has the principle benefit of reducing computational time for dynamic simulation of the system by considering multiple smaller problems rather than a single global one. In this context, dynamic substructuring methods may form an essential component of hybrid simulation, wherein they can be used to couple physical and numerical substructures at reduced computational cost. Since most engineered systems are inherently nonlinear in nature, particular potential lies in incorporating nonlinear treatment in existing sub-structring schemes, which are mostly developed within a linear problem scope.
The most widely used and studied dynamic substructuring methods are classical linear techniques such as the Craig-Bampton and MacNeal-Rubin methods. These methods are widely studied and implemented in many commercial FE packages. However, as linear methods they naturally break down in the presence of nonlinearities. Recent advancements in substructuring have involved the development of enrichments to the linear substructuring methods, which allow for some degree of nonlinearity to be captured. The use of mode shape derivatives has been shown to be able to capture geometrically non-linear effects as an extension to the Craig-Bampton method. Other candidates include the method of Finite Element Tearing and Interconnecting. Linear substructuring methods have been used in several cases for hybrid simulation, rendering additional benefits in removing non-physical high frequency modes which may be excited due to controller tracking errors in hybrid simulations.
In this work, a virtual hybrid simulation is presented in which a linear elastic vehicle frame supported on four nonlinear spring damper isolators is decomposed into separate domains. One domain consisting of the finite element model of the vehicle frame, which is reduced using the Craig-Bampton method to only include modes relevant to the structural response. The second domain consists of the nonlinear isolators whose restoring forces are characterised by nonlinear spring and damper forces. Coupling between the models is carried out using a Lagrange multiplier method and time series simulations of the system are conducted and compared to the full global system with regards to simulation time and accuracy.
Keywords: Dynamic Substructuring, Hybrid Simulation, Nonlinearity
In this study we investigate implementation of dynamic substructuring schemes with hybrid simulation. Dynamic substructuring (DS) methods are used to decompose dynamic models into several separate regimes or substructures. The global model of the system is always recoverable via use of coupling constraints enforced between substructures. The key benefit of DS methods lies in the capability of analyzing each substructure separately and possibly in parallel, eventually assembling the global solution; this can greatly reduce the computational burden of a given analysis. However, further benefits are generated from the implementation of different analysis schemes on different regions of a structure. A key motivation of the above being the ability to individually treat nonlinear regions, without the requirement for adoption of a global nonlinear solution method [Voormeeren2011].
Dynamic substructuring methods are typically combined with reduction methods, such as Component Mode Synthesis (CMS), whereby reduced order models of the individual substructures are determined, which further decreases the computational load of the solution process. Most traditional techniques such as the Craig-Bampton and Rubin-Macneal methods [VanderValk2010], rely on a modal decomposition of the substructures whereby only a subset of mode shapes is retained -those relevant to the structural response- whilst the remaining modes are ignored. The commonly adopted linear schemes are naturally not suitable for the treatment of nonlinear systems. Recent advancements in substructuring have involved the development of enrichments to the linear substructuring methods, [Wu2018], [Kim2015], [Kim2017], which allow for some degree of nonlinearity to be captured. The use of mode shape derivatives has been shown to capture geometrically non-linear effects as an extension to the Craig-Bampton method [Wu2018].
An alternative to these methods is the so-called finite element tearing and interconnecting (FETI) technique [Farhat1991]
in which a coupling method making use of Lagrange multipliers is used, whilst the reduction of the substructures can be carried out separately. This allows for the use of any generic meta-modelling technique for the reduction of the non-linear substructures, such as a non-linear autoregressive models or a Gaussian process regression. These techniques can be shown to provide accurate representation of a complex model over a given input range of interest[ChatziCostas], [Mai2016], [inproceedings]. The FETI method can also allow for the coupling of nonlinear substructures or linear to nonlinear substructures. This has previously been used in hybrid simulations [Abbiati2019] and was shown to perform well.
Hybrid testing, which forms a focus of the European ITN Project DyVirt, on Dynamic Virtualization, involves the setup of coupled dynamic tests in which a numerical system in a computer is coupled to a physical component in a dynamic testing rig. The system is then solved as a whole by integrating the numerical system, calculating predicted forces at the interface between the structures and applying these to the physical system. The response of the physical system is then measured and fed back into the equations of the numerical system [BensonShing2008]. The benefit of hybrid testing is that simple or well characterised components of a system can be represented numerically, whilst more complex or critical components can be physically tested. This allows for a reduction in cost and complexity of experimentation against full scale testing with an increase in fidelity and reduction in uncertainty with respect to purely numerical modelling. The coupling between the physical and numerical components within hybrid testing is a parallel of the DS problem, which attempts to couple numerical systems; as such, similar methods to those in DS can be used [Wagg].
Within the domain of hybrid testing the“gold standard” is considered to be real time hybrid testing . Real time testing implies that the time scales of the numerical and physical system are the same; and that this time scale matches that of “wall time”. The benefit of real time testing is that it allows for rate dependent behaviour and nonlinearities of the physical system to be accounted for [BensonShing2008]. However as we attempt to approach real time testing the problem of the time required to integrate the numerical substructure becomes an issue. In order to carry out a real time hybrid simulation, each step of the numerical integration must be able to be carried out within a smaller time increment than the integration time step. This becomes a problem when numerical substructures of increasing complexity are considered and leaves the options of either increasing the time step of the integration which can have negative effects on accuracy and stability [Pegon2008], or finding a reduced order model of the numerical system which can be solved more quickly. However, there remains considerable potential for the use of model order reduction techniques in combination with hybrid simulation.
Application of CMS methods in hybrid testing have been made as in [Abbiati2019], although in this case the key purpose of the dynamic substructuring was to remove high frequency modes from the numerical system to yield improved stability characteristics, rather than to improve the speed of calculation.
Another issue which arises in hybrid simulation is that as the integration time step is increased to allow for real time simulation, the time step for control of the physical system is also increased to match it in order to couple the systems. This increased time step of the physical system can result in a non-smooth operation of the actuator due to control commands being sent at such a large interval. This can cause issues with regards to the fidelity of the applied loads and issues of stability in the solution [doi:10.1002/eqe.743]. As such, when considering the solution time step not only must the numerical integration be considered but rather also the physical control system. This problem has inspired algorithms which allow for adoption of individual time steps in the physical and numerical subsystems [Abbiati2018], [Abbiati2019a].
The key elements used in the virtual hybrid simulation we present herein are the model order reduction carried out on the linear ”numerical” system, and the integration algorithm which solves the systems and imposes coupling between them.
2.1 Craig-Bampton Reduction
The Craig-Bampton is one of the most prominent and long standing techniques in CMS and is already implemented in numerous commercial finite element softwares [Craig1968].
In the Craig-Bampton method each sub-structure can be reduced individually and then usually the reduced sub-structures are coupled together using a primal assembly. However, in this case only the linear system will be reduced whilst the coupling procedure will be carried out within the integration algorithm.
The key aspects of the reduction are that the degrees of freedom (DOFs) of each substructure are partitioned into internal DOFs,, and external DOFs, , i.e., those which lie on the interface between substructures.
The internal DOFs, which in most cases are the majority, can be significantly reduced using a modal decomposition from the eigenvalue problem in Equation2. Given that for structural dynamics the lowest eigenmodes are dominant, the majority of these modes can be discarded, with only the lower modes, retained. These are known as the fixed interface normal modes.
In addition to the fixed interface normal modes, constraint modes, , are also used to characterise the boundary DOFs and the static response. Constraint modes correspond to the static deformation shape due to a unit displacement applied at a boundary DOF. As such there exist as many, as there exist boundary DOFs. These are calculated as follows.
Using both the fixed interface and constraint modes, the original high dimensional coordinate set can be approximated on a significantly lower coordinate set, with the reduction of DOFs dependant upon the number of fixed interface modes that are discarded. The reduced coordinate set consists of the generalised internal coordinates and the retained boundary coordinates .
The mass and stiffness matrices of the sub-structures can then be reduced by projection onto these reduced coordinate sets.
In the Craig-Bampton method the reduction basis consists of the generalised internal coordinates, whilst the boundary coordinates are retained as physical; this simplifies the coupling procedure.
2.2 Integration and Coupling
The coupled integration scheme used herein is inspired by that presented in [Abbiati2019a], this comprises a technique designed for hybrid simulation in which each of the substructures is solved as a free solution over the time interval with coupling enforced at the end of the time interval using a method similar to the FETI method [Farhat1991]. The partitioned design of this integration method allows it to be applied in a hybrid test in that, the 2 systems are solved in parallel with coupling only enforced at the end of each time step. This means that in the context of a hybrid test, displacements can be enforced on the physical system, with the resulting restoring force being measured and being used to solve the free problem in the physical system and finally coupling being enforced at the end of the time step.
The coupled equations of motion for the numerical and physical systems to be integrated are shown in equations 6 and 7. These present the typical equations of motion of a nonlinear mechanical system in state space form. However, these equations are augmented with the term , this represents the coupling forces between the two systems with being the Lagrange multipliers giving the interface force intensities and a boolean matrix which locates the interface forces on the overall system.
Each of the physical and numerical subsystems is solved over a given time interval. This is initially carried out as a free problem, i.e., ignoring the coupling of the other substructures within this interval. For the free solution of the substructures a trapezoidal integration rule was used which corresponds to a parameter =0.5. The trapezium rule integration uses a prediction and a correction step in order to step forward in the integration. The prediction step as indicated in equation 8 gives the predicted state at the next integration step based on the previous state and state rate values, the integration parameter and the time step .
Within equation 8, indicates the predicted state of the free solution at the next time step which will be corrected to find the final free solution . The Predicted state is made based on the solutions from the previous coupled time step: and .
In the above equations
is the external force vector,evaluates the potentially nonlinear restoring force vector, consisting of the elastic and damping forces, evaluated at state and D is a matrix as formed in equation 10 wherein is the mass matrix, and the tangential restoring force vector linearised at zero state. This tangential restoring force vector was found using the ADIGATOR toolbox for automatic differentiation [Patterson2013Aeo]. This toolbox allows the tangential stiffness matrix to be found numerically inclusive of systems with nonlinear restoring forces.
The correction step is then made as in equation 11 making use of the predicted state and state rates found previously.
Having found the free solution for each of the substructures individually and in parallel, the coupling can then be enforced between them. The coupling is enforced using a dual assembly procedure, by enforcing equal and opposite forcing at the boundary DOFs using Lagrange multipliers. Firstly, a matrix is assembled, as in equation 12, which forms the so called Steklov-Poincaré operator [Quarteroni]. This relates the discrepancy in state at the boundary nodes to the Lagrange multiplier intensities.
In equation 12 the matrices and are boolean matrices, which enforce compatibility of the substructures [VanderValk2010] for the numerical and physical systems respectively, whilst and are boolean matrices, which locate the interface forces for the numerical and physical structures [VanderValk2010].
Having derived the Steklov-Poincaré operator, which is assembled only once at the beginning of the integration procedure, the Lagrange multipliers can be identified based on the free solutions of each of the substructures as in equation 13.
The link solutions for each of the substructures are then found, these represent the effect on the state of the substructure of the interface forces due to coupling. These link solutions are found by first finding the state rate due to the interface forces as in equation 14 and then the link state using this state rate as in equation 15.
Finally the global solution for each of the substructures taking into account the coupling between them is found as in equation 16 by summing the free and link solutions for each of the substructures.
2.3 With Sub-cycling
In order to ameliorate the problem of actuator smoothness in hybrid simulations, wherein the numerical integration time step is necessarily larger than the control frequency of the actuator, a set of methods exist in which the physical system is controlled at a smaller time step than the numerical system, with coupling enforced at the larger numerical time step [Abbiati2019a]. In connection with the above described algorithm, the sub-cycling method can be implemented in the calculation of the free solution of the physical substructure.
Firstly, a number of sub-cycles is decided upon, defined as the ratio of the numerical time step to the physical time step. The integration scheme with sub-cycling operates similarly as described previously except that the free solution for the physical substructure, as described by equations 8-11 above, is now replaced by the following procedure.
The sub-cycling loop runs over the number of sub-cycles . It is initiated by the coupled solution at the previous numerical system time step. The loop ends by finding the free solution at the next numerical time step. Equation 17 shows the prediction step of the trapezium rule in which is the time step of the physical system.
Equation 18 calculates the state rate prediction in which it is worth noting that the Lagrange multipliers from the last coupling time step are now included in the solution and that the D matrix is formed as in equation 10 albeit with the finer time step of the physical system used in the assembly.
Finally, the corrected free solution is calculated as in equation 19.
The sub-cycling loop then ends after sub-cycles whereby the prediction of the free state of the physical system is found and the coupling procedure can then be followed as described in the previous section.
3 Case Study
The case study considered herein is a steel vehicle-like frame considered to be linear elastic, which is mounted on four nonlinear suspension structures. Computer models of the frame and suspension structures were created based on physical specimens. The frame structure model, along with the mounting points are presented in figure 1.
The frame structure is considered to be linear elastic and is meshed using a combination of quadrilateral shell elements and hexahedral solid elements with a total of 45,564 DOFs. The material properties of the frame are a Youngs modulus of , a Poissons’ ratio of 0.3 and a density of . This model has previously been calibrated to a physical structure, details of which along with dimensions of the frame can be seen in [Giagopulos2006]. The nonlinear suspension systems consist of a base excited mass representing the wheel and its interaction with the road surface, alongside a linear spring and nonlinear damping element which connect the wheel mass with the frame structure. These nonlinear suspension elements have also been identified from a physical component connected to the physical frame in a previous study [Tsotalou2017]. The best model found to approximate the non linearity in the damper was found to be a dry friction type damping force [Giagopulos2006]. The restoring force equations of the nonlinear suspension system are shown in equation 20, in which represents the linear stiffness coefficient, the linear damping coefficient and, and the nonlinear damping coefficients. The coefficients are summarised below.
The form of the nonlinearity present in the damper is illustrated in figure 2 where it is seen that the frictional nonlinearity delivers a more pronounced effect around the zero point, whilst further from zero more linear behaviour is observed.
Figure 3 presents the linear FE mesh of the frame structure and a diagrammatic representation of the nonlinear suspension systems. In the case of this virtual hybrid simulation, the linear frame is considered to be the numerical system which is to be reduced in order, whilst the suspension systems are considered to be the physical systems to be solved in full order.
4 Virtual Hybrid Simulation
A virtual hybrid simulation was carried out in which the case study described above was considered. The linear frame structure was considered to be the numerical system with the four suspension substructures considered as the physical systems. The linear system was first reduced using a Craig-Bampton reduction, and then the partitioned integration procedure was implemented on the substructures.
4.1 Reduction of the Linear Frame
The linear substructure was reduced according to the methodology described in section 2.1. As described earlier, a subset of fixed interface mode shapes were retained, which represent the dynamic behaviour of the structure within the frequency region of interest. The reduction basis was formed on the 30 lowest fixed interface normal mode shapes of the frame, along with 4 constraint modes each corresponding to the vertical DOFs at the 4 mounting points of the suspension systems.
Figures 4 and 5 demonstrate the fidelity of the reduction basis relative to the full model of the frame structure. In figure 4 a comparison is made between the 20 lowest natural frequencies and modal shapes between the full model and the CB reduction. Figure 3(a) shows the cross modal assurance criterion [Pastor2012] between the mass orthogonalised mode shapes of the full system and the CB reduced system. This gives a metric for the orthogonality of the mode shapes to one another and is calculated as in equation 22.
For identical mode shapes the diagonal of the plot should show one whilst all off diagonal elements should be zero. The plot demonstrates the fidelity of the mode shapes, as all diagonal elements are very close to one whilst all off diagonals are close to zero. Figure 3(b) compares the first 20 natural frequencies of the reduced and full systems, along with the normalised mean square error (NMSE) between them. The natural frequencies of the systems are very similar with none of the first 20 natural frequencies exhibiting an error larger than 0.1%
compares the time history response of both the reduced and full order systems to the same input. The input used in this case was a band limited white noise signal between 0 and 200 Hz with a variance of 1 which was applied in the vertical direction at the four nodes to which the suspension systems connect. This frequency region was chosen as the region of interest as ISO 8608, detailing the classification of road surface profiles, considers the region of interest to be up to a spatial frequency of, at an assumed vehicle speed of . This corresponds to a maximum frequency of interest of [Mucka2017]. The frequency range of interest was then selected to be comfortably above this value so as to ensure fidelity within the region.
The reconstruction of the time series response is near perfect for this input with a mean square error value of . According to both the eigenvalue and time series analyses the Craig-Bampton reduced system shows excellent fidelity to the full order model when excited. In both cases a Newmark integration scheme was used in the time series simulation.
4.2 Comparison of Monolithic and Partitioned Solutions
To assess the performance of the partitioned solution of the reduced order system, it was compared to a monolithic solution of the full system with regards to both fidelity and simulation time. The inputs chosen for the simulation were all multi sine excitations corrupted by white noise. The inputs were applied as base motion to the masses in the suspension substructure. Whilst each of the suspension substructures were excited with the same multi-sine content, they differed in the white noise corruption. Figure 6 shows an exemplar signal applied to one of the four suspension systems.
Figure 7 compares the response of the DS reduced system to the monolithic solution produced in Abaqus. The solutions show good agreement although some discrepancy can be seen particularly around the peaks. The mean squared error between the signals was showing the good fidelity. Some discrepancy is expected between the responses firstly due to the reduced linear system not performing identically to the full linear system. Additionally, the tangential stiffness matrix which is used for the coupling procedure introduces errors as it is effectively a first order linearisation of the coupling process.
In both the monolithic and partitioned solutions displayed, a time step of was used and 1 second of time history was simulated. This simulation from the monolithic solver required 1557 s compared to 0.15 s for the DS reduced system. This demonstrates the massive time savings possible using the reduction scheme. Additionally, it is noteworthy that the DS solution was found in faster than real time hence demonstrating its potential for real time hybrid simulation. It is worth noting that there are additional offline calculations required for the DS solution such as, the CB reduction and the forming of tangential stiffness matrices, the total offline time for these operations was 109 seconds. As such the overall time was still considerably less for the DS reduced system and furthermore the offline times are not of great concern with regards to hybrid simulation.
Figure 8 shows the solution of one of the boundary DOF of the numerical system and the boundary DOF of the physical system to which it is coupled. There is very close agreement between the two solutions showing that the coupling procedure was successful. However, figure 7(b) shows that there is a slight discrepancy indicating that the coupling procedure does not enforce hard compatibility between the two substructures but rather that a slight discrepancy is allowed as a result of the dual assembly [Rixen2004].
The effect of the addition of the sub-cycling of the physical system in the algorithm was also investigated. Firstly a comparison was made between the non sub-cycling algorithm with a time step of and the sub-cycling algorithm with 10 sub-cycles and a time step of this comparison is presented in figure 9. The overall response of the two algorithms are similar however, it is clear that there is some loss in fidelity. Overall, the sub-cycling algorithm exhibits stability issues. Previous implementations of the method in [Abbiati2019a] and [Abbiati2018] adopted an algorithm with numerical damping, as opposed to the undamped trapezium rule method used here. It is thought that the use of a damped algorithm could ameliorate these issues.
The algorithm shows potential however given that the sub-cycling algorithm can allow the time step of the physical system to be decreased by order 10 with an increase of computational time of approximately 2 times. It is also evident from figure 8(b) that the sub-cycled system did exhibit a smoother response in the physical system, which should improve the control of the actuator.
Figure 10 presents the response of a coupled boundary DOF on the numerical and physical substructures. Compared to the case without sub-cycling there is a clearly discrepancy, this discrepancy again highlights the capacity of the dual coupling procedure to allow for a soft coupling procedure. As seen in figure 9(b) the response of the physical system is initially very uneven which seems to cause the discrepancy which is maintained throughout. This unevenness is once again thought to be a result of the instability in the sub-cycling procedure, as highlighted previously.
A dynamic sub-structuring process has been presented, which can be used to couple linear and nonlinear systems in a manner suitable to hybrid simulation. It was demonstrated that this method can considerably decrease the required computational time and was also shown to reduce the integration time sufficiently to implement a real time hybrid simulation. The reduced system was also shown to deliver good fidelity when compared to the full monolithic system. The dual assembly procedure carried out was shown to enforce coupling between the substructures without enforcing hard compatibility and allowed an amount of discrepancy.
The DS integration algorithm was also implemented using a sub-cycling procedure whereby the physical system was integrated at a time step smaller than that of the numerical system. This method had considerable problems with instability of the physical system but it is thought that with the use of an algorithm with numerical damping these problems could be improved. It was shown that the use of the sub-cycling method resulted in a smoother response in the physical system which has advantages for the control of the actuator related to the physical system.
*2*B This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 764547