Abstract
The evaluation of the response function of the structural pounding problem is generally time-consuming if high-order systems are applied. The well-known modal truncation strategy is outstandingly efficient for a single linear ground-accelerated structure. However, for the analysis of the structural pounding problem, the classical modal truncation technique turns out to be ineffective as additional higher frequency motion due to possible contact impact occurs. This makes the determination of how many modes should be taken into account in order to obtain a required level of accuracy more difficult. Therefore, in this paper, a systematically controlled modal truncation strategy adapted to the seismic pounding formulation under consideration of high nonlinearity and nonsmoothness of contact problems is introduced. A comparative study of the classical and the controlled modal truncation technique is presented and a comparison with the commercial software package ABAQUS© is provided. It is shown that the computational accuracy is significantly improved when applying the new systematically controlled modal truncation strategy.
1. Introduction
It has been observed that the motion of buildings due to earthquake excitation may not cause as much damage to buildings as the result of possible pounding between them [1–4]. Pounding occurs if the separation distance between adjacent buildings is not sufficient in order to compensate for their relative movements. It is dynamically complex and even chaotic behavior that occurs if the eigenfrequencies of the single buildings differ considerably [5]. Therefore, numerical simulations of the dynamic pounding problem are necessary in order to understand the characteristics of structures based on which guidelines for structural design can be provided [6–8].
Due to the high number of degrees of freedom and the nonlinearity as well as nonsmoothness of the contact problem, the evaluation of the multiple structure pounding is time-consuming [9, 10]. As a result, the evaluation of the response to a whole frequency range can be computationally disproportionately expensive. Thus, low-order models are necessary that guarantee a certain level of accuracy.
Some studies have investigated model order reduction techniques based on modal truncation strategies and the proper orthogonal decomposition for structures with nonlinear materials. However, they applied those strategies only to dynamical single body systems without consideration of multibody pounding; see [11–17]. Only limited research studies have been conducted considering model reduction techniques for structural pounding problems. In the work of Zucca [18], another model order reduction technique based on dynamic substructuring technique is applied to the forced response of structures with contact interfaces. In the work of Géradin and Rixen [19], a nodeless dual-superelement is formulated to reduce the degrees of freedom of the unilateral contact problems. Most of the researches focus on the modeling and mitigation of the constitutive pounding law [1, 3, 8]; few of them attempt to apply model reduction techniques for the pounding problems.
The aim of this paper is to provide a low-order representation of structural pounding formulation by applying the classical and a systematically controlled modal truncation technique considering dynamic contact impact interaction. The theoretical principles and mathematical formulation as well as the solution methods in time domain are introduced in Section 2. A representative academic model, an ensemble of three adjacent frame structures, is demonstrated in Section 2.2. In Section 3, the numerical investigations of the dynamic model are presented. In order to verify the in-house software, a comparison with the commercial software ABAQUS© for the solution is given. Convergence studies and frequency response analyses as well as error estimations are carried out. The full and the reduced solutions are compared concerning efficiency and error of the displacements, acceleration, and contact forces. The effectiveness of the improved method is verified by performing a frequency response analysis applying different number of considered modes for the reduced computation. Additionally, the proposed controlled modal truncation strategy is examined on three real earthquake records. Conclusions and outlook of the demonstrated methods are presented in Section 4.
2. Mechanical Formulation
2.1. Initial Boundary Value Problem (IBVP) for Multiple Pounding
2.1.1. Strong Form
The problem statement associated with the kinematic linear frictionless contact problem for the three demonstrated deformable bodies is summarized in the following equations. The IBVP for a dynamic system with finite deformation must be satisfied for each body (cf. [20]):The boundary of each body is divided into three disjoint parts, shown in Figure 1:where the Dirichlet, Neumann, and contact boundary conditions are defined on , respectively. It should be pointed out that one body can have more contact surfaces; see Figure 1(b). Conventionally, one defines slave and master surfaces for the contact problem with superscripts and , respectively; see [20, 21]. Generally, the definition of master and slave is arbitrary. In the demonstrated formulation, Figure 1(b) plays a role as master with respect to both Figures 1(a) and 1(c). This convention will be kept to describe the academic test example and contact constraints in Sections 2.2 and 2.3.
Taking small deformation into consideration, the linearized Green-Lagrangian strain can be expressed asLinear-elastic material is introduced; therefore, the constitutive relation between strain and displacement under small deformation is given aswhere is the Cauchy stress tensor and the tensor is the forth-order elasticity tensor. Karush-Kuhn-Tucker (KKT) conditions are applied to describe the contact constraints in a strong sense, which are in case of normal frictionless contact [21]The scalar magnitude is the gap function between two participants, which is determined by using the closest point projection. The magnitude is the contact traction acting on the slave surfaces and the vector is the tangential traction describing the friction and other viscous tractions. In this paper, friction is neglected.
2.1.2. Weak Form
Based on the strong form of contact problems for multiple deformable bodies, the weak statement can be achieved in a variational sense using the principle of virtual displacement. By multiplying the virtual displacement as test functions with the local momentum balance in (1), one gets the weak form of the pounding system with frictionless contactwhere stands for the virtual work of the inertia column forces, is the virtual work of the internal stresses, is the virtual work of the external forces, and stands for the virtual contact work. They are given asIt should be pointed out that the weak form is also valid for each separate body . There are many possible ways to contribute to the contact virtual work such as Lagrange multiplier method, augmented Lagrange multiplier method, or mortar method [20, 21]. In this work, the contact constraints are regularized by the penalty method, which is introduced in Section 2.3.
2.2. An Academic Multiple Pounding Benchmark Example
A system of two three-story frames and one higher four-story adjacent frame is selected to represent the academic pounding benchmark example. The selected frames are modeled by column-beam elements. The frames are arranged with an initial gap of m. The widths and the heights of the frames are m, m, and m and m, m, and m. This is depicted in Figure 2. Different cross sections as well as different stiffness parameters are chosen for the columns and beams. The measurements and the material parameters of the cross sections are presented in Table 1.
The structures are discretized by two node structural beam/truss elements with three degrees of freedom per node, that is, axial and vertical local displacements and rotations (in total 978 degrees of freedom).
In order to keep the convention of indices for slave and master surfaces in Section 2.1, the contact domain between frame and frame will be marked as (slave side) and (master side) and, furthermore, (slave side) and (master side) for the contact domain between frame and frame , rather than , , and ; see Figure 2. An a priori master-slave separation is possible in this academic case.
2.3. Method of Contact Constraint Enforcement
Pounding between two adjacent objects is usually modeled by introducing the impact elements, which are widely used to describe the contact force if collision occurs (cf. [1, 3, 8, 22–24]). Mostly pounding elements with a linear or a nonlinear spring damper system, following the Hertz law, are applied (cf. [2, 25–29]). In this paper, a pounding element with a linear spring is chosen to describe the impact forces in both horizontal and vertical directions by interpolating the node tractions. This is equivalent to the penalty method, using the penalty parameter for the stiffness of the spring. This method enjoys the particular advantage of removing the constraints explicitly from the variational formulation; see [21]. The contact constraints, regularized by the penalty method, are briefly introduced.
As mentioned in Section 2.1, the slave side is labeled by index and the master side by index , respectively. Based on the node-to-segment formulation, the gap function mentioned in Section 2.1 is defined by solving the problem of the minimum distance between the slave point and the opposite point on the master surface of the discretized system; see Figure 3. The possible element of the master surface can be defined through the equationThe gap function is then defined aswhere is the parameter of the intersection point coordinate. The contact traction can be described as the product of the penalty parameter and the gap function . The nonconform meshes between two bodies are combined in this sense. Inserting the gap function into the weak form (7), the virtual contact work can be formulated aswhere the bracket is the Macaulay bracket, which meansEvaluation of the volume integration of the individual term in (7) delivers the explicit formulation of the set of equations of motion describing the discretized dynamic systemwhere , , and are the global mass and damping and stiffness matrices, respectively. Rayleigh-Damping is chosen for the damping matrix , which is proportional to the mass and stiffness matrix:where the second and the fifth modes are used to evaluate the coefficients and :with a damping ratio , rad/s, and rad/s. The contact force is then directly expressed as a function of the displacement field by the penalty method, which couples all degrees of freedom of all frames. The initial conditions are assumed to be zero for both the displacements and velocities.
2.4. Modal Truncation Strategies
The global set of equations of motion (12) is solved numerically. Dependent on the number of degrees of freedom and the performance of the contact search algorithms, this is generally time-consuming. Thus, the simulation of the system response for the whole range of frequencies is even much more challenging. Therefore, solution strategies applying the classical and a systematically controlled modal truncation technique are introduced.
2.4.1. Classical Modal Truncation
The undamped homogeneous set of equations of motion leads after inserting the solution ansatz function to the eigenvalue problem . The solution of the eigenvalue problem leads then to the spectral matrix , with the eigenfrequencies in ascending order and the corresponding eigenmatrix (cf. [12]). The transformation into the reduced space is physically motivated by the fact that, dependent on the frequency range of the excitation force, the whole deformation energy of the vibrating system is associated with a low number of vibration modes. In case of seismic ground excitation, it can a priori be assumed that the whole deformation energy is stored in the first couple of low frequency modes of vibration (cf. [13]). Based on modal analysis, the transformation matrix can be extracted from the sorted modal matrix according to correspondent eigenfrequencies of the full system with , which meansThe physical coordinate is transformed into the reduced space byAdditionally, left transformation by leads toConcerning the multiple structural pounding problem, the matrices and are block-wise diagonal and only the contact force is responsible for possible interaction. The solution of the eigenvalue problem including these matrices provides eigenvalues, which can be sorted in increasing order, but it is not possible to assign these eigenvalues to the frame system which they belong to without visible check. Therefore, concerning the multiple pounding problem, this is an uncontrolled modal truncation strategy. The motivation is to divide the whole mechanical problem into subsystem coupled by a contact term and investigate each subsystem separately.
2.4.2. Controlled Modal Truncation
The observation of the afore presented classical uncontrolled modal truncation method is that the system modes do not necessarily cover the essential modes of all bodies within the overall number of considered modes. Applying a systematic modal truncation strategy, the system modes are picked in a controlled way, which means that the main modes of each body and the unessential modes are truncated. In other words, the truncation will be carried out for each separate body. Therefore, the system is decoupled into subsystems and for each subsystem the following set of equations of motion must hold:For each body, the modal truncation technique is now separately executed with ; that is,The set of equations of motion of one body is rewritten by left multiplication of the reduced modal matrix asAssembling all of the controlled modal truncation matrices, one obtains one controlled modal truncation matrixand the controlled truncated set of equations of motion:Through this transformation, the global movement of the system is captured and, therefore, a more reliable low-order system is ensured.
2.5. Time Integration
The nonlinear dynamic system is solved using the central difference integration scheme (cf. [12]). Considering second-order Taylor forward and backwards expansion, the approximation of the velocity is obtained by the sum of the two expansions, and the approximation of the acceleration is obtained by the difference of the two expansions. Thus, the approximations of the velocity and acceleration by the central difference method are described byInserting the approximations in (12) leads to an explicit formulation for the displacement at the subsequent time step:As presented in (24), the displacement vector at time instant is evaluated based on the already known displacement vectors at time instants and . Therefore, an additional value for must be calculated before the start of the integration loop. This is realized by the second-order backwards Taylor extrapolation:The initial acceleration is obtained by inserting the initial conditions and into the global equation of motion (12). The coefficient matrix is positive definite and can be LU factorized a priori, which accelerates the solution procedure. The big advantage of explicit time integration schemes is that for each step no iteration is required. However, the disadvantage is that relatively small time steps must be chosen in order to ensure convergence. In other words, if the chosen time step exceeds a critical value, the solution function increases exponentially to infinity. The critical time step of the central difference integration scheme depends on the highest eigenfrequency of the dynamic system:where defines the highest eigenfrequency of the dynamic system. Applying the frequency based controlled modal truncation strategy, the highest eigenfrequency of the low-order system is significantly reduced. Thus, the critical time step is considerably enlarged. Applying this strategy, a much higher iteration time step can be chosen, which leads to a considerably smaller number of integration loops to be processed and a significant decrease of the total computational time.
The new strategy, presented in this paper, is a remarkably efficient explicit time integration scheme, as it is able to overcome the general downside of these types of methods, that is, computational instability. Implicit time integration schemes, as, for example, the Newmark method (cf. [12]), require, on the one hand, an iterative procedure within one time integration step, and, therefore, more computational effort within this integration step. On the other hand, in case of the right choice of certain parameters, they are unconditionally stable. Choosing an arbitrary large time step, generally, a converged solution function can be observed. However, applying very large time steps, a solution function is obtained, which converges to the static solution. This is caused by the numerical damping effect. Thus, in order to ensure response functions with sufficient accuracy, additional calculations can ensure a sufficiently accurate solution. In summary, a smaller number of integration steps must be applied, but the computational effort and storage per integration step are considerably higher compared to explicit time integration schemes.
Within one test calculation, the standard central difference scheme is compared with the new reduced integration strategy and, additionally, with the implicit Newmark method, using the commercial software package ABAQUS©. Through this approach, not only is the overall correctness of the implemented in-house code verified but also the computational efficiency is tested. The comparison of the three methods can be found in Section 3.5. A more detailed discussion on the advantages and disadvantages of the different time integration methods for the dynamic Signorini problem is given in the study of Doyen et al. [30].
3. Comparison of Methods
3.1. System Modal Analysis
The uncontrolled classical modes of vibration of the whole system are calculated. The first modes are depicted in Figure 4. It is observed that, due to different stiffness and characteristics of frames, the system modes alter from one frame to another frame in an irregular pattern. The disadvantage of the classical modal truncation technique is obvious: if one takes, for example, the first modes and cuts the rest modes, the distribution of modes will be concerning the substructures . A consistent distribution is therefore not achieved, which leads to an unsatisfactory low-order model.
The idea of the controlled modal truncation is motivated by the above introduced disorder of system modes. The aim of the systematic modal truncation is to decouple the whole nonlinear contact problem and to keep the important modes of each frame by executing the modal truncation for each frame subsystem separately. The choice of the modes can be empirically selected according to the number of stories of the frames or according to the frequency range of excitation. For instance, if the first ten modes are used and the rest of the modes are truncated, one can choose three modes for the first frame, four modes for the second frame, and three modes for the third frame. This is depicted in Figure 5. A consistent distribution of modes over the subsystem is achieved and a sufficient number of deformation patterns within the contact area is ensured in order to describe global contact behavior. These two methods are compared in detail applying harmonic and transient seismic excitations in the following sections.
3.2. Sinusoidal Seismic Excitation
The level of accuracy with respect to the number of considered modes for both the classical and the controlled modal truncation technique is studied under sinusoidal excitation. The results are focused on the displacement, the velocity, and the acceleration as well as the contact force of the horizontal degrees of freedom of the output node depicted in Figure 2. The ground acceleration is defined as with m/s2 and rad/s. In order to examine the accuracy of the strategies, six, ten, and fifteen modes of vibration are chosen and compared. We obtain unsatisfactory results for all values of interest if only the first six modes (with a mode distribution of concerning for the controlled modal truncation strategy) are chosen; see Figure 6. The deviation of the acceleration and the contact force is considerably large. Consideration of only six modes for both classical and systematic modal truncation is obviously insufficient for a reliable approximation of the response outputs.
(a) Displacement
(b) Velocity
(c) Acceleration
(d) Contact force
The results improve significantly if ten or more modes are taken into account. Regarding the classical modal truncation, if the first ten modes, illustrated in Figure 4, are chosen, the results improve considerably as shown in Figure 6 (red curves). However, noticeable errors for the acceleration and the contact forces are observed. By means of controlled modal truncation technique (with a mode distribution of concerning ), which is depicted in Figure 5, the displacements and velocities are almost the same compared to the full system. It can be concluded that the systematically controlled modal truncation technique delivers much better results than the uncontrolled modal truncation for every quantity of interest; see Figure 7.
(a) Displacement
(b) Velocity
(c) Acceleration
(d) Contact force
If 15 system modes (with a mode distribution of concerning for the controlled modal truncation strategy) are taken into consideration, the results between the full and the reduced systems agree almost; see Figure 8. The critical time step is smaller than the truncated system with ten modes. The convergence study with respect to the number of modes is carried out in the next section. Error analysis of the two techniques is given in Section 3.4.
(a) Displacement
(b) Velocity
(c) Acceleration
(d) Contact force
3.3. Frequency Response Analysis
The frequency response analysis is essential for a better understanding of system characteristics in earthquake engineering. It is also the most commonly used method regarding the design of structures subject to earthquakes [31]. The frequency response analysis investigates how the maximum response of the displacement, the velocity, and the acceleration as well as the contact force of the structure behaviors subjected to harmonic excitations with different frequencies. For the demonstrated model, the driving angular frequencies vary from rad/s to rad/s. The results regarding the maximum displacements, the velocities, and the acceleration as well as the impact forces versus driving angular frequencies applying the first ten modes are shown in Figure 9. It is observed that the extremal values are reached near rad/s. The numerical calculations of the first five eigenfrequencies are , , , , and rad/s. The extremal values are not located near these frequencies as complex pounding interaction influences the overall system behavior.
(a) Displacement
(b) Velocity
(c) Acceleration
(d) Contact force
It is observed that the controlled modal truncation technique approximates the results significantly better than the standard modal truncation technique. The maximum displacements, velocities, and contact forces for the whole range of interest of the excitation frequencies are well approximated by the controlled modal truncation method. The computational effort for each test frequency of the full system is about 315 seconds. However, the truncated system needs only 15 seconds. As the whole evaluation is based on the representative node, it is necessary to carry out a global error estimation, considering all degrees of freedom.
3.4. Error Estimation
As shown in Section 3.1, the results are mainly dependent on how many modes are chosen. The root-mean-square deviation (RMSD) of the node output histories is evaluated and plotted in Figure 10. The RMSD is defined aswhere is the reference solution of the full system and is a reduced solution at the time step . The RMSD errors of the standard method are larger than the controlled modal truncation if the applied number of modes is smaller than 20. If more than 20 modes are taken into consideration, the RMSD of both methods are almost the same.
(a) Displacement
(b) Velocity
(c) Acceleration
(d) Contact force
3.5. Comparison of the New Strategy with the Commercial Software Package ABAQUS©
In order to verify the correctness and efficiency of the implemented in-house software and the new strategy presented this paper, the full system response is compared with the commercial software package ABAQUS© and the reduced system using the new reduction strategy. Applying ABAQUS©, the same discretization compared to the in-house simulation is applied and the same penalty parameter, that is, stiffness parameter of the linear pounding element, is chosen within the contact module. The test is carried out in two steps.
In the first step, a statical analysis is performed. A horizontal force is applied at the top left node of the first frame system, where the magnitude of the force is chosen to be large enough in order to ensure contact between frames one and two and contact between frames two and three. The displacements of numerous representative nodes of the first, the second, and the last frame are evaluated by using ABAQUS©. The results of the commercial software package agree perfectly with the in-house software tool. This verifies the correctness of the stiffness matrix and the contact formulation of the in-house software.
In the second step, the dynamic response of the structural system is evaluated, using the in-house implementation and the commercial software package ABAQUS©. The comparison of the displacement of the representative node, using ABAQUS© and the in-house software, is presented in Figure 11 by the blue dashed and black lines, respectively. Additionally, the solution of the new proposed controlled modal truncation strategy is depicted by the red dotted line in this figure.
In Figure 11, it can be seen that the in-house software agrees well with the result of the commercial software tool. Thus, the correctness of the calculation is verified. However, small deviations are observed, which are caused by the fact that, for the ABAQUS© simulation, the implicit Newmark integration is chosen, applying a time step of 0.001 seconds, and, for the in-house implementation, the explicit central difference method is used. It is also observed that, for the implicit integration scheme of the commercial software package, a maximum time step of 0.001 seconds must be chosen; otherwise, numerical damping is considerably high, and an inaccurate solution is observed. Detailed comments about the advantages and disadvantages of the two methods are discussed in Section 2.5. The computational cost of the in-house software, using an integration time step of 0.0001 seconds for the explicit integration scheme, is comparable with the one of commercial software ABAQUS©. On the same machine, the calculation time using the in-house tool is about 315 seconds and the calculation time using the commercial software is about 297 seconds. However, the calculation time, applying the novel strategy, is only 15 seconds. It is observed that the new method shows an outstanding efficiency by preserving a required level of accuracy. Thus, the new method provides the possibility of evaluating high dimensional dynamic multiple pounding systems, which demand generally very high computational effort, using the standard integration scheme.
3.6. Excitation under Different Earthquake Measurements
Regarding general problems in earthquake engineering, ground excitations are not harmonic. They are transient. Therefore, the presented methods are verified applying different earthquake excitations. Two earthquake test samples depicted in Figure 12 are chosen. The calculation is performed taking into account ten modes (with a mode distribution of concerning for the controlled modal truncation strategy). Data of the earthquake measurements are presented in Table 2.
(a) Kobe earthquake
(b) Bam earthquake
The numerical results of the representative node, presented in Figures 13 and 14, reveal that the systematically controlled modal truncation technique provides reliable approximations of the system response. Taking into account the displacements and velocities compared to the full benchmark response, a higher accuracy is achieved by the controlled modal truncation technique. Concerning the acceleration and contact forces, especially at the contact time instants, the low-order model must be improved by taking more modes for both the classical and the controlled modal truncation methods.
(a) Displacement
(b) Velocity
(c) Acceleration
(d) Contact force
(a) Displacement
(b) Velocity
(c) Acceleration
(d) Contact force
4. Conclusion and Outline
In this paper, the classical and a systematic modal truncation technique are studied, taking into account computational contact impact problem. Pounding is realized applying linear pounding elements within a node-to-surface contact formulation. Due to the pounding impact the classical modal truncation technique is not sufficient to approximate the system response since the unconsidered higher frequency vibration modes play also an important role for the system response. The computational results show that the systematically controlled modal truncation is able to deliver a more reliable approximation of the full systems, especially for the displacements and velocities. This is also verified applying real earthquake excitations. The systematically controlled modal truncation technique makes the frequency response analysis for complex structural systems with a high number of degrees of freedom more feasible with consideration of the computational effort. It should be pointed out that the stated method is of course not limited to the seismic pounding problem only, but to a much broader field of dynamical contact and impact problems.
It must be added that the systematically controlled modal reduction technique is effective for materially linear systems. In the present paper, it is shown that it can also be efficiently applied for dynamic pounding problems. However, it is still necessary to increase the accuracy of the presented method, especially, with respect to the acceleration and the contact force functions. The next step concerning future research should also concentrate more on local vibration near the contact area. Here, it should be more focused on substructuring techniques, performed by, for example, the Craig-Bampton method, where in this paper defined subsystems are again subdivided into substructures.
In case of material nonlinearity, further extensions of the above-mentioned strategies are necessary. Generally, the method of proper orthogonal decomposition was found to be useful for creating model order reduced descriptions, which are able to provide accurate approximations of nonlinear dynamic systems. The objective for upcoming studies, regarding the accurate low-order description of pounding problems, is to combine the proposed strategy with the Craig-Bampton method (cf. [18]) and with the proper orthogonal decomposition (POD), which has been presented in recent studies [11, 12]. In particular, the frame subsystems, introduced in this paper, can be described by constraint modes, concerning the global motion of the structure, and internal modes, describing the local behavior of the structure within the contact area. Then, the nonlinear system behavior can be approximated using proper orthogonal modes, which are evaluated from an information matrix (snapshot matrix [13]), which contains already existing solution vectors of the materially nonlinear system and the constraint as well as the internal modes of the frame subsystems.
Conflicts of Interest
The authors declare that they have no conflicts of interest.