Substructure Hybrid Simulation Boundary Technique Based on Beam/Column Inflection Points

Compatibility among substructures is an issue for hybrid simulation. Traditionally, the structure model is regarded as the idealized shear model. The equilibrium and compatibility of the axial and rotational direction at the substructure boundary are neglected. To improve the traditional boundary technique, this paper presents a novel substructure hybrid simulation boundary technique based on beam/column inflection points, which can effectively avoid the complex operation for realizing the bending moment at the boundary by using the features of the inflection point where the bending moment need not be simulated in the physical substructure. An axial displacement prediction technique and the equivalent force control method are used to realize the proposed method. The numerical simulation test scheme for the different boundary techniques was designed to consider three factors: (i) the different structural layers; (ii) the line stiffness ratio of the beam to column; and (iii) the peak acceleration. The simulation results for a variety of numerical tests show that the proposed technique shows better performance than the traditional technique, demonstrating its potential in improving HS test accuracy. Finally, the accuracy and feasibility of the proposed boundary technique is verified experimentally through the substructure hybrid simulation tests of a six-story steel frame model.


Introduction
Hybrid simulation (HS), also known as pseudo-dynamic testing, was initially proposed in 1969 [1]. This approach [2][3][4], which is a combination of physical testing and numerical simulation, has been widely used to evaluate the seismic performance of a structure, because it is efficient and cost-effective. HS takes the critical parts of a structure, which cannot be simulated adequately with the numerical simulation, as the physical substructure to test in the lab. The rest of the structure is taken as the numerical substructure and simulated using appropriate finite element software. Simulation of the laboratory and computer parts of the structure are carried out simultaneously, exchanging information about the responses at the boundary at each time step. Furthermore, the HS method transforms a real-time dynamic test under earthquake loads into a static loading test by simulating the dynamic components of the structural response in the computer, thus reducing the requirement for the loading devices and making large-scale or even full-scale model tests possible.
Initial research on the HS was overwhelmingly focused in Japan and the United States. Beginning in the 1990s, development and application of HS in Europe, Asia, and other regions grew significantly [5][6][7]. For structural systems with rate-dependent components, real-time HS testing was proposed and developed [8][9][10][11][12][13][14][15]. Moreover, HS also allows testing to be geographically distributed in multiple testing facilities. Many countries have been begun to build networked structures laboratories and developed a variety of HS platforms [16][17][18][19] to allow performance assessment of modern building structure systems. The largest is the Network for Earthquake Engineering Simulation (NEES) supported by the United States National Science Foundation from 1999 to 2014 [20][21][22][23]. NEES featured 14 geographically distributed and shared-use laboratories to test modern structural system.
Although substantial advances have been made in HS, achieving equilibrium and compatibility of the boundaries between the two substructures remains problematic, possibly limiting the accuracy of HS tests. There are currently four Multi-Axis Substructure Testing MAST systems in the world including MAST at the University of Illinois at Urbana-Champaign (UIUC), the University of Minnesota (USA), Swinburne University of Technology (Australia) and Ecole Polytechnique Montreal (Canada) that have the capabilities of six-DOF switched and mixed displacement and force modes of control [24][25][26]. The MAST at Swinburne has already been used in such hybrid simulation tests. However, controlling and measuring the force and displacement in all DOF is challenging. The method we present is relatively simple and can be widely implemented.
The fact that it cannot consider the effect of bending moments has always been a flaw for the HS method. However, this method has also been widely used in various structures due to its effectiveness, not just limited to shear buildings, as described in the literature [27][28][29][30][31][32][33][34][35]. In these tests, the conventional method of placing the loading position at each layer was widely adopted without considering the effect of the bending moment.
This paper proposes and improves a novel boundary technique based on a beam's inflection point that can avoid the complex operation for realizing bending moments at the substructure boundary, based on previous a series of experimental techniques [34,36,37]. Following a description of the proposed boundary condition technique, a discussion of various implementation challenges is presented. The proposed method was validated through virtual HS, and experimental validation was carried out using a six-story steel frame model. The results demonstrate the efficacy of the proposed approach.

Boundary Technique for the HS
The compatibility and equilibrium conditions should be satisfied at the boundary between the physical and the numerical substructures during the HS. Consider the five-story example shown in Figure 1. Figure 1a shows the full model of the structure. Figure 1b illustrates the traditional boundary between two substructures, which is located at the top of the first story. This traditional boundary strategy mainly takes the structure model as an idealized shear model. The equilibrium of shear forces and the compatibility of the horizontal or shear displacement is satisfied by the hydraulic actuator. The equilibrium of the bending moment and axial force, as well as the compatibility axial and rotational deformations are typically neglected. Consequently, only two kinds of loading device are required at the boundary between the numerical substructure and the physical substructure: (i) one is the horizontal loading device to impose the shear force or displacement; and (ii) the other is the vertical loading device to simulate the gravity load due to upper stories, which is usually taken to be constant during the test. Therefore, several inevitable errors are introduced by not satisfying compatibility and equilibrium conditions.
To meet the equilibrium and compatibility in the axial and rotational direction at the boundary, more than two hydraulic actuators are required, as well as appropriate sensors to measure the rotational angle and the axial deflection, which can be challenging. To address the problems, a new boundary technique based on the inflection point is proposed, as shown in Figure 1c. The bending moments at the section of the inflection point in a structure is zero; therefore, this new boundary technique can effectively avoid the complex operation for realizing the bending moment at the boundary. However, two new problems arise in the proposed boundary technique: (i) how to calculate the horizontal displacement of the inflection point; and (ii) how to meet equilibrium and compatibility in the axial direction. technique; and (c) proposed technique. di denotes the displacements of each level on the full models. i ∈ [1,5]. db denotes the displacements at the inflection point. δ1 and δ2 denote the axial deformation on interfaces. NG denotes the constant representative value of gravity load calculated by the numerical substructure. Nb1 and Nb2 denote the axial forces resulted from the overturning moments on physical substructure transmitted from numerical substructure on each interface) For the first problem, assuming a uniform column between two floors, then the displacement with respect to the lower floor at the inflection point is a linear interpolation value of the interstory drift for the corresponding floor. Adding a new concentrated mass or a degree of freedom (DOF) on the inflection point of each story is the direct method to solve the first problem for the proposed boundary technique. However, the DOF of the equation of motion (EOM) will increase significantly so that the stability and efficiency of the HS is influenced. A linear interpolation technique is used to determine the displacement at the inflection point. The calculated displacement at the inflection point is from the interstory drift of the corresponding floor. Therefore, the solution of the EOM is not different from the traditional method.
For the second problem, the axial displacement prediction technique at the boundary is presented, because the axial displacement is usually too small to measure accurately using displacement sensors [25]. Moreover, the axial displacement at the boundary must be obtained to start the simulation of the numerical substructure for the next step, so that the axial forces can be obtained from the numerical substructure simulation results to apply to the boundary of the physical substructure by the hydraulic actuators. Traditionally, an iteration strategy should be introduced to calculate the axial displacement. However, such iteration strategies will result in inevitable measure errors, because the physical substructure will be loaded many times during the iteration. To avoid iteration, an axial displacement prediction technique is presented to predict the axial displacement using the loading axial force, the assumption of the elastic axial deformation during the test, and the equilibrium equation.
The proposed prediction method is given by the following equation i ∈ [1,5]. d b denotes the displacements at the inflection point. δ 1 and δ 2 denote the axial deformation on interfaces. N G denotes the constant representative value of gravity load calculated by the numerical substructure. N b1 and N b2 denote the axial forces resulted from the overturning moments on physical substructure transmitted from numerical substructure on each interface.
For the first problem, assuming a uniform column between two floors, then the displacement with respect to the lower floor at the inflection point is a linear interpolation value of the interstory drift for the corresponding floor. Adding a new concentrated mass or a degree of freedom (DOF) on the inflection point of each story is the direct method to solve the first problem for the proposed boundary technique. However, the DOF of the equation of motion (EOM) will increase significantly so that the stability and efficiency of the HS is influenced. A linear interpolation technique is used to determine the displacement at the inflection point. The calculated displacement at the inflection point is from the interstory drift of the corresponding floor. Therefore, the solution of the EOM is not different from the traditional method.
For the second problem, the axial displacement prediction technique at the boundary is presented, because the axial displacement is usually too small to measure accurately using displacement sensors [25]. Moreover, the axial displacement at the boundary must be obtained to start the simulation of the numerical substructure for the next step, so that the axial forces can be obtained from the numerical substructure simulation results to apply to the boundary of the physical substructure by the hydraulic actuators. Traditionally, an iteration strategy should be introduced to calculate the axial displacement. However, such iteration strategies will result in inevitable measure errors, because the physical substructure will be loaded many times during the iteration. To avoid iteration, an axial displacement prediction technique is presented to predict the axial displacement using the loading axial force, the assumption of the elastic axial deformation during the test, and the equilibrium equation.
The proposed prediction method is given by the following equation in which δ N (i) is the prediction axial deformation at the step i; m is the number of stories in the physical substructure; N j (i) is the calculated axial force of the j-story at the step i; N b (i) is the calculated axial force of the boundary story, which is calculated from the numerical substructure; h j is the height of the j-story of the physical substructure; h I,b is the height from the inflection point to the floor of the boundary story; E j A j is the axial stiffness of the jth-story of the physical substructure, which is constant during the test due to the assumption of the elastic axial deformation; and E b A b is the axial stiffness of the boundary story of the physical substructure. Here, N j (i) can be calculated by the equilibrium equation. For example, consider the top story of the physical substructure shown in Figure 2. Here, L is the span of the frame; point O m,1 and O m,2 are inflection points of the m story; h I,m is the height from the inflection point to the roof of the m story; R b is the applied force by the actuator at the boundary; and R m is the applied force by the actuator at the mth-story. Noting that the dynamics of the frame are represented in the computer and using the moment equilibrium equation, the axial force can be calculated by taking moments about point O 1 or O 2 . In the same way, the axial force of each story of the physical substructure can be calculated from the results of the upper story. Then, the axial forces of the physical substructure can be obtained. Subsequently, the prediction axial displacement at the boundary can be calculated using Equation (1), i.e., () j Ni is the calculated axial force of the j-story at the step i; b () Ni is the calculated axial force of the boundary story, which is calculated from the numerical substructure; j h is the height of the j-story of the physical substructure; Ib , h is the height from the inflection point to the floor of the boundary story; j j A E is the axial stiffness of the jth-story of the physical substructure, which is constant during the test due to the assumption of the elastic axial deformation; and bb EA is the axial stiffness of the boundary story of the physical substructure.

HS Component Interaction
Typically, the HS test includes three main components: the physical substructure, the numerical substructure, and the numerical integration routine. Early HS tests usually used an idealized mass-spring-damper model whose performance is known a priori to simulate the numerical substructure, allowing the software controlling the entire HS test to be realized as in a single program (e.g., written in FORTRAN or C). As HS testing has developed, researchers have sought to use more and more complicated models, in terms of both experimental and computational components. Computational components tailored for a particular strength are combined with general purpose finite element programs to achieve more accurate simulation. For example, a

HS Component Interaction
Typically, the HS test includes three main components: the physical substructure, the numerical substructure, and the numerical integration routine. Early HS tests usually used an idealized mass-spring-damper model whose performance is known a priori to simulate the numerical substructure, allowing the software controlling the entire HS test to be realized as in a single program (e.g., written in FORTRAN or C). As HS testing has developed, researchers have sought to use more and more complicated models, in terms of both experimental and computational components. Computational components tailored for a particular strength are combined with general purpose finite element programs to achieve more accurate simulation. For example, a high-fidelity ABAQUS model was used [38], combined with the Matlab-based UI-SimCor [17], to investigate the seismic performance of semirigid steel joints. To address current research needs, HS testing software must accommodate different computational components. Communication between the computational and experimental components must be considered to realize a HS test. Figure 3 shows the sketch of the HS component interaction. First, OpenSEES is used to simulate both the physical and the numerical substructure in the virtual HS test, while MATLAB is employed to integrate the EOM. To realize the required data transfer, a communication bridge between the two different computational components is constructed. Data communication will take place at each step; each independent unit writes data to the other units and waits for the data written by other units. Moreover, the waiting or calculating time at different step is not uniform. This paper employs the network socket technique as a bridge to communicate the data between the two different computational components. In the MATLAB component, the TCPIP and ECHOTCPIP commands are used to construct an ICPIP object associated with OpenSEES. Concurrently, the command SOCKET in OpenSEES program is used to communicate with the Matlab program. Subsequently, a data acquisition card of PCI 1712 manufactured by the Advantech is employed as the communication bridge between the MATLAB and the MTS loading device. The PCI-1712 card provides a total of 16 single-ended or eight differential A/D input channels or a mixed combination, two 12-bit D/A output channels, 16 digital input/output channels, and three 10 MHz 16-bit multifunction counter channels. high-fidelity ABAQUS model was used [38], combined with the Matlab-based UI-SimCor [17], to investigate the seismic performance of semirigid steel joints. To address current research needs, HS testing software must accommodate different computational components. Communication between the computational and experimental components must be considered to realize a HS test. Figure 3 shows the sketch of the HS component interaction. First, OpenSEES is used to simulate both the physical and the numerical substructure in the virtual HS test, while MATLAB is employed to integrate the EOM. To realize the required data transfer, a communication bridge between the two different computational components is constructed. Data communication will take place at each step; each independent unit writes data to the other units and waits for the data written by other units. Moreover, the waiting or calculating time at different step is not uniform. This paper employs the network socket technique as a bridge to communicate the data between the two different computational components. In the MATLAB component, the TCPIP and ECHOTCPIP commands are used to construct an ICPIP object associated with OpenSEES. Concurrently, the command SOCKET in OpenSEES program is used to communicate with the Matlab program. Subsequently, a data acquisition card of PCI 1712 manufactured by the Advantech is employed as the communication bridge between the MATLAB and the MTS loading device. The PCI-1712 card provides a total of 16 single-ended or eight differential A/D input channels or a mixed combination, two 12-bit D/A output channels, 16 digital input/output channels, and three 10 MHz 16-bit multifunction counter channels.

Calculated program
Physical Substructure Numerical substructure Entire model Figure 3. Sketch of the HS component interaction.

Equivalent Force Control Method
Research on various integration algorithms employed for HS testing is very active. For example, the CR and KR-α methods [10,35] which have unconditional stability for linear elastic and stiffness softening nonlinear systems, were presented. However, the implicit algorithms generally require iteration, which can result in inconsistent HS test results. To avoid iteration, an equivalent force control (EFC) method [39], which uses feedback control to replace numerical iteration for the real-time HS testing, was proposed. EFC method was used to a full-scale six-story frame-supported reinforced concrete masonry shear wall [39].
This proposed approach combines the EFC method with the proposed inflection point boundary condition technique to conduct HS tests. The EFC method is based on the feedback control

Equivalent Force Control Method
Research on various integration algorithms employed for HS testing is very active. For example, the CR and KR-α methods [10,35] which have unconditional stability for linear elastic and stiffness softening nonlinear systems, were presented. However, the implicit algorithms generally require iteration, which can result in inconsistent HS test results. To avoid iteration, an equivalent force control (EFC) method [39], which uses feedback control to replace numerical iteration for the real-time HS testing, was proposed. EFC method was used to a full-scale six-story frame-supported reinforced concrete masonry shear wall [39].
This proposed approach combines the EFC method with the proposed inflection point boundary condition technique to conduct HS tests. The EFC method is based on the feedback control method shown in Figure 4. Here, the F EQ,i+1 , which is used for the equivalent force (EF) command, can be calculated by the information before the i + 1 step. Therefore, the F EQ,i+1 is available and remains unchanged at i + 1 step. The F EQ,i+1 , which is used for equivalent force (EF) feedback, is obtained by the three terms: (i) the reaction of the experimental substructure R E d c E,i+1 (t) , which is measured by the load cells; (ii) the reaction of the numerical substructure R N d c E,i+1 (t) , which is obtained from the results of the numerical simulation; and (iii) the pseudo-dynamic force, which is linearly related to displacement by K PD . The displacement command, d c i+1 (t), is determined from the force signal using the force-displacement conversion matrix C F . The equivalent force controller is used to make certain that the EF feedback can track the EF command quickly without overshoot. This paper uses a PD controller as the equivalent force controller. The parameters of the PD controller are obtained as discussed in Reference [39].
Sustainability 2018, 10, x FOR PEER REVIEW 6 of 20 method shown in Figure 4. Here, the FEQ,i+1, which is used for the equivalent force (EF) command, can be calculated by the information before the i + 1 step. Therefore, the FEQ,i+1 is available and remains unchanged at i + 1 step. The FEQ,i+1, which is used for equivalent force (EF) feedback, is obtained by the three terms: (i) the reaction of the experimental substructure , which is measured by the load cells; (ii) the reaction of the numerical substructure , which is obtained from the results of the numerical simulation; and (iii) the pseudo-dynamic force, which is linearly related to displacement by KPD. The displacement command, , is determined from the force signal using the force-displacement conversion matrix CF. The equivalent force controller is used to make certain that the EF feedback can track the EF command quickly without overshoot. This paper uses a PD controller as the equivalent force controller. The parameters of the PD controller are obtained as discussed in Reference [39].  Figure 5 shows the flow diagram for the proposed HS for both the virtual and the physical HS. The only difference between the virtual and the physical HS is that the physical substructure is modeled computationally. For the physical HS, the results of the physical substructure can be measured by the sensors from the model in the lab. For the virtual HS, the results of the physical substructure can be obtained from the results of an additional numerical simulation using the OpenSEES. In Figure 5, istep is the ith step of the earthquake records, and n corresponds the total number of the earthquake records. Matlab is used to integrate the EOM for calculating the responses of the structure. OpenSEES is used for simulation of the numerical substructure. The basic detailed process is as follows:


Step 1: The initial stiffness matrix, the initial mass matrix, the initial damping matrix, other initial parameters, and the ground motion record are obtained. The ground motion record may be scaled in amplitude to accommodate the goal of different HS tests.  Step 2: The EF command is calculated. The displacement command is then calculated using the force-displacement conversion matrix and the EF controller.  Step 3: The pseudo-dynamic force is calculated using the pseudo-dynamic stiffness KPD, which can be determined from the initial parameters. Concurrently, the predicted axial deformation is obtained using Equation (1). Before the next step, all calculations are realized using Matlab. After this step, the calculated displacement command of the numerical substructure and the prediction axial deformation at the boundary is applied to the numerical substructure modeled by OpenSEES.  Step 4: The calculated displacement command of the numerical substructure and the prediction axial deformation at the boundary are sent to the numerical substructure. Subsequently, the  Figure 5 shows the flow diagram for the proposed HS for both the virtual and the physical HS. The only difference between the virtual and the physical HS is that the physical substructure is modeled computationally. For the physical HS, the results of the physical substructure can be measured by the sensors from the model in the lab. For the virtual HS, the results of the physical substructure can be obtained from the results of an additional numerical simulation using the OpenSEES. In Figure 5, istep is the ith step of the earthquake records, and n corresponds the total number of the earthquake records. Matlab is used to integrate the EOM for calculating the responses of the structure. OpenSEES is used for simulation of the numerical substructure. The basic detailed process is as follows:

•
Step 1: The initial stiffness matrix, the initial mass matrix, the initial damping matrix, other initial parameters, and the ground motion record are obtained. The ground motion record may be scaled in amplitude to accommodate the goal of different HS tests.

•
Step 2: The EF command is calculated. The displacement command is then calculated using the force-displacement conversion matrix and the EF controller.

•
Step 3: The pseudo-dynamic force is calculated using the pseudo-dynamic stiffness K PD , which can be determined from the initial parameters. Concurrently, the predicted axial deformation is obtained using Equation (1). Before the next step, all calculations are realized using Matlab. After this step, the calculated displacement command of the numerical substructure and the prediction axial deformation at the boundary is applied to the numerical substructure modeled by OpenSEES.

•
Step 4: The calculated displacement command of the numerical substructure and the prediction axial deformation at the boundary are sent to the numerical substructure. Subsequently, the reaction force at the boundary is calculated. Similarly, the data calculated by OpenSEES are sent to the Matlab program and subsequently wait to receive the data calculated by Matlab for the next step.

•
Step 5: The calculated reaction force in the vertical direction and the displacement command in the horizontal direction are sent to the actuator controllers for the physical substructure. Then, the reaction of the physical substructure is obtained by: (i) the sensors in the associated loading devices; or (ii) OpenSEES for the virtual test. • Step 6: The EF feedback is obtained by summing up the pseudo-dynamic force, and the reaction of both the numerical substructure and the physical substructure. Then, compare the EF feedback with the EF command to obtain the error between them. If the error between the EF feedback with the EF command is less than the specified tolerance, go the next step. If not, repeat the steps from the third step to the sixth step.

•
Step 7: The responses of the entire structure are obtained. Then, repeat the above-mentioned steps until the entire ground motion record has been processed.
Sustainability 2018, 10, x FOR PEER REVIEW 7 of 20 reaction force at the boundary is calculated. Similarly, the data calculated by OpenSEES are sent to the Matlab program and subsequently wait to receive the data calculated by Matlab for the next step.  Step 5: The calculated reaction force in the vertical direction and the displacement command in the horizontal direction are sent to the actuator controllers for the physical substructure. Then, the reaction of the physical substructure is obtained by: (i) the sensors in the associated loading devices; or (ii) OpenSEES for the virtual test.  Step 6: The EF feedback is obtained by summing up the pseudo-dynamic force, and the reaction of both the numerical substructure and the physical substructure. Then, compare the EF feedback with the EF command to obtain the error between them. If the error between the EF feedback with the EF command is less than the specified tolerance, go the next step. If not, repeat the steps from the third step to the sixth step.  Step 7: The responses of the entire structure are obtained. Then, repeat the above-mentioned steps until the entire ground motion record has been processed. After the details of proposed HS, a series of numerical simulation is carried out though the virtual HS in the next section. Then, experimental validation is accomplished using a six-story steel frame model. The El Centro (N-S, 1940) is taken as the ground motion record for both the virtual and the physical HS tests.

Numerical Simulation Model Based on the Uniform Design
To illustrate the proposed approach, a series of steel frame structures with 3-17 stories were considered. A nonlinear Beam-column element, which has a hardening bilinear constitutive model, was chosen in OpenSEES to model the frame. The yield strength is 235 Mpa. The elastic modulus is 200,000 Mpa. The post-yield to pre-yield stiffness ratio is 0.02. The mass matrix was calculated using provision 5.1.3 of the Chinese Code for Seismic Building Design GB50011-2010. Rayleigh Damping was adopted with a damping ratio of 0.05 in the first two modes. The number of the stories in the After the details of proposed HS, a series of numerical simulation is carried out though the virtual HS in the next section. Then, experimental validation is accomplished using a six-story steel frame model. The El Centro (N-S, 1940) is taken as the ground motion record for both the virtual and the physical HS tests.

Numerical Simulation Model Based on the Uniform Design
To illustrate the proposed approach, a series of steel frame structures with 3-17 stories were considered. A nonlinear Beam-column element, which has a hardening bilinear constitutive model, was chosen in OpenSEES to model the frame. The yield strength is 235 Mpa. The elastic modulus is 200,000 Mpa. The post-yield to pre-yield stiffness ratio is 0.02. The mass matrix was calculated using provision 5.1.3 of the Chinese Code for Seismic Building Design GB50011-2010. Rayleigh Damping was adopted with a damping ratio of 0.05 in the first two modes. The number of the stories in the building, the beam-column's linear stiffness ratio (LSR), and the peak ground acceleration (PGA) of the ground motion record were taken as parameters for the numerical study. Each story height is 4 m, and the span of the steel frame structure is 6 m.
This experimental design method proposed by Wang and Fang [40] was used to increase efficiency and reduce the required number of tests by focusing on the uniformity of the test matrix. For example, assume that the m parameters are chosen in an experiment. Let the number of values of each parameter to be considered be n. The number of experiments is m n for a full factorial design. The number of samples is about m 2 for an orthogonal design, which takes both the orthogonality and uniformity of the design space into consideration. However, the number of experiments is about m for the uniform design, which only highlights the uniformity of the design space.
The uniform design table, Un(q s ), which is similar to the orthogonal design table, was constructed in advance, where the symbol U is the label for the uniform design, n is the number of test, s is the maximum number of parameters, and q is the number of levels for each parameter. The deviation D of the uniform design table is taken as the index to evaluate the uniformity of the table. The smaller is the value of D, the better is the uniformity of the design table. By considering the number of stories, the beam-column's linear stiffness ratio, and the peak ground acceleration (PGA) of the ground motion record, the uniform design table U* 15 (15 7 ) was chosen. The symbol * means that the design table has more uniformity [40]. The deviation D of the selected uniform table is 0.1361, which meets the precision required. More detailed parameters are listed in Table 1.

Numerical Simulation Results Analysis
This section provides a comparison of the numerical simulation for each model listed in Table 1 for the traditional technique and the proposed elastic inflection point technique. The simulation result of entire structure using OpenSEES is considered to be the "exact" solution and used as a point of reference. The errors in the numerical simulation for the different boundary simulation techniques at the ith step, as compared to the reference solution, are designated as e i . Then, the cumulative error in the numerical simulation for the different boundary simulation techniques is given by Due to space limitation, only representative results are given for three of these buildings: tests No. 2, No. 8, and No. 13. Figure 6 shows the time history of the errors e i for tests No. 2, No. 8, and No. 13 for the traditional and proposed boundary simulation techniques. In Figure 6, we can see clearly that the error in the displacement for the proposed boundary technique is smaller than the traditional technique. The simulation results for the other numerical tests shown in Table 1 followed similar trends. Additionally, the maximum displacement errors and the cumulative errors for the different boundary techniques were calculated for the 15 tests shown in Table 1. The mean of the maximum displacement error and the mean of the cumulative error over all 15 tests was calculated for both the traditional technique and the proposed technique. Both error measures were found to be an order of magnitude smaller for the proposed approach, as compared with the traditional technique. The maximum displacement error ratio between two different techniques had an average of 9.18 with a maximum of 13.00 and a minimum of 2.66. The cumulative error ratio between two different techniques had an average of 9.87 with a maximum of 18.04 and a minimum of 3.12.  Figure 6 shows the time history of the errors ei for tests No. 2, No. 8, and No. 13 for the traditional and proposed boundary simulation techniques. In Figure 6, we can see clearly that the error in the displacement for the proposed boundary technique is smaller than the traditional technique. The simulation results for the other numerical tests shown in Table 1 followed similar trends. Additionally, the maximum displacement errors and the cumulative errors for the different boundary techniques were calculated for the 15 tests shown in Table 1. The mean of the maximum displacement error and the mean of the cumulative error over all 15 tests was calculated for both the traditional technique and the proposed technique. Both error measures were found to be an order of magnitude smaller for the proposed approach, as compared with the traditional technique. The maximum displacement error ratio between two different techniques had an average of 9.18 with a maximum of 13.00 and a minimum of 2.66. The cumulative error ratio between two different techniques had an average of 9.87 with a maximum of 18.04 and a minimum of 3.12.    The maximum internal forces (the bending moment, the shear force and the axial force) errors and the cumulative errors of different boundary techniques at the top of the first story column follow a similar tendency with the displacement. The means of maximum errors ratio and the cumulative errors ratio of the bending moment between the traditional and proposed technique are 4.86 and 9.42, respectively. Both the shear force and the axial force have a similar tendency but the means of the corresponding ratio are larger than the bending moment. The means of maximum errors ratio and the cumulative errors ratio of the shear force are 9.44 and 12.90, respectively. The means of maximum errors ratio and the cumulative errors ratio of the axial force are 17.62 and 23.65, respectively. Therefore, these results demonstrate that the proposed technique based on the inflection point can improve substantially the test accuracy and show the better performance than the traditional method.    4.86 and 9.42, respectively. Both the shear force and the axial force have a similar tendency but the means of the corresponding ratio are larger than the bending moment. The means of maximum errors ratio and the cumulative errors ratio of the shear force are 9.44 and 12.90, respectively. The means of maximum errors ratio and the cumulative errors ratio of the axial force are 17.62 and 23.65, respectively. Therefore, these results demonstrate that the proposed technique based on the inflection point can improve substantially the test accuracy and show the better performance than the traditional method.

Experimental Validation
To validate the proposed boundary method, a six-story steel frame structure, in which the first story was taken as the physical substructure and the rest of the structure as the numerical substructure, was designed as the HS testing model. Ground motions with PGAs of 310, 620, 800, 1200, and 1600 cm/s 2 were sequentially applied to the test structure.

Test Model
The six-story frame structure was designed based on both Chinese Code GB50017-2003 for Design of Steel Structures and Chinese Code for Seismic Building Design GB50011-2010. The height of the first story is 4200 mm and others are 4000 mm. The span of the bay for the test model is 8000 mm. The sections of the columns are HW600 × 600 × 20 × 30, and the section of the beams are HW500 × 300 × 14 × 22 (with units of mm). Q235 steel is used. Because of laboratory limitations, a 1/2-scale test model was manufactured in the Seismic Experimental Center of Civil Engineering at the Harbin Institute of Technology.
The first story, along with the columns of the second story up to the inflection point, was taken as the physical substructure for the proposed boundary technique, as illustrated in Figure 10. A transfer-structure, which is hinge connected with the inflection point of both columns by the high-strength bolts M10.9, was designed to realize the synchronous horizontal displacement by

Experimental Validation
To validate the proposed boundary method, a six-story steel frame structure, in which the first story was taken as the physical substructure and the rest of the structure as the numerical substructure, was designed as the HS testing model. Ground motions with PGAs of 310, 620, 800, 1200, and 1600 cm/s 2 were sequentially applied to the test structure.

Test Model
The six-story frame structure was designed based on both Chinese Code GB50017-2003 for Design of Steel Structures and Chinese Code for Seismic Building Design GB50011-2010. The height of the first story is 4200 mm and others are 4000 mm. The span of the bay for the test model is 8000 mm. The sections of the columns are HW600 × 600 × 20 × 30, and the section of the beams are HW500 × 300 × 14 × 22 (with units of mm). Q235 steel is used. Because of laboratory limitations, a 1/2-scale test model was manufactured in the Seismic Experimental Center of Civil Engineering at the Harbin Institute of Technology.
The first story, along with the columns of the second story up to the inflection point, was taken as the physical substructure for the proposed boundary technique, as illustrated in Figure 10. A transfer-structure, which is hinge connected with the inflection point of both columns by the high-strength bolts M10.9, was designed to realize the synchronous horizontal displacement by using the two-equilateral angle steel 100 × 10 (mm). The physical experiment is shown in Figure 11, and Table 2 provides the test results for the performance of the steel material. The mass of the first-fifth stories m 1-5 = 88,000 kg, and the mass of the topper story m 6 = 84,000 kg. The record of the ground motion is the El Centro. Rayleigh Damping is chosen, with parameters a is 0.115 and b is 0.0024, corresponding to 5% damping in the first two modes.
using the two-equilateral angle steel 100 × 10 (mm). The physical experiment is shown in Figure 11, and Table 2 provides the test results for the performance of the steel material. The mass of the firstfifth stories m1-5 = 88,000 kg, and the mass of the topper story m6 = 84,000 kg. The record of the ground motion is the El Centro. Rayleigh Damping is chosen, with parameters a is 0.115 and b is 0.0024, corresponding to 5% damping in the first two modes.  using the two-equilateral angle steel 100 × 10 (mm). The physical experiment is shown in Figure 11, and Table 2 provides the test results for the performance of the steel material. The mass of the firstfifth stories m1-5 = 88,000 kg, and the mass of the topper story m6 = 84,000 kg. The record of the ground motion is the El Centro. Rayleigh Damping is chosen, with parameters a is 0.115 and b is 0.0024, corresponding to 5% damping in the first two modes.

Displacement-Force Mixed Control Technique
The traditional HS test, which usually focuses on the compatibility and equilibrium only in the horizontal direction, is conducted using the displacement control. The calculated displacements obtained by solving the EOM were applied to the physical substructure test model. Then, the reaction forces corresponding to the command displacements were measured and fed back to the EOM for the calculation of the next displacement increment. The stiffness of the columns in the vertical (axial) direction is typically high, consequently making displacement control challenging to implement, primarily because the resolution of the displacement transducer used for feedback is inadequate. Therefore, force control was used to realize the axial loading in the proposed approach. To address the problem of measuring the axial displacement of the column, the technique proposed in Section 3.2 for predicting the axial displacement was employed.
To realize the displacement-force mixed control, the arrangement of the vertical and horizontal loading devices was designed as shown in Figure 12. Two hydraulic jacks were chosen as the vertical loading devices while two hydraulic servo actuators as the horizontal loading devices.

Displacement-Force Mixed Control Technique
The traditional HS test, which usually focuses on the compatibility and equilibrium only in the horizontal direction, is conducted using the displacement control. The calculated displacements obtained by solving the EOM were applied to the physical substructure test model. Then, the reaction forces corresponding to the command displacements were measured and fed back to the EOM for the calculation of the next displacement increment. The stiffness of the columns in the vertical (axial) direction is typically high, consequently making displacement control challenging to implement, primarily because the resolution of the displacement transducer used for feedback is inadequate. Therefore, force control was used to realize the axial loading in the proposed approach.
To address the problem of measuring the axial displacement of the column, the technique proposed in Section 3.2 for predicting the axial displacement was employed.
To realize the displacement-force mixed control, the arrangement of the vertical and horizontal loading devices was designed as shown in Figure 12. Two hydraulic jacks were chosen as the vertical loading devices while two hydraulic servo actuators as the horizontal loading devices.  Figure 13a shows the vertical loading devices which are two hydraulic jacks. The hydraulic jacks were arranged at the boundary of the physical substructure model where the self-balancing device (Figure 13a) was designed to realize the axial load at the boundary. The self-balancing device includes the 20 mm-thick steel plate which is placed on the top of the hydraulic jack and four 25 mm-diameter steel bars which connects the 20 mm-thick steel plate and the mudsill of the physical substructure model. The loading capacity of the hydraulic jacks is 1000 kN. The loading errors are 1% of the loading capacity.  Figure 13a shows the vertical loading devices which are two hydraulic jacks. The hydraulic jacks were arranged at the boundary of the physical substructure model where the self-balancing device (Figure 13a) was designed to realize the axial load at the boundary. The self-balancing device includes the 20 mm-thick steel plate which is placed on the top of the hydraulic jack and four 25 mm-diameter steel bars which connects the 20 mm-thick steel plate and the mudsill of the physical substructure model. The loading capacity of the hydraulic jacks is 1000 kN. The loading errors are 1% of the loading capacity.  Figure 13b shows the horizontal loading devices which are two hydraulic servo actuators. One end of the hydraulic servo actuators is fixed into the actuator wall by four 100 mm-diameter bolt bars. The other end of the hydraulic servo actuators connects with the physical substructure model, in which one is located at the beam-column joint of the first story while the other at the boundary of the physical substructure model. The force capacity of two hydraulic servo actuators is ±1000 kN, and the displacement capacity is ±100 mm. The errors of both force and displacement are the 0.1% of the corresponding capacity. To avoid the horizontal slip of the physical substructure model, two measures were designed: one is that the mudsill of the physical substructure model is fixed on the floor of the lab by using 100 mm-diameter bolts spaced on 1 m centers; and the other is that two hydraulic jacks are used on both sides of the mudsill to apply the equal and opposite forces ( Figure  14). In addition, four steel rollers were used on the side of both the first beam and the truss at the boundary, to avoid the out-of-plane deformation of the steel frame ( Figure 15).    Figure 13b shows the horizontal loading devices which are two hydraulic servo actuators. One end of the hydraulic servo actuators is fixed into the actuator wall by four 100 mm-diameter bolt bars. The other end of the hydraulic servo actuators connects with the physical substructure model, in which one is located at the beam-column joint of the first story while the other at the boundary of the physical substructure model. The force capacity of two hydraulic servo actuators is ±1000 kN, and the displacement capacity is ±100 mm. The errors of both force and displacement are the 0.1% of the corresponding capacity. To avoid the horizontal slip of the physical substructure model, two measures were designed: one is that the mudsill of the physical substructure model is fixed on the floor of the lab by using 100 mm-diameter bolts spaced on 1 m centers; and the other is that two hydraulic jacks are used on both sides of the mudsill to apply the equal and opposite forces ( Figure 14). In addition, four steel rollers were used on the side of both the first beam and the truss at the boundary, to avoid the out-of-plane deformation of the steel frame ( Figure 15).  Figure 13b shows the horizontal loading devices which are two hydraulic servo actuators. One end of the hydraulic servo actuators is fixed into the actuator wall by four 100 mm-diameter bolt bars. The other end of the hydraulic servo actuators connects with the physical substructure model, in which one is located at the beam-column joint of the first story while the other at the boundary of the physical substructure model. The force capacity of two hydraulic servo actuators is ±1000 kN, and the displacement capacity is ±100 mm. The errors of both force and displacement are the 0.1% of the corresponding capacity. To avoid the horizontal slip of the physical substructure model, two measures were designed: one is that the mudsill of the physical substructure model is fixed on the floor of the lab by using 100 mm-diameter bolts spaced on 1 m centers; and the other is that two hydraulic jacks are used on both sides of the mudsill to apply the equal and opposite forces ( Figure  14). In addition, four steel rollers were used on the side of both the first beam and the truss at the boundary, to avoid the out-of-plane deformation of the steel frame ( Figure 15).    Figure 13b shows the horizontal loading devices which are two hydraulic servo actuators. One end of the hydraulic servo actuators is fixed into the actuator wall by four 100 mm-diameter bolt bars. The other end of the hydraulic servo actuators connects with the physical substructure model, in which one is located at the beam-column joint of the first story while the other at the boundary of the physical substructure model. The force capacity of two hydraulic servo actuators is ±1000 kN, and the displacement capacity is ±100 mm. The errors of both force and displacement are the 0.1% of the corresponding capacity. To avoid the horizontal slip of the physical substructure model, two measures were designed: one is that the mudsill of the physical substructure model is fixed on the floor of the lab by using 100 mm-diameter bolts spaced on 1 m centers; and the other is that two hydraulic jacks are used on both sides of the mudsill to apply the equal and opposite forces ( Figure  14). In addition, four steel rollers were used on the side of both the first beam and the truss at the boundary, to avoid the out-of-plane deformation of the steel frame ( Figure 15).

Arrangement of Displacement Measuring Points
Five linear variable differential transformers (LVDTs) were employed to measure the physical displacements of the substructure model. One of the LVDTs is placed at the mudsill to measure the slip of the physical substructure model. Two of the LVDTs, which are 300 and 1050 mm above the top of the mudsill, are placed on the column. One of the LVDTs is placed at the beam-column joint of the first story, and one of the LVDTs is placed at the boundary of the physical substructure model. The stroke of all the used LVDTs is ±150 mm. The errors of the LVDTs are 0.1% of the corresponding capacity. Figure 16a shows the arrangement of strain gauges. Each column has five sections which are chosen to paste the strain gauges. Both the top and bottom of the first column has two sections, respectively. The bottom of the second column has one section. Each end of the first beam has two sections which are chosen to paste the strain gauges. Each section has eight strain gauges to calculate the bending moment and axial force of the section (see Figure 16b).

Arrangement of Displacement Measuring Points
Five linear variable differential transformers (LVDTs) were employed to measure the physical displacements of the substructure model. One of the LVDTs is placed at the mudsill to measure the slip of the physical substructure model. Two of the LVDTs, which are 300 and 1050 mm above the top of the mudsill, are placed on the column. One of the LVDTs is placed at the beam-column joint of the first story, and one of the LVDTs is placed at the boundary of the physical substructure model. The stroke of all the used LVDTs is ±150 mm. The errors of the LVDTs are 0.1% of the corresponding capacity. Figure 16a shows the arrangement of strain gauges. Each column has five sections which are chosen to paste the strain gauges. Both the top and bottom of the first column has two sections, respectively. The bottom of the second column has one section. Each end of the first beam has two sections which are chosen to paste the strain gauges. Each section has eight strain gauges to calculate the bending moment and axial force of the section (see Figure 16b).  Figure 17 shows the comparison curves of time history between the test results and the simulation results at the first story for different levels of PGA. The results in Figure 17 were obtained by simulating the entire structure in OpenSEES, based on the parameters in Section 5.1. The HS test results are close to the simulation results at the different PGA levels. Figures 18 and 19 show the comparison curves of the internal force time history between the test results and the simulation results with the 310 and 620 cm/s 2 , respectively. Note that, for PGA levels above 620 cm/s 2 , the response of the HS test model is significantly nonlinear, resulting in failure of most of the strain gauges; as a result, the internal forces calculated using the stain response are not available above a PGA of 620 cm/s 2 .  Figure 17 shows the comparison curves of time history between the test results and the simulation results at the first story for different levels of PGA. The results in Figure 17 were obtained by simulating the entire structure in OpenSEES, based on the parameters in Section 5.1. The HS test results are close to the simulation results at the different PGA levels. Figures 18 and 19 show the comparison curves of the internal force time history between the test results and the simulation results with the 310 and 620 cm/s 2 , respectively. Note that, for PGA levels above 620 cm/s 2 , the response of the HS test model is significantly nonlinear, resulting in failure of most of the strain gauges; as a result, the internal forces calculated using the stain response are not available above a PGA of 620 cm/s 2 . The shear force in Figures 18 and 19 can be obtained conveniently from the force sensor of the horizontal hydraulic actuators. The bending moment and axial force can be calculated from the stain gauges of the calculated section by the following equation: The shear force in Figures 18 and 19 can be obtained conveniently from the force sensor of the horizontal hydraulic actuators. The bending moment and axial force can be calculated from the stain gauges of the calculated section by the following equation:

Test Results and Analysis
where ξ i is measured by the stain gauge at the section, ξ 1 and ξ 8 are two compressed strain gauges at the flange of the section, ξ 4 and ξ 5 are two tensile strain gauges at the flange of the section, E is the elasticity modulus, I is the sectional inertia moment, A is the sectional area, y is the distance between the point for the stain gauge to the center of the section, N is the axial force on the section, and M is the bending moment. Because the measure of the stain gauge is only recorded at the step when the vertical force changes, Figures 18 and 19 give the calculated values of the bending moment and the axial force at the recorded step. Figures 18 and 19 also show that the test results of internal forces are close to the simulation one. Table 3 gives the relative errors at the peak of the internal force. In Table 3, the maximum relative error is less than 10%, showing good performance of the proposed boundary technique based on the inflection point. where ξi is measured by the stain gauge at the section, ξ1 and ξ8 are two compressed strain gauges at the flange of the section, ξ4 and ξ5 are two tensile strain gauges at the flange of the section, E is the elasticity modulus, I is the sectional inertia moment, A is the sectional area, y is the distance between the point for the stain gauge to the center of the section, N is the axial force on the section, and M is the bending moment. Because the measure of the stain gauge is only recorded at the step when the vertical force changes, Figures 18 and 19 give the calculated values of the bending moment and the axial force at the recorded step. Figures 18 and 19 also show that the test results of internal forces are close to the simulation one. Table 3 gives the relative errors at the peak of the internal force. In Table 3, the maximum relative error is less than 10%, showing good performance of the proposed boundary technique based on the inflection point.   where ξi is measured by the stain gauge at the section, ξ1 and ξ8 are two compressed strain gauges at the flange of the section, ξ4 and ξ5 are two tensile strain gauges at the flange of the section, E is the elasticity modulus, I is the sectional inertia moment, A is the sectional area, y is the distance between the point for the stain gauge to the center of the section, N is the axial force on the section, and M is the bending moment. Because the measure of the stain gauge is only recorded at the step when the vertical force changes, Figures 18 and 19 give the calculated values of the bending moment and the axial force at the recorded step. Figures 18 and 19 also show that the test results of internal forces are close to the simulation one. Table 3 gives the relative errors at the peak of the internal force. In Table 3, the maximum relative error is less than 10%, showing good performance of the proposed boundary technique based on the inflection point.

Conclusions
This paper proposes a new boundary technique based on the inflection point, where the bending moment is zero, to conveniently realize the equilibrium and compatibility at the boundary. The axial displacement prediction technique is used to meet the equilibrium and compatibility of the axial direction at the boundary. To realize the process of the HS based on the proposed boundary technique, the HS component interaction is used to build a bridge between the Matlab and the OpenSEES, and the equivalent force control method is employed for solution of the EOM. Comparison of the results of the numerical simulation and the physical test show that the proposed method has the better performance than the traditional boundary technique.
Fifteen numerical examples are presented to verify the accuracy and feasibility of the proposed boundary technique. The numerical simulation results show that the proposed technique based on the inflection point can improve tremendously the test accuracy over the traditional technique.
The systematic laboratory hybrid simulation tests of the six-story steel frame structure model were designed to verify the accuracy and feasibility of the proposed boundary technique. The displacement-force mixed control technique is proposed to realize the real hybrid simulating test based on the proposed boundary technique. The tests results show reliable and accurate performance of proposed boundary technique.
Author Contributions: All authors made substantial contributions to this study. Z.C. provided the concept and design of the study. X.Y. and X.Z. performed the experiment and collected the data. H.W. and B.F.S. wrote and revised the manuscript. Z.C., X.Y. and X.Z. analyzed the experimental data.

Funding:
The research is financially supported by the Natural Science Foundation of China (51678199, 51578151, and 51161120360).