Abstract
Based on the principle of real-time substructure shaking table test, an interactive numerical computation method which built the calculation model of each substructure in different software programs was proposed for seismic analysis of equipment–high rise structure–soil systems in order to account the effect of nonlinear soil. Considering that the response of soil under strong earthquakes does not totally enter the nonlinear stage, a locally nonlinear soil model was introduced as the numerical substructure, and the equipment–high rise structure subsystem was treated as an experimental substructure in this method. The equation of motion for the equipment–high rise structure–soil system was derived through a combination of the branch modal substructure and linear–nonlinear hybrid constraint modal substructure approaches. A 13-layer steel framework system model is used as an analytical example that the equipment–high rise structure system and local nonlinear soil computing model are built by MATLAB and ANSYS, respectively. The time histories of the system dynamic responses were obtained by interactive numerical computation, to investigate the effects of equipment–high rise structure–soil interaction on the seismic performance of the equipment and structure.
Keywords
Introduction
With the rapid economic growth, large numbers of high rise structures (HRSs) have been built in the world. Besides, the proportion of equipment and other nonstructural components in buildings increases. In this context, the soil–structure interaction (SSI) becomes an increasingly important issue in seismic analysis and design of HRSs. Research revealed that the total lateral displacement of an HRS tended to increase when the behavior of soil was taken into account and that nonlinear soil properties can increase or decrease the base shear of structures compared to linear soil response, depending largely on the type of structure, seismic input frequency, and soil type.1,2 This suggests that SSI may have an adverse effect on the structure’s response in some cases. There are many computational and experimental methods that take SSI into consideration, including finite element calculation,3,4 on-site model excitation testing, 5 centrifuge testing, 6 and shaking table testing.7–12 Shaking table testing, currently one of the most widely used methods, is capable of effective simulation real earthquake excitations. In this method, the shaking table platform applies seismic excitations to the structural models. To take into account SSI in shaking table tests, a soil-filled box model is usually used to represent on-site soil, and the structural system is then embedded in the ground for testing. This method requires a construction of complete test models and can be used to study the seismic responses of a structure and soil simultaneously. However, the modeling process is usually costly and technically difficult, and the size of shaking table and other limitations prevent large-scale model tests. With the development of real-time substructure testing and control techniques,13,14 the real-time substructure shaking table test (RTSSTT) testing method has been applied to the seismic field. And the RTSSTT that includes the effect of soil15–17 is gaining much attention. It is a huge challenge that nonlinear soil model as the numerical substructure in RTSSTT. Therefore, there is an urgent need for a more efficient and accurate computational method that considers the effect of nonlinear soil.
Considering that the response of soil under strong earthquakes does not totally enter the nonlinear stage, Wang and Jiang 18 systematically analyzed SSI using the combination of the branch modal substructure and linear–nonlinear hybrid constraint modal substructure approaches and then validated the proposed computational method by a case study. The advantages of linear–nonlinear hybrid constraint modal substructure method were also demonstrated in the dynamic analysis of three-dimensional models. 19
Based on the principle of RTSSTT, an interactive numerical computation method was proposed to achieve a simulated process of RTSSTT. This method used the linear–nonlinear hybrid constraint modal substructure approach to construct a locally nonlinear soil (LNS) model considering that the soil is not completely in nonlinear status under strong earthquakes. Then the equation of motion for an equipment–HRS–LNS system was derived using a combination of the branch modal substructure and linear–nonlinear hybrid constraint modal substructure approaches. Then this equation was transformed to obtain the equations of motion for the equipment–HRS subsystem and soil substructure for use in RTSSTT. A computational model of an equipment–HRS–soil system was selected. In this model, the HRS was a 13-story, 4-span steel frame structure, the equipment was simulated with a simplified model with a single degree of freedom (DOF), and the soil was Class III ground. The equipment–HRS subsystem was constructed using the programming software called MATLAB, and the LNS model was created using ANSYS. A comparison of the results of interactive numerical computation and integrated finite element analysis demonstrates that the analysis method proposed in this study is reliable. After that, the time histories of the equipment–HRS–soil system dynamic response were obtained by interactive numerical computation to investigate the effects of equipment–HRS–soil interaction on the seismic performance of the equipment and structure.
Derivation of the equation of motion for the equipment–HRS–LNS system
In this study, the linear–nonlinear hybrid constraint modal substructure approach was incorporated to create an LNS model, with an attempt to overcome problems such as redundancy of DOFs of soil.
The equipment–HRS–soil system shown in Figure 1 consists of equipment, HRS, and a soil. Branch modal substructure method can be employed to derive the equation of motion for the equipment–HRS–nonlinear soil substructure β system, from which the coupling terms between the equipment–HRS subsystem and nonlinear soil substructure β can be obtained. And the equipment–HRS–nonlinear soil substructure β system was divided into branch d (equipment–HRS rigid system on nonlinear soil substructure β), which is shown in Figure 1(b), and branch s (equipment–HRS deformed system on rigid nonlinear soil substructure β), which is shown in Figure 1(c). Unlike common branch modal substructure approaches, the computational matrices for the equipment–HRS subsystem and nonlinear soil substructure β were totally preserved in this study to consider their nonlinearity separately. Later, the number of DOFs of the linear soil substructure α was reduced through the constraint modal substructure method. After coordinate transformations based on the deformation compatibility at the interface between the linear soil substructure α and nonlinear soil substructure β, the equation of motion was derived for the equipment–HRS–LNS system. The derivation process is detailed below.

Derivation of the equation of motion for the equipment–HRS–LNS system. (a) Equipment–HRS–soil system. (b) Equipment–HRS rigid system on nonlinear soil substructure β. (c) Equipment–HRS deformed system on rigid nonlinear soil substructure β. (d) Reduce DOFs of linear substructure α.
As the nonlinearity of both the equipment–HRS subsystem and nonlinear soil substructure β during a strong earthquake must be considered, the computational matrices for them were preserved. The equipment–HRS subsystem displacement, us, can be divided into two components: the rigid body displacement caused by the motion of the nonlinear soil substructure β (Figure 1(b)) and displacement due to the motion of the subsystem itself, qs (Figure 1(c)). us is calculated using equation (1)
R in the equation can be determined from the equipment–HRS subsystem’s rigid body displacement produced by the deformation of nonlinear soil substructure β. The equation governing the relationship between the displacements of the equipment–HRS subsystem and nonlinear soil substructure β can be described in matrix form given by equation (2)
According to the principle of inertia coupling in a branch modal substructure method, rearranging the equation above gives the equation of motion for the equipment–HRS–nonlinear soil substructure β system (equation (3)). In the resultant equation, only the mass matrix contains coupling terms, while the stiffness and damping matrices are decoupled
where ms, cs, ks, and fs represent the mass matrix, damping matrix, stiffness matrix, and load matrix, respectively, for the equipment–HRS subsystem in branch s; mβ, cβ, kβ, and fβ denote the mass, damping, stiffness, and load matrices, respectively, for the nonlinear soil substructure β in branch d. The computational matrices for the nonlinear soil substructure β were then repartitioned based on the types of internal DOFs i and boundary DOFs b, yielding a new equation of motion (equation (4)). In the resultant equation,
The computational matrices for linear soil substructure α were also rearranged based on the types of internal DOFs i and boundary DOFs b. After that, the modal reduction was performed through the constraint substructure method. The corresponding coordinate transformation matrix T is given by equation (5), which consists of the dominant mode
After the coordinate transformation based on the constraint modal substructure method, the equation of motion for the linear soil substructure α has the following form (equation (6))
The interface between the nonlinear soil substructure β and linear soil substructure α must satisfy the following compatibility condition (equation (7))
The equation of motion for the equipment–HRS–soil system was finally derived by combining the equations of motion for the linear and nonlinear soil substructures based on the condition for displacement compatibility. It is given by equation (8)
where msRi and msRb are two coupling terms in the mass matrix, representing connections between the equipment–HRS subsystem and LNS. To enable substructure testing, loads associated with the coupling terms were moved to the right-hand side of the equation, giving the equations of motion for the equipment–HRS subsystem and the LNS (equations (9) and (10))
Then the interactive numerical computation process for the equipment–HRS–soil system was inferred from the obtained equations of motion. It involves the following steps: (1) assume that the external load and the loads associated with the coupling terms acting on the LNS at the ith step are known; (2) calculate the soil’s acceleration response at (i + Δt)th step; (3) calculate the external load and the loads associated with the coupling terms to be applied to the equipment–HRS subsystem at this step and then transmit them to it; (3) obtain the acceleration response of the equipment–HRS subsystem and the force on the LNS produced by the seismic excitation and then transmit the force to the soil for the next numerical computation. In this way, data communication between the equipment–HRS subsystem and LNS can be achieved at every step and continue until the process ends.
And the equipment–HRS subsystem model and LNS model could be created using MATLAB and ANSYS, respectively, and then they participated in the interactive numerical computation described above.
Since every step of the computational process involves data transmission, it is necessary to create a text document in each simulation software working directory to store the data transmitted. This can help the simulation software easily read the computational data during the whole computational process. For example, MATLAB uses the command loop to continuously read data in its text document, where the first line of numbers will tell whether the transmission of the computational data required is completed. If the first number is “0,” ANSYS is still computing, and thus MATLAB will continue waiting and reading the text document until the number in the first line becomes “1,” which indicates that ANSYS completed computation. After reading the data, MATLAB starts computation. The text document in ANSYS working directory works on a similar principle.
Case study
A case study of a 13-story, 4-span equipment–HRS–soil system was performed. In this system, the structure was a 13-story, 4-span steel framed building, and the equipment, a communication device placed on the floor of the top story, was simulated with a single DOF model. The soil was Class III ground.

Acceleration time history curves of ground motions.
Three ground motions that are suitable for Class III ground, including El Centro, TianJin, and PerSon ground motions, were used as the seismic inputs (Figure 2). Table 2 shows the acceleration time history curves of ground motions. Then time histories of the system’s dynamic responses to the ground motions of level VIII intensity during the minor earthquake, moderate earthquake, and major earthquake stages were calculated with the acceleration amplitude adjusted to 0.7, 2.0, and 4.0 m/s2, respectively. The equipment–HRS subsystem model was constructed by MATLAB programming, and the LNS model was created with ANSYS through the linear–nonlinear hybrid constraint modal substructure approach. Then both models were used in interactive numerical computations.
Equipment–HRS subsystem model
The computational model of the equipment–HRS subsystem is shown in Figure 3. The structure model was constructed of H-beam, and the cross section of each beam or column is shown in Table 1. Each story of the structure was 3.56 m high, longitudinal and transverse spans were 7.32 m. The total mass of each story was 5.86 × 105 kg, except for the 13th story which weighed 5.32 × 105 kg. All the beams and columns had an elastic modulus of 2.01 × 1011 Pa. The yield strengths of each column and beam were 345 and 245 MPa, respectively. The damping ratio of the structure was 0.02. The equipment was simulated with a single DOF model and placed on the floor of the 13th story, and its mass was 10% of the story’s mass. Its damping ratio was 0.02 and the stiffness was 3.88 × 105 N/m. The material models for the equipment and structure were both bilinear kinematic hardening models, and their post-yield stiffness were 1% of the initial stiffness values.

The computational model of the equipment–HRS subsystem.
Cross section of beams and columns.
Based on the given parameters, a finite element model was constructed for the equipment–HRS subsystem by MATLAB programming. The beams and columns in the HRS were modeled with Euler–Bernoulli beam elements. Each beam element has two nodes, and each node has six DOFs, which are denoted ux, uy, uz, urotx, uroty, and urotz, respectively. For the structure, the mass matrix was created using the lumped mass method, and the damping matrix was based on Rayleigh damping. The single DOF equipment model was constructed of a spring-damper element. The equipment and structure were connected into an equipment–HRS subsystem model based on the displacement compatibility condition.
LNS model
Given that the soil is not totally nonlinear during strong earthquakes, an LNS model was created using the linear–nonlinear hybrid constraint modal substructure approach. The material parameters of the soil model are provided in Table 2. The LNS model has overall dimensions of 200 m × 200 m × 60 m, and a 2.50-m-thick foundation with a 31.5 m × 31.5 m cross section embedded in it.
Material parameters of the soil and foundation model.
The modeling of the soil and mesh discretization of the LNS model was done using ANSYS. The finite element model of the soil is made up of SOLID45, an 8-node hexahedral solid element, and each node has three translational DOFs, that is, translation along the x-axis, y-axis, and z-axis, respectively. Later, the constraint equations of a rigid foundation were defined with O representing the center of the foundation and S is an arbitrary point on the foundation excepting O. Since point O is a solid element node, it has no rotational DOFs. So, the Mass21 element was introduced to this point, and the real constants were set to the smaller values. Then point O has both translational and rotational DOFs, including translation along and rotation about the x-axis, y-axis, and z-axis. The DOFs of point O and S satisfy the following relations:
The range of nonlinear soil will be changed under earthquakes. In order to solve the problem, Wang and Jiang 19 implemented detailed study of the range of nonlinear soil and found that the nonlinear soil scope kept within three times the size of the foundation. This conclusion was also confirmed by trial calculation. So the nonlinear soil substructure β was designed to be three times the size of the foundation which retained 11,925 DOFs, and linear soil substructure α occupied the rest of the soil model, as illustrated in Figure 4. The DOFs of the linear soil substructure α were reduced through the constraint modal substructure approach, and its 2000-order dominant mode was retained based on the principle of potential energy truncation. 20 Besides, the nonlinear soil material model satisfied the Druck–Prager yield criterion. So an LNS model was created using the linear–nonlinear hybrid constraint modal substructure approach using ANSYS.

Substructure division of the soil model.
Then the equipment–HRS subsystem and LNS models were used for interactive numerical computation simulating the RTSSTT of the equipment–HRS–soil system. The three ground motions mentioned above were separately inputted as required by minor, moderate, and major earthquakes with level VIII intensity to obtain seismic responses of the system.
Feasibility assessment
The results of interactive numerical computation were compared with those of the integrated finite element analysis for the purpose of assessing the feasibility of the testing method proposed in this study. In the integrated finite element analysis, an integrated model of the system was created using ANSYS. The assessment focused on the acceleration responses at the center of soil during the minor earthquake, moderate earthquake, and major earthquake stages under TianJin ground motion of level VIII intensity. Figure 5 compares the acceleration responses at the center of soil output by the two methods. During the minor earthquake, moderate earthquake, and major earthquake stages, the peak acceleration at the center of soil in the interactive numerical computation increased 4.3%, 5.5%, and 7.8%, respectively, compared with those in the integrated finite element analysis. Besides, the time histories of acceleration obtained by the two methods are generally consistent during a minor earthquake and moderate earthquake, while they differ slightly during a major earthquake.

Comparison of acceleration responses at the center of soil in different computational methods.
The errors during a minor earthquake and moderate earthquake were small primarily because the soil response was linear or it experienced plastic deformation in a very limited zone in these two stages. The LNS model constructed using the linear–nonlinear hybrid constraint modal substructure approach proves to be able to capture the inherent dynamic behavior of the soil. During the major earthquake, the soil underwent relatively significant plastic deformation, and the LNS model could affect the actual seismic response of nonlinear soil behavior. Overall, the computational method incorporating the LNS model can offer the desired accuracy.
Computational efficiency is another important criterion used to assess the computational model. In this study, the number of DOFs of each model and the computation time were extracted for the two methods considered. Table 3 presents the number of DOFs of the models in the two methods. The system model created in the interactive numerical computation has 15,876 DOFs in total, a 72.5% decrease from the integrated finite element model’s DOFs 57,821. This suggests that the method proposed in this study can substantially reduce the amount of computation. Table 4 summarizes the computation time spent by each model during different earthquake stages under TianJin ground motion. MATLAB and ANSYS used the method “tic/toc” and “*get, B, active, 0, time, cpu,” respectively, to extract the computation time.
Number of DOFs of the models in the two computational methods.
DOFs: degrees of freedom.
Computation time for different earthquake stages in the two computational methods.
As can be seen in Table 4, the computation time spent by each software product significantly increased from minor earthquake to major earthquake in both the interactive numerical computation and integrated finite element analysis primarily because iterative computation needed to be performed after the equipment–HRS subsystem and the soil reached the nonlinear state. Besides, the total computation time spent in the interactive numerical computation method is significantly shorter than that in the other method. The computation time for the minor earthquake, moderate earthquake, and major earthquake stages were only 34.8%, 32.6%, and 29.5% of the corresponding computation time in the integrated finite element analysis. As the earthquake intensity increased, the decrease in computation time gradually accelerated. This is mainly because the LNS model used in interactive numerical computation has only a smaller number of vibration modes involved in the computation and thus avoids redundancy of DOFs. This method can significantly enhance the computational efficiency and lower the cost of computation by reducing the number of equations required for the computation and the number of iterations required to solve the nonlinear equations. Therefore, it can be employed to solve dynamic models of nonlinear systems efficiently and accurately.
Overall, the use of the LNS model can significantly improve the computational efficiency while ensuring accuracy. This confirms that the application of LNS model can provide powerful technical support for nonlinear soil participating in analysis. Meanwhile, as individual models are constructed for all substructures during interactive numerical computation, the difficulty of establishing integrated analysis model is considerably reduced. Moreover, the linear–nonlinear hybrid constraint modal substructure approach enables efficient and accurate modeling the nonlinearity of soil. The interactive numerical computation method proposed also offers an effective channel to compute and analyze complex equipment–HRS–soil systems.
The advantage of the proposed method lies in that the use of the LNS model can significantly improve the computational efficiency while ensuring accuracy. And interactive numerical computation can reduce the difficulty of modeling. However, the shortcoming of the proposed method lies in that it may be hard for the engineers for its demands of program skills and expertise.
Computation and analysis
To analyze the impact of equipment frequency changes in different systems, an equipment–HRS subsystem on a rigid ground and an equipment–HRS–soil system on Class III ground were constructed using the interactive numerical computation method. The frequency range of equipment was from 0 to 3 Hz. For the equipment–HRS subsystem on a rigid ground, ground motions were directly inputted for computation. The seismic responses of the equipment and structure were assessed based on the displacement of the equipment top relative to the floor and displacement of the structure top relative to the ground surface.
Figure 6 illustrates the peak displacements of the equipment top in the cases of the rigid ground and Class III ground. As shown in the figure, when the equipment’s frequency was fe≈0.43 Hz, (i.e. close to the fundamental frequency of the HRS), the peak displacement of the equipment top reached a noticeable maximum under each ground motion. After the SSI effect was taken into consideration (Class III ground), the maximum peak displacement of the equipment top shifted toward lower frequencies. At this maximum point, fe≈0.36 Hz, which means that the equipment’s frequency was approximately equal to the fundamental frequency of the HRS–soil subsystem in the case of Class III ground. The maximum peak displacements of the equipment top in different earthquake stages were considered and compared between the two cases. The results show that the maximum values during the minor earthquake, moderate earthquake, and major earthquake stages under El Centro ground motion in the case of Class III ground decreased by 25.7%, 18.3%, and 11.3%, respectively, compared with the corresponding values in the case of the rigid ground. The corresponding decreases under TianJin ground motion were 28.5%, 18.0%, and 9.0%, respectively. During the three earthquake stages under PerSon ground motion, the equipment top in the case of Class III ground experienced maximum peak displacements increased on 37.9%, 23.8%, and 13.4%, respectively, compared with those in the case of a rigid ground.

Peak displacements of the equipment top in the cases of rigid ground and Class III ground.
The peak displacements of the structure top in the cases of the rigid ground and Class III ground are plotted in Figure 7. The peak displacement of the structure top reached an apparent minimum under each ground motion. These minimum points correspond to the maximum points shown in Figure 6. In other words, the peak displacement of the structure top reached its minimum when the equipment’s frequency was close to the HRS’s fundamental frequency in the case of rigid ground, or when it was approximately equal to the HRS–soil subsystem fundamental frequency after the SSI effect was considered. The minimum peak displacements of the structure top in the three earthquake stages were considered and compared between the two cases. Under El Centro ground motion, the minimum values during a minor earthquake, moderate earthquake, and a major earthquake in the case of Class III ground decreased by 12.7%, 9.6%, and 6.7%, respectively, compared to those in the case of the rigid ground. Under TianJin ground motion, the structure top showed minimum peak displacement decreases of 18.1%, 15.7%, and 10.8%, respectively, compared to those in the case of the rigid ground. Under PerSon ground motion, the corresponding decreases in the minimum peak displacement of the structure top were 8.7%, 6.3%, and 5.8%.

Peak displacements of the structure top in the cases of rigid ground and Class III ground.
An analysis of the seismic responses reveals that after the effect of equipment–HRS–soil interaction was taken into account, the peak displacement of the equipment top exhibited a noticeable maximum while the peak displacement of the structure top had an obvious minimum. The equipment’s frequency at this point was close to the HRS’s fundamental frequency in the case of rigid ground or approximately equal to the HRS–soil subsystem’s frequency after the SSI effect was considered. The analysis above concludes that, after the effect of SSI on the equipment was considered, the peak displacements of the equipment top in the case of Class III ground tended to decrease from those in the case of a rigid ground under El Centro and TianJin ground motions but that tended to increase under PerSon ground motion. By contrast, the peak displacement of the structure top tended to decrease significantly when considered the SSI effect, indicating that the SSI effect enhances the structure’s safety in the event of an earthquake. As the earthquake intensity increased, the magnitude of a change in the structure’s minimum response tended to decrease, suggesting that the positive effect of SSI on the structure weakens with increasing earthquake intensity.
Conclusion
The interactive numerical computation was proposed to achieve a simulated process of RTSSTT which introduced an LNS model considering that the soil is not completely in nonlinear status under strong earthquakes. This study leads to the following conclusions:
The equation of motion for the equipment–HRS–LNS system was derived through a combination of the branch modal substructure and linear–nonlinear hybrid constraint modal substructure approaches, which is characterized by clear definitions, ease of use, and so on.
The interactive numerical computation process was successfully simulated by interactive numerical computation using MATLAB and ANSYS. The continuous reading and writing of data by the software provided computational data communication.
Compared to integrated finite element analysis, the interactive numerical computation method can significantly improve the computational efficiency while ensuring accuracy. This confirms that the application of LNS model can provide powerful technical support for nonlinear soil participating in analysis. And the interactive numerical computation method also offers an effective channel to compute and analyze complex equipment–HRS–soil systems.
When the equipment’s frequency was close to the fundamental frequency of the HRS in the case of rigid ground, or when it was close to the fundamental frequency of the HRS–soil subsystem after the SSI effect was considered, the peak displacement of the equipment top reached a noticeable maximum while the peak displacement of the structure top showed an obvious minimum. Given that the seismic responses of equipment and HRS will change substantially when the SSI effect is taken into account, it is necessary to investigate the seismic response of the whole equipment–HRS–soil system in the seismic design of high buildings.
Footnotes
Handling Editor: Elsa de Sa Caetano
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was financially supported by the National Natural Science Foundation of China (Research Project Nos 51478312, 51208356, and 51278335).
