A new model order reduction strategy adapted to nonlinear problems in earthquake engineering

Summary Earthquake dynamic response analysis of large complex structures, especially in the presence of nonlinearities, usually turns out to be computationally expensive. In this paper, the methodical developments of a new model order reduction strategy (MOR) based on the proper orthogonal decomposition (POD) method as well as its practical applicability to a realistic building structure are presented. The seismic performance of the building structure, a medical complex, is to be improved by means of base isolation realized by frictional pendulum bearings. According to the new introduced MOR strategy, a set of deterministic POD modes (transformation matrix) is assembled, which is derived based on the information of parts of the response history, so‐called snapshots, of the structure under a representative earthquake excitation. Subsequently, this transformation matrix is utilized to create reduced‐order models of the structure subjected to different earthquake excitations. These sets of nonlinear low‐order representations are now solved in a fractional amount of time in comparison with the computations of the full (non‐reduced) systems. The results demonstrate accurate approximations of the physical (full) responses by means of this new MOR strategy if the probable behavior of the structure has already been captured in the POD snapshots. Copyright © 2016 The Authors. Earthquake Engineering & Structural Dynamics Published by John Wiley & Sons Ltd.


INTRODUCTION
The evaluation of the response history of a structure in the time domain is one of the main topics in earthquake engineering and structural dynamics. It is a common practice to create simple structural models, for example, multistory shear frames, which should be able to describe the structural behavior and peculiarities of the real structure. This approach leads to useful results for the investigation of rather simple and uniform structures in order to come to meaningful engineering decisions regarding structural resistance. On the contrary, the analysis of complicated systems can require the application of nonlinear high-order systems, as a characterization by a low dimensional structural model could lead to an oversimplification, that is, important motion patterns could be ignored. Therefore, an effective strategy is to obtain a set of a low number of 'important' equations of motion that approximates the high-dimensional nonlinear dynamical system as accurately as possible, that is, model order reduction (MOR). 539 of nonlinearities. The method is applied to the dynamic model of a realistic building structure. The building is erected on friction pendulum bearings for the sake of seismic isolation to minimize the transferred acceleration to the building during an earthquake. A three-dimensional dynamic model of the base-isolated structure is derived by implementing the FE model of the structure and the bidirectional friction pendulum systems. The paper then deals with the nonlinear dynamic model of the base-isolated structure. Thereafter, the new POD-based MOR strategy and the example of its practical application is provided. Finally, the results and conclusions are outlined.

EARTHQUAKE EXCITATIONS
Within the presentation of the new strategy and the application to a realistic building structure, a set of six earthquake excitation records is chosen for the numerical demonstrations. The earthquake records are applied in fault-parallel and fault-normal directions. The excitation set includes the Bam earthquake (2003) in Iran and the following five representative events in California, USA: Northridge Rinaldi (1994), Imperial Valley (1979), Landers (1992), Loma Prieta (1989), and North Palm Springs (1986). Table I presents a list of the events taken from the Pacific Earthquake Engineering Research Center (PEER) [22]. Fault-parallel is defined in x-direction and fault-normal in y-direction. Concerning the Bam event, only a one-dimensional record is available; therefore, an excitation attack angle of 30 ı with respect to the x-axis is chosen.

NONLINEAR MODEL ORDER REDUCTION AND THE POD-MATHEMATICAL FORMULATION
The n-dimensional set of equations of motion of a structure with nonlinear material behavior excited by horizontal components of ground acceleration is expressed as (cf. Chopra [23]) where M and C are mass and damping square matrices of order n and R.x/ is the nonlinear internal restoring force vector dependent on the displacement x with the dimension n 1. The right hand side of the set of equations of motion describes the earthquake excitation term, while R x g and R y g denote the ground acceleration in x and y direction, and f j ; .j D x; y) are the influence vectors in the corresponding direction, that is, f x .x i / D 1 ; f y .y i / D 1 ; i D 1 : : : n; (2) at the global x and y DOF of all nodes, whereas the other components of f x and f y are zero. Thus, i describes the number of nodes of the FE discretized structure. This general approach indicates that in this paper, the ground acceleration in the corresponding direction, that is, x and y direction, is equal in all structural support points. In the following equations, the term on the right hand side of the set of equations of motion (1) is denominated by F.t/, which has the unit of a force. Nonlinear systems, as they are depicted in Equation (1), have generally to be solved by the application of a numerical algorithm, that is, a step by step procedure in the time domain in order to obtain the response of the structure as a function of time. The application of a numerical method inevitably leads to the existence of computational effort if n is a large number. Therefore, the approximation by a low-dimensional description of the system seems to be useful, namely the application of MOR. The main goal of MOR techniques is primarily to define a transformation matrix T 2 R n m ; m n to approximate the displacement vector x 2 R n through a reduced coordinate vector q r 2 R m by the relation (cf. Koutsovasilis and Beitelschmidt [24]) x D Tq r ; ( such that the dynamic properties of the system are preserved, and the error is small. The notation of the variable m 2 N is the dimension of the reduced subspace. The projection of the nonlinear system defined by Equation (1) onto that subspace leads to another second-order ordinary differential equation (cf. Koutsovasilis and Beitelschmidt [24]) where m r D T T MT and c r D T T CT 2 R m m are mass and damping matrices, and f r D T T F.t/ 2 R m 1 is the force vector in the reduced subspace. It should be noted that the reduced system matrices m r and c r are generally not diagonal. The vector of the restoring forces in the reduced subspace is r D T T R.x/ D T T R.Tq r /: Consequently, one necessity of nonlinear MOR is the evaluation of the vector of the restoring forces in the physical (full) coordinate at every time step. Modal truncation is a widely-used tool and an effective method for order reduction of linear systems in the field of earthquake engineering. An accurate approximation of the response history is achieved by applying a small number of lower structural modes proportional to the number of DOF. An early and successful application of modal truncation to problems involving small local nonlinearities is presented in [25] and [26]. In this work, the objective is to find a new strategy that is applicable to nonlinear systems in a similar manner to modal truncation. The approach is to define a set of deterministic modes that can be evaluated from the information of an existing response history of the structure. Consequently, this set of modes contains nonlinear motion patterns if the structure shows nonlinear response behavior to the excitation.
The proposed strategy is established based on the theory of the POD method. Generally, the POD (cf. [11], [27], [28]) is a straightforward approach to obtain a low-dimensional uncorrelated process from a correlated high dimensional or even infinite-dimensional process. Holmes et al. [28] examined the theoretical background of the POD and its properties profoundly. In the following sections, the mathematical basics of the POD are discussed shortly, but as the paper is more targeted to the strategic approach in earthquake engineering, the attention to the mathematical background and the numerical problems are limited to an essential minimum.
The aim of the POD is to find a set of ordered orthonormal basis vectors in a subspace so that samples in a sample space are expanded in terms of l basis vectors in an optimal form. This means that the POD is able to find an orthonormal basis, which describes an observation vector in a subspace better than any other orthonormal basis can do. A measure for this problem is the mean square error (cf. Qu [29]) where x 2 R n 1 is the random vector, x.l/ is the approximation of this random vector in an ldimensional POD subspace, and O x.l/ is the approximation of the random vector by any other possible orthonormal basis. Therefore, the random vector can be expressed as (cf. Qu [29]) 541 2 .l; t/ D E ® kx x.l/k 2¯! mi n (8) subject to the orthonormality condition (cf. Qu [29]) ' T p;i ' p;j D ı ij .i; j D 1; 2; : : : ; s/ : The transformation into the l-dimensional POD subspace is a truncation of the first l lower POD modes (cf. Qu [29]) x.l/ ˆpq p ;ˆp D OE' p;1 ; ' p;2 ; : : : ; ' p;l ; l < s n:

THE NEW APPROACH
In this section, the reader's attention is drawn to the methodical approach. Therefore, a simple and illustrative two-dimensional structural example subjected to a set of ground motions goes along with the presentation of the new MOR strategy. Earthquake excitation records are presented in Section 2, respectively in Table I. The academic example concerns a two-story frame system as an academic example modeled by nonlinear elasto-plastic beam elements. The geometrical discretization is shown in Figure 1. The frame span is l D 6OEm, and the height of it is h D 4OEm; consequently, the total height is eight 8 m. The columns as well as the beams of the frame structure are discretized by five elements, which leads under consideration of the boundary conditions shown in Figure 1, to a total number of 86 DOFs. Elasto-pasticity is realized by a bilinear stress-strain curve with kinematic hardening in axial beam direction as shown in Figure 2. Fictitious material parameters are assumed in this academic example to induce large plastic deformations in order to illustrate a clear visualization of the nonlinear behavior and accordingly a good insight into the nonlinear MOR procedure. The paramters of the elasto-plastic  In addition to the presented nonlinear system with hysteretic material behavior, the equivalent structure with linear material parameters is presented. The linearized material parameter is equal to the initial stiffness of the elasto-plastic material model presented in Figure 2, that is, Young's modulus E D 2: 1 10 11 h N m 2 i . The extended approach of the new MOR strategy is based on the mathematical formulations of the POD in Equations (6)- (10). In structural dynamics, systems are discretized in space and time; therefore, the numerical realization of the proposed strategy is implemented based on the following straightforward algorithm.
Firstly, an a priori response identification, that is, the realization of a random vector x.t/, is evaluated, where t is in a limited time interval t 0 6 t 6 t 1 . Therefore, the response to one excitation in the specific time period is computed, that is, Equation (1) is solved numerically in the defined time period. This analysis is performed in the physical coordinate x. Consequently, this constitutes dependent on the size of the time period t 1 t 0 (snapshot time period) and the number of DOFs of the system, the most time consuming part of the procedure. According to the actual example, for the linear as well as the nonlinear system, the Bam earthquake in Table I is chosen as the representative event. The snapshot time period is here limited by t 1 D 0 OEs and t 2 D 12 OEs. This means that the largest possible time window for the snapshot response identification is chosen; therefore, the maximum possible response information is captured here. Although in order to increase effectivity of the algorithm, a representative deterministic basis can be obtained by concentrating on the strong-motion phase of the earthquake because the relevant nonlinearities will be activated during that phase, but then it is more likely that insufficient POD response bases for MOR transformation matrix for the upcoming analyses are acquired. The time history of the excitation is shown in the left subplot of Figure 3.
The snapshot response function x s .t/, that is, the random vector, is realized by s observations (snapshots) at different time instances (cf. Han and Feeny [30]) within the snapshot time period X s D OEx t 1 ; x t 2 ; : : : ; x t s D 0 @ It is of utmost importance to capture nonlinear deformation patterns in order to create the possibility to depict possible nonlinear reactions in forthcoming calculations. In the specific example, 40 observations (snapshots) in equidistant time instances are taken into account from the snapshot response function for the a priori identification of the overall response behavior. Those observations of the nonlinear system are specified by dots in the right subplot of Figure 3, to be distinguished from the overall response within the time period of analysis. According to the linear system, the computation of the snapshots is realized by the same procedure. As a result, it is also necessary here to capture the main motion patterns according to possible future response histories. However, it can be observed that for the system identification of the linear response, the required number of snapshots is considerably smaller than for the nonlinear system response.

543
If is the expectation of all observations, then the sample covariance matrix † s of this random vector, which is realized by the observation matrix, is defined by (cf. Kerschen et al. [31] The POD modes and the POD values are defined by the eigensolution of the sample covariance matrix. If the data have zero mean, the covariance matrix is (cf. Kerschen et al. and the POD is realized by the SVD of the observation matrix X s . The POD modes ' p;i are equal to the left singular vectors and the POD values p;i to the singular values of X s , which are all real and positive and arranged in a rectangular diagonal matrix in descending order. The energy, which is contained by the snapshot matrix, is defined by the summation of the POD values, that is, V D P s i D0 p;i . As a consequence, the energy ratio of the i th POD mode is (cf. Kerschen et al. [31]) In structural dynamics applications, the sum of only a few POD values often captures the main part of the total energy included in the observation matrix, which reflects the big advantage of the POD, that is, the property of optimality with respect to energy in a least square sense. In the present demonstration, 99% of the total energy is captured by only two deterministic modes. The first two evaluated POD modes of the nonlinear system are depicted in Figure 4. In comparison, the natural modes of vibration, which create the basis for the classical modal truncation method are shown in Figure 5. They are computed according to the eigenvalue problem K ! 2 M ˆD 0, where K is the linearized initial stiffness matrix according to the aforementioned assumption. As depicted in Figure 4, the POD modes can obviously represent nonlinear displacement patterns, which are expected to occur in this illustrative example, for example, plastic deformation areas around the joints of the frame system. But rather these deformation patterns are not representable by the first two linear natural modes of vibration depicted in Figure 5. Therefore, according to the nonlinear  response, an accurate representation of the snapshot response function x s .t/ is expected by the representation of only the two deterministic POD modes, but a significant error by the representation of the first two natural modes of vibration is expected. Following, the transformation into the reduced subspace is performed in the same manner, as if the classical method of modal truncation would be applied to linear systems. The low-order set of equations of motion is then where Q M DˆT P MˆP and Q C DˆT P CˆP are mass and damping matrices, and Q F DˆPF is the excitation vector in the POD reduced subspace. The reduced vector of the inner restoring forces Q R is still dependent on the displacement in the physical coordinate x, Therefore, the vector of the inner restoring forces R.x/ has to be evaluated from the physical model in the full-order coordinates in every calculation time step. According to the linear POD reduced system, the vector of the inner restoring forces is described by the relation Q R D Q Kq P , where Q K DˆT P Kˆp is the stiffness matrix in the POD reduced subspace. In this case, no updating in every time step of the inner restoring forces has to be performed. Furthermore, the equations of motion in the reduced-order set are not decoupled and have to be solved numerically. Finally, after the time integration is finished, the solution vector, q P , dependent on time is transformed back into the physical coordinate x. Concerning the representative excitation, a comparison of the full benchmark solution (snapshot response function x s .t/) makes sense in order to verify the quality of the low-dimensional representation of the POD modes. If the approximation of the benchmark solution (in this example, Bam earthquake response) is satisfactory, then the reduced solution to the whole set of the remaining earthquakes is calculated by application of these deterministic modes. If the approximation is not sufficient, one possibility is to increase the number of snapshots within the snapshot time period; another one is to change the start and end time instances t 0 and t 1 of the snapshot response function in order to improve the response information quality of the snapshot matrix (remark: in this paper t 0 and t 1 are already the beginning and end time step of the Bam excitation record). The visualization of the novel approach is depicted in Figure 6.
All computational results according to the whole excitation set presented in Table 2 are depicted in Figures 7-12. Within each of the response computations of the linear, as well as the nonlinear system, two different integration algorithms were chosen for the integration procedure in the physical (full) coordinate (86 DOFs), that is, the implicit Newmark integration (NEW) and the explicit central difference scheme (CD). They serve as reference solutions to the suggested new POD strategy, which is compared with the classical modal truncation method (MT). All time response functions are compared according to the horizontal displacement in the right corner of the first floor of the frame structure,  According to the proposed strategy in Figure 6, the Bam earthquake serves as the representative event ( Figure 3). A test run of the reduced order model over the time period of the representative Bam earthquake shows that the deterministic POD modes can represent the nonlinear (elasto-plastic) response time history. This is shown in the top right and the bottom right subplot in Figure 7. According to Figure 7, both the linear and the nonlinear POD response functions approximate the full reference solutions (NEW and CD) accurately. Therefore, the deterministic set of POD modes is able to represent nonlinear (plastic) deformation patterns. Consequently, the chosen number of 40 snapshots within the chosen snapshot time period (Figure 3) is apparently sufficient for this example. Figures 8-12 present the numerical evaluation of the rest of the earthquake set shown in Table 2, that is, the comparison of the reduced order models with the full benchmark solutions. In each of the linear and nonlinear response evaluations, it is shown that the POD responses provide reliable approximations of the benchmark solutions (response functions are overlapping with each other). This underlines the robustness of the new strategy according to changes in the excitation if the deterministic POD modes are evaluated based on a different excitation history ( Figure 3).
According to the linear computations, both strategies show accurate results (see the top right subplots of . In this special linear elastic case, the method of modal truncation shows slight advantages, which is not surprising as no snapshot response function in the physical coordinate must be evaluated. However, a distinct advantage is that the POD strategy does not require an investigation of the Fourier transform of the excitation in order to determine the cut-off frequency for estimating the truncation of the natural basis modes. On the contrary, the POD method fulfills the requirement in Equations (6) and (7), where an optimally truncated basis is defined. Therefore, dependent on the required level of accuracy, an automatic truncation based on the energy content per mode is performed (Equation (14)). According to the nonlinear system (see the bottom left subplots of Figures 7-12), the modal truncation responses show a severe deficit with respect to the exact solution, because the natural modes (as earlier presented within the Figures 4 and 5) cannot detect the nonlinear deformation patterns. On the contrary, the POD strategy, as already discussed previously, shows reliable approximations. Therefore, in the presence of strong nonlinearities, the new POD strategy has sizeable advantages over the classical method of modal truncation. However, within a limited range, the modal truncation method is reliable. In the top right subplots of Figures 7-12, the separation point (sep. point, black mark) defines the time instant, when, according to Figure 2, the first nonlinear effect occurs (yielding of the material). Before this time instant, each response line (CD, NEW, POD, and MT) is the same. This is not surprising, because this is the simple linear-elasic solution. Afterwards, as presented in the top right subplots of Figures 7-12, the nonlinear yielding effect is shown by the drift-off of the light gray nonlinear reference responses in these subplots (nonl. resp.), which is equivalent to the nonlinear Newmark response of the bottom right subplots of these figures. Here, the black mark also shows the separation point of the linear and the nonlinear response, and additionally, the qualitative separation point of the modal truncation solution (blue mark), that is, where the modal truncation solution drifts off. Therefore, it is clearly shown that there is a limited range of nonlinearity, where the modal truncation response also presents useful results -after the modal truncation separation point, a severe drift-off of this solution function is observed.
The first substantial numerical benefit of the proposed strategy arises from the combination of the MOR and the application of an explicit time integration method, such as the second order central difference scheme. The time consuming part of the computation is the evaluation of the inner restoring force vector, which has to be evaluated in every time step on Gauss integration point level in the physical (full) coordinate according to Equation (16). This recalculation cannot be avoided. The application of an explicit time integrator to the full system leads inevitably to considerably small integration time steps, which have to be smaller than a certain value (critical time step for the central difference scheme t cr D 2 ! n , where ! n is the highest eigenfrequency). Therefore, a huge number of time steps has to be processed within one integration run and consequently, a high number of the expensive recalculations of the inner restoring force vector R.x/. For the integration in the POD reduced subspace the critical time step is considerably larger (according to this example in the paper even larger than the measurement time step of the earthquake data), which leads to a stabilization of the procedure and as a consequence, to approximately a 500-fold increase of speed compared with the full central difference integration method (CD). These numerical issues are discussed in a similar manner by Gutierrez and Zaldivar [19] applied to modal truncation. For the numerical benefit of the combination of the basic POD method with explicit numeric time integration, the reader is referred to Bamer and Bucher [21]. It would also be possible to implement an implicit integration scheme (i.e., the Newmark method with constant acceleration assumption) in the reduced subspace, which is under special circumstances convergent for all possible time steps. However, the application of this integration scheme requires the iterative Newton Raphson procedure, following, in each timestep of the reduced system, the restoring force vector R must be computed several times (modified Newton Raphson method). If the classical (not modified) Newton Raphson method is applied, the number of iterations is minimized, but then in every iteration procedure the tangential stiffness matrix in the full space is to be evaluated, which has then again transformed into a tangential stiffness matrix in the reduced subspace. Therefore, the explicit central difference scheme is chosen, which is stabelized by the model order reduction process. Generally, each integration scheme can be applied for the realization of the POD reduced strategy, and of course, the choice of the integration scheme is also dependent on the level of high dimensionality, type of nonlinearity, and so on.
The second big numerical advantage of the new strategy is that the actual time consuming process, which is the evaluation of the snapshot matrix, is only executed once at the beginning of the whole calculation procedure. This a priori assumption of nonlinear mode patterns makes sense if the excitations show physical 'similarities', which is the case in earthquake analysis, where a considerably small number of low frequency modes is mainly affected.
The practical application of this novel approach in earthquake engineering analyses can be realized in a straightforward manner. Earthquake design codes frequently specify the earthquake loading in terms of response spectra rather than actual acceleration records, but this is directly usable only in the case of linear system behavior. For the analysis of nonlinear systems, several response spectrum compatible artificial earthquake accelerations (as specified, e.g., in Eurocode 8) can be generated, see for example, Ref. [32]. Based on these artificial accelerations, the proposed model order reduction approach is carried out exactly as described in this section.

Error evaluation
For the sake of accuracy evaluation a posteriori error analysis is presented. As reference solution to the reduced response functions, that is, solutions by POD reduction and modal truncation, and the solution of the central difference algorithm in the physical coordinate is chosen. The error indicators are  Figure 1. Figure 13 represents the absolute error of the reduction methods corresponding to analyses provided in Figures 7-12 as defined in Equation (17). As depicted in Figure 13, the error by application of the classical modal truncation method is much more considerable than the error produced by the approximation through the POD strategy. Especially, if the structure responds nonlinearly, that is, plastic deformation patterns occur, a huge difference by application of the modal truncation method is observed. Figure 13 provides evidence of the inability of the linear modes to represent elasto-plastic displacement patterns shown in Figures 4 and 5. Additionally, the applicability and accuracy of the proposed POD strategy seem not to be remarkably sensitive to differences of the peak acceleration response spectra of the excitation time histories in the presented examples. This is seen in the peak acceleration response spectra in the bottom left plot of Figures 7-12 and the related error functions in Figure 13, where a correlation between those parameters cannot be detected directly. One noticeable point here is that the error functions of the Landers and the Northridge response in Figure 13 seem to be in similar ranges, but the characteristics of the response spectra show considerable differences. Therefore, the proposed strategy demonstrates a high robustness with regard to time history analysis in earthquake engineering and should be of great value in this field. However, a more extensive evaluation of direct correlations between changes in the spectra of the excitation histories and the quality of the outcome of the POD response is beyond the scope of this paper. Additionally, future research should also include a priori error estimations based on wavelet transformations of the excitation functions. The underlying idea is that earthquake records with similar intensities and frequency contents activate the same nonlinearities and consequently lead to the same POD basis. A publication concerning this strategy appeared quite recently [33].

PRACTICAL APPLICATION OF THE NEW APPROACH
In addition to dealing with the development of the introduced novel POD-based MOR approach, it is within this section to represent the application of the new proposed MOR strategy on a realistic example. For this purpose, a dynamic structural model of a medical complex, according to its constructional plan, was derived. A schematic three-dimensional sketch of the building is depicted in Figure 14.
As shown in Figure 14, the building structure exhibits complex geometries. As a result, it seems to make sense to discretize the geometry by a FE model in order to capture the main dynamic specifications.
If such a structure with medical function is located in an earthquake prone region, one way to improve its seismic performance can be realized through base isolation by means of frictional pendulum bearings. Consequently, the analytical simulations demand large computational time and storage because of the presence of nonlinearity imposed by those frictional isolators. In the following, firstly, the structural system specifications and implementation of frictional bearings are presented. Then, the displacement responses to the set of the six earthquake events presented in Table 2 are evaluated. The numerical evaluations compare the new introduced strategy, as an alternative means, with the iterative Newmark integration scheme, which is known as an efficient and exact method.

Structural system and model specifications
The building structure consists of three wings, referred to as wing I, II, and III. Figure 15 shows a schematic sketch of the ground plan of the building containing the basic dimensions. The floor slabs of each wing are separate from the others except for the basement slab, which is indiscrete over all three wings. This means that all three wings are coupled through this slab, and they work all together during the earthquake excitations. However, the distance between the wings, which are connected by the basement slab, is about 1.5 m; consequently, contact problems induced by ground motion are not considered in the computations. The grid indicates the location of the columns and the binding beams, and the red lines indicate the location of the shear walls, which are responsible for the lateral reinforcement. The regular distance between the columns is 6.5 m. The building structure has three stories below the ground level, while the highest parts of the building above ground level have 13 stories and the remaining parts have eight stories including the basement levels. Therefore, the plan of the structure is irregular along its height along with the irregularities in the horizontal area. The dashed lines define the area, where the building is only located below the ground. The height of one story is 3m; this leads to a total construction height of 42 m.
Below the three stories at the basement level, there is the indiscrete slab on the top of the isolators at level of 9.00 [m]. Below this slab, along each of the columns, a single friction pendulum (FP) bearing system is attached. Figure 16 depicts a part of the cross section A-A of the basement level shown in Figure 15. The horizontal diameter of the FP system is 2.00 m. Thus, the dimension of the quadratic cross section of the columns in the basement and FP story is 2:00 2:00 OEm 2 , while in the remaining stories the columns are modeled as quadratic cross sections with the dimensions 0:40 0:40 OEm 2 . All FP bearings have the same radius of the concave surface, which is equal to 3.00 m.   A representative full-scale FE model of the building structure was created in the software package slangTNG [34]. The shear walls and slabs were modeled by shell/plate elements and the columns and beams by beam elements. A linear elastic material was considered for the modeling of the superstructure (Young's modulus E D 3:5 10 10 h ). Nonlinear FP elements, whose implementation in slangTNG is presented in Section 5.2, are assigned below the lowest basement plate of the structure. The total number of DOFs is 33000.
As a result, the superstructure behaves linearly, which is the major reason for implementing base isolator systems for earthquake vibration protection.

Dynamic model of the frictional pendulum element
This is to present how the FP element in the finite element model of the structure behaves. The geometrical diagram of the FP element, which is realized as a spherical shell, is defined in Figure 17. As depicted, R denotes the radius of the concave spherical surface, and the origin of the local coordinate system is chosen to be in the center of the sphere. The position vector of the slider is described by U D OEu; v; w T . Since the desired behavior of the FP element is an in-plane elasto-plastic bidirectional action, the change of the vertical position w can be neglected. Accordingly, the displacement of the FP element is reduced to an in-plane motion defined only by the components u and v, that is, U D OEu; v T . This simplification makes sense as the radius R is much larger relative to the horizontal displacement 552 F. BAMER, A. KAZEMI AMIRI AND C. BUCHER Figure 18. Internal specifications of the friction pendulum element; F x , F y , N , recentering force F k , frictional force F Fr . [Colour figure can be viewed at wileyonlinelibrary.com] jUj p u 2 C v 2 . The equivalent representation of such an element together with the acting forces on it is represented in Figure 18.
The horizontal force equilibrium of the dynamical system is where F Fr and F k are the elasto-plastic frictional and centring force, and F ex D OEF x ; F y T accounts for the interacting horizontal force, which couples the FP element to the super structure. Following, the force equilibrium is split into two parts as two dynamic situations can occur: situation stick and situation slide. The force equilibrium during the situation stick yields to This relation renders a linearly-elastic system, where the friction coefficient must be a value between 0 and 1. For the computations, this value was taken to be once 4%, and for the second demonstration run of the method it was set equal to 8%. The normal contact force N acts orthogonal to the contact area of the slider and the concave surface. The vector U D OEu; v T defines the radial distance with respect to the current sticking point of the slider if the sticking condition is true. During the situation slide, the FP element is described by the following horizontal force equilibrium: where P U D OE P u; P v T is the velocity vector. In both relations, that is, Equations (19) and (20), the centring force, jF k j D k 1 r D k 1 p u 2 C v 2 , acts linearly orthogonal to the vertical axis through the deepest point of the surface and the center of the sphere. The fact that the centring force is linear indicates that the spherical shell of the real system is approximated by the paraboloid, whose potential energy increases with W R r in radial distance from the deepest point, that is, the stiffness is inversely proportional to the radius of the sphere k 1 D W R . The frictional force F Fr is modeled either linearly elastic or elastic-perfectly plastic as presented in Equations (19) and (20), respectively. Note that the force corresponding to U accounts for the elastic behavior of the bearing coating material, in a small elastic range (situation stick, Equation (19)) and acts towards the current location of the slider (not the center of the concave sphere). Generally, the implementation of a realistic model requires k 2 to be much larger than k 1 , that is, k 2 k 1 . During the situation slide, the frictional force acts in opposite direction to the velocity with the magnitude (perfectly plastic) N . This is discussed in Equation (20). Another point regarding Equation (19) is that the reacting force N is assumed to be constant throughout the calculation procedure. This is justified by the following reasons: firstly, just x and y components of the exciting ground motion are taken into account for the computations. Secondly, the motion has already been simplified to be planar and therefore no additional force component due to vertical motion is generated. Finally, in our preliminary analysis, the uplift force on the isolator slap was observed to be extremely small in comparison with the downward force because of the weight of the structure. Considering the abovementioned fact together with the force diagram given in in Figure 18, follows that the normal contact 553 Figure 19. Proper orthogonal decomposition values; friction coefficient D 0:04 OE (logarithmic scale).
force N is approximately constant and equal to the weight induced force of the super structure, W , that is, N D W in Equation (20).
The FP bearing element governed by Equations (18)- (20) has been implemented in the software package slangTNG [34]. For a comparable study on this implemented friction pendulum system, the experimental work of Mosqueda et al. [35] is suggested. Concerning the computations in this paper, the friction coefficient of the FP slider is assumed to be a constant value, as discussed previously. However, is indeed a value dependent on velocity, pressure, and temperature. Recent publications dealing with this issue are authored by Castaldo and Tubaldi [36], as well as Kumar et al. [37]. Regarding the proposed new strategy, it is expected that the implementation of a friction coefficient, which depends on velocity, pressure, and temperature, would lead to comparable results concerning accuracy and speed of the reduced order model. In other words, it is envisiged that a specific nonlinear model of the FP bearing element with variable friction coefficient does not influence the overall behavior of the reduced order model in comparison with the full model, that is, the effectiveness of the model order reduction procedure. In every time step, the vector of the inner restoring forces in Equation (16) has to be evaluated in the full (physical) coordinate anyway.
For additional information about friction pendulum systems, the reader is referred to the literature (e.g., [38][39][40][41][42][43]). More detailed examination of this topic would lead beyond the scope of this paper, which should focus more on the methodical extension of the new MOR strategy as well as the application on a complex realistic system.

Numerical evaluation
The evaluation of the MOR strategy applied to the realistic building structure is dealt with displacement response calculations to the six different earthquake excitations presented in Table I. The calculation outputs are presented by the in-plane motion response of the slider at the red marked node (node 1: coordinates OE19:5; 0:0; 9:5 T OEm) in Figure 15 in x and y direction. As this node defines the location of a moving friction pendulum, it shows directly the nonlinear response behavior of the system. Additionally, a second output node (node 2: coordinates OE32:5; 0:0; h T OEm) is chosen. Here, the coordinate h stands for every possible story of the building structure. According to h D 32:5 OEm; the maximum acceleration of the roof is presented and according to each story .h D 9:5; 6:5; 3:5; : : : ; 32:5 OEm/, the maximum drift responses are shown.
The numerical demonstrations are performed by assuming two different friction coefficients for the FP-isolators, namely 4% and 8%. The methodology applied to the realistic example is equivalent to the academic example, which is presented in Section 4. Therefore, the main focus in the next sections is on the presentation of the numerical results.
Integration over the whole Bam earthquake by the Newmark method in the physical coordinate leads to the snapshot response function x S .t/ and consequently to the snapshot matrix X S , which contains 400 snapshots in equidistant time intervals. Following, the evaluation of the left singular vectors of the snapshot matrix leads to the POD modes and its singular values to the POD values in descending order. The number of POD modes that have to be taken into account in order to capture 99.99% of the total energy was evaluated to be 31   the reduced response must be compared with the benchmark solution (full Newmark response). The response motion of the slider (output node 1) of the full and the POD-reduced calculation is depicted in Figures 21(a) (friction coefficient of 4%) and 22(a) (friction coefficient of 8%). As clearly shown, a reliable approximation of the benchmark solution is achieved; therefore, nonlinear effects can be represented. In the next step, the responses to the remaining earthquakes are evaluated in the reduced subspaces, as well as the full reference solutions by application of the standard Newmark algorithm. The first 30 modes are computed according to a linearization of the slider element. The boundary conditions below each slider element are fixed. The required linearization of the slider element is realized by taking the initial stiffness matrix of the system into account, similarly to the linearization assumption of the academic example. This inevitably leads to an over-estimation of the slider stiffness, although the vector of the inner restoring forces is calculated.
The red lines in Figures 21 and 22 show the POD reduced responses and the blue lines show the full reference solutions obtained by the standard Newmark method. As depicted in these figures, accurate approximations are achieved for both assumptions of the friction coefficients (note: the lines cover each other). Additionally, it is clearly shown in these figures that the modal truncation method (green lines) fails to approximate the full solution.
The maximum accelerations at the roof in output node 1 are presented in the left and right subplot of Figure 23. Additionally, the maximum drifts between all floors above output node 2 are presented in Figures 24 and 25. As for both friction coefficients, the slider drift itself (drift in the floor 0) is generally much larger (which is actually the purpose of base isolation). Therefore, for each of the It is expected that the proposed new strategy is also applicable to different types of nonlinearities, for example, large deformations, viscoelasticity, viscoplasticity, hyperelasticity, and so on. In addition, concerning seismic protection, it is expected that the proposed new strategy is also applicable to  smart structures with shape memory alloy-based seismic damping and isolation tools. The new strategy should produce response approximations with comparable accuracy to the nonlinearites already presented in this paper. However, it is also of utmost importance here to emphasize the fact that the nonlinear deformation patterns, namely the special hysteretic behavior dependent on superelasticity, temperature, and memory effect, must be captured in the snapshot matrix.

CONCLUSION
In this paper, a MOR strategy, which is applicable to the dynamic response analysis of linear and nonlinear structural systems is presented. Usually, the analysis of building structures with complex geometries makes the engineer to create a FE model with a large number of DOFs, which is associated with computational effort in the response analysis. Therefore, the goal of this paper is to provide a new MOR strategy that is simple in application, but also very effective even in the presence of nonlinearities for problems in the field of earthquake engineering and structural dynamics. This strategy is extended based on the POD method to derive a proper transformation matrix in order to transform the nonlinear systems into another low-dimensional subspace, which demands considerably less computational effort for the response calculation. Once the transformation matrix is derived, the approach of the strategy is similar to the method of modal truncation for linear systems. The presentation of the novel approach comes along with a simple and illustrative example that points out the benefit compared with standard methods as modal truncation.
In addition to the development of the MOR strategy, its application for the response calculation of a realistic numerical nonlinear example is demonstrated. The example is the displacement response calculation of a building structure serving as a medical complex, which is base-isolated by FP bearing systems excited by six earthquake excitations. In order to evaluate the accuracy of the introduced approach, the exact structural responses are also calculated. Numerical evaluations show that reliable approximations can be achieved if nonlinear response patterns of the structure are already captured in the POD snapshots to extract the transformation matrix. The advantage of this strategy is obvious that the transformation matrix is derived just once, and it can be used for response calculation of the structure under different earthquake excitations.
Another substantial advantage of the introduced MOR concerns the speed of the response calculations. Firstly, compared with the basic central difference algorithm, the new introduced strategy has a much larger critical time step. Secondly, compared with the Newmark method, which allows usually larger time steps, no iteration procedure is required.