Numerical Algorithm for Dynamic Impedance of Bridge Pile-Group Foundation and Its Validation

The characteristics of bridge pile-group foundation have a significant influence on the dynamic performance of the superstructure. Most of the existing analysis methods for the pile-group foundation impedance take the trait of strong specialty, which cannot be generalized in practical projects. Therefore, a project-oriented numerical solution algorithm is proposed to compute the dynamic impedance of bridge pile-group foundation. Based on the theory of viscous-spring artificial boundary, the derivation and solution of the impedance function are transferred to numerical modeling and harmonic analysis, which can be carried out through the finite element method. By taking a typical pile-group foundation as a case study, the results based on the algorithm are compared with those from existing literature. Moreover, an impact experiment of a real pile-group foundation was implemented, the results of which are also compared with those resulting from the proposed numerical algorithm. Both comparisons show that the proposed numerical algorithm satisfies engineering precision, thus showing good effectiveness in application.


Introduction
Pile-group foundation has been widely used in bridge engineering, and pile-soil interaction has a significant effect on the dynamic performance of the superstructure [1][2][3][4]. Mylonakis et al. [5] investigated the influence of soil-structure interaction (SSI) on the seismic response of structures by using recorded motions, and concluded that some simplifications about SSI in design may lead to detrimental consequences. Soyluk et al. [6] studied the effects of SSI and spatially varied ground motions on the seismic response of a cable-stayed bridge, and the results indicated that softer soils amplify the displacement response of the bridge. Shirgir et al. [7] developed an analytical model to obtain free vibration frequencies of bridges with soil-structure interaction taken into consideration, and found that the effect of SSI is significant for soils with a shear wave velocity less than 175 m/s. Gara et al. [8] discussed the contribution of SSI to the bridge dynamics, and found that SSI may significantly affect the structural dynamics, especially in the case of soft or medium soil deposits. Qiao et al. [9] established an analytical model of the train-bridge system considering SSI based on the substructure method, and found that it is important to consider SSI in the dynamic analysis of train-bridge systems, but such influence becomes weaker with the increase in the soil shear wave velocity. Tochaei et al. [10] investigated the effects of soil-structure interaction and near-field seismic ground motions on cable-stayed bridges based on experiments and numerical analysis; their results indicated that the effect of foundation soil stiffness on the response of the bridge was influenced by the type of input ground motion and that the bridge response was amplified when their foundation was in softer soils. All of the above studies indicate the importance of considering SSI in the dynamic analysis of bridges.
There are two main approaches for soil-structure analysis during the dynamic analysis of bridges-namely, the direct approach, and the substructure approach [11]. Compared with the direct approach-which studies the structure, foundation, and soil simultaneously (usually by means of a finite element model)-the substructure approach, which divides the whole system into two subsystems (superstructure subsystem and soil-foundation subsystem), is less time consuming and more efficient, for the two subsystems can be treated separately with their respective optimal method. Based on the substructure method, the coupling of two interacting subsystems is realized through the concept of the soilfoundation dynamic impedance functions. Foundation impedance functions represent the dynamic stiffness of the soil medium surrounding the foundation, which characterize the compliant base restraints for the superstructure and the input motion experienced by the foundation, thus providing a simple means of accounting for pile-soil interaction [12]. Therefore, the study of the solution of dynamic impedance functions of pile-group foundation has drawn more and more attention.
The generally accepted standard solutions for the dynamic response of pile groups in viscoelastic media were obtained by Kaynia and Kausell, by means of the boundary integral method [13]. Afterwards, Gazetas et al. [14] proposed a method to calculate the impedance functions of pile-group foundations using a simplified beam-on-dynamic-Winkler-foundation model and Green's function-based rigorous method. Based on the thin layer method and the flexible volume method, Jiang et al. [15] calculated the impedance functions of pile-group foundations from literature as well as in a practical layered site, and compared the results with those in the literature [16] and from experimental data. In recent years, more and more methods have been proposed for the impedance functions of pile-group foundations under various conditions [17][18][19][20][21][22][23]. For example, Padrón et al. [24] presented a BEM-FEM model for the time-harmonic dynamic analysis of piles and pile groups embedded in an elastic half-space, and successfully reduced the number of degrees of freedom. Emani and Maheshwari [25] developed a general three-dimensional finite element procedure for the dynamic impedances of pile groups with caps embedded in isotropic homogeneous elastic soils. Ai et al. [26] investigated the dynamic behavior of a laterally loaded fixed-head pile group in a transversely isotropic multilayered half-space by BEM-FEM coupling method. With the porous medium model and the Rayleigh-Love rod model describing the soil and the pile groups, respectively, Zhang et al. [27] developed a simplified analytical solution for the vertical impedance function of viscoelastic pile groups partially embedded in a layered, transversely isotropic, liquid-saturated, viscoelastic soil medium. Dezi et al. [28,29] proposed a numerical model for the 3D kinematic interaction analysis of vertical and inclined pile groups in horizontally layered soils, based on BEM-FEM coupling, in which piles are modelled with finite beam elements while the soil is schematized with infinite independent horizontal layers. Moghaddasi et al. [30] presented a hybrid analytical-numerical method for the dynamic analysis of piles and pile groups embedded in nonhomogeneous transversely isotropic media subjected to lateral excitations, in which the piles are simulated using finite elements and the interactions between piles and soil are captured through a system of rigid, massless elements. Varghese et al. [31] proposed a substructure-method-based soil-structure interaction analysis model to study the vertical and horizontal dynamic impedances of piled rafts in homogeneous and layered soil conditions. Wang and Ai [32] developed a method to compute the dynamic response of pile groups in layered soils with compressible constituents by combining the analytical layer element method and finite element method. Zhang et al. [33] presented an analytical method to study the horizontal dynamic response of pile groups embedded in partially saturated soil based on the motion theory of three-phase poroelastic media, with the pilesoil system modeled using the Timoshenko beam theory and soil resistance incorporated using Novak's plane strain model. Almost all of the abovementioned methods possess the trait of strong specialty, the theoretical derivation procedures of which are complex and even involve programming processes. In view of the complexity of actual soil layers and the parametric diversity of bridge group foundations, it can be found that the above calculation methods are too complex and specialized to be generalized in practical projects. Therefore, this paper presents a project-oriented numerical algorithm for the dynamic impedance of bridge pile-group foundation. Based on the viscous-spring artificial boundary theory, the derivation and solution of the impedance function are transferred to numerical modeling and harmonic analysis, which can be carried out using the finite element method. The applicability of the proposed numerical algorithm is verified by comparing its calculation results with those obtained through methods in existing literature and a field experiment. Both comparisons show that the numerical algorithm satisfies the project demand of analysis precision, thus showing good effectiveness in application.

Definition of the Impedance Functions
For a generic bridge founded on pile groups, the soil-structure interaction problem can be conveniently studied in the frequency domain by assuming a linear behavior for soil and structure. By adopting a finite element method, the following system of complex linear equations governs the dynamics of the whole system [3]: where Z is the dynamic stiffness matrix of the soil-foundation-superstructure system; d is the vector of the nodal displacements; and f is the vector of nodal forces. With the assumption that the caps connecting the pile heads are rigid, and assigning different master nodes to them, Equation (1) can be partitioned as: in which the subscripts s, f, and e denote the superstructure, the rigid pile cap, and the piles, respectively. It is worth noting that no displacement component for the soil appears in Equation (2) because the pile-soil interaction is accounted for in the definition of the dynamic stiffness matrices relevant to the rigid caps and piles [3].
With the linear assumption, and according to the substructure approach, the inertial interaction for the superstructure is governed by: and for the soil-foundation system, the kinematic interaction problem is governed by: where d e,kin are the displacements of the embedded piles and d FIM are those of the master modes of the rigid caps due to the kinematic interaction. Then, the impedance matrix of the soil-structure system S can be expressed as: It can be found that foundation impedance functions can be defined as the ratio of a harmonic force (or a moment) applied on the foundation-soil interface to the corresponding harmonic displacement (or rotation), dependent not only on the frequency of the excitation, but also on the foundation geometry and the dynamic characteristics of the soil medium [12]. Based on this definition, a project-oriented numerical algorithm for impedance functions of bridge pile-group foundation is proposed. First, according to the actual situation of the bridge site, a refined finite element analysis model of the cap, pile group, and surrounding soil is established; second, based on the viscous-spring artificial boundary theory, a series of springs and dampers are set at the boundaries of the aforementioned model so as to simulate the influence of the infinite soil; third, a unit resonant force in each direction within the concerned frequency band is applied at the center of the cap, and the corresponding displacement response can then be obtained-namely, the dynamic flexibility matrix; finally, the impedance functions of the pile-group foundation can be obtained in virtue of the inverse properties of the dynamic flexibility matrix and dynamic stiffness matrix.

Refined FE Model of the Cap-Pile Group-Surrounding Soil System
The pile-group foundation of the bridge can be regarded as a dynamic interaction system composed of the pile cap, pile group, and layered soil, the pile group effect of which is prominent. Therefore, to better consider the influence of the pile cap, pile layout, pile spacing, layered soil, soil properties, and other factors on the dynamic performance of pile-group foundations, a three-dimensional analysis model is necessary.
Based on the finite element method, which is widely recognized and mastered by structural engineers, a refined FE model of the cap-pile group-surrounding soil system is established, as shown in Figure 1. Here, the ANSYS software [34] is used for establishing the model, in which the cap, piles, and soil medium are all modelled with Solid45 elements. The connection between the cap and the piles, and that between the piles and the soil, are realized by master-slave nodes or constraint equations. It is assumed that the interface between the piles and the surrounding soil is always in an ideal bonding state. element analysis model of the cap, pile group, and surrounding soil is est second, based on the viscous-spring artificial boundary theory, a series of spr dampers are set at the boundaries of the aforementioned model so as to sim influence of the infinite soil; third, a unit resonant force in each direction w concerned frequency band is applied at the center of the cap, and the corre displacement response can then be obtained-namely, the dynamic flexibilit finally, the impedance functions of the pile-group foundation can be obtained in the inverse properties of the dynamic flexibility matrix and dynamic stiffness m

Refined FE Model of the Cap-Pile Group-Surrounding Soil System
The pile-group foundation of the bridge can be regarded as a dynamic in system composed of the pile cap, pile group, and layered soil, the pile group which is prominent. Therefore, to better consider the influence of the pile cap, pi pile spacing, layered soil, soil properties, and other factors on the dynamic per of pile-group foundations, a three-dimensional analysis model is necessary.
Based on the finite element method, which is widely recognized and ma structural engineers, a refined FE model of the cap-pile group-surrounding soil established, as shown in Figure 1. Here, the ANSYS software [34] is used for est the model, in which the cap, piles, and soil medium are all modelled with elements. The connection between the cap and the piles, and that between the the soil, are realized by master-slave nodes or constraint equations. It is assume interface between the piles and the surrounding soil is always in an ideal bondin

Boundary Conditions
To establish the FE model of the cap-pile group-surrounding soil sys necessary to extract a finite domain from the unbounded soil medium. In accurately model the originally continuous medium after extraction, the pro characteristics of the soil should be kept unchanged, which necessitates an boundary capable of absorbing the energy of the scattering waves [35].

Boundary Conditions
To establish the FE model of the cap-pile group-surrounding soil system, it is necessary to extract a finite domain from the unbounded soil medium. In order to accurately model the originally continuous medium after extraction, the propagation characteristics of the soil should be kept unchanged, which necessitates an artificial boundary capable of absorbing the energy of the scattering waves [35]. For this purpose, the modified viscous-spring artificial boundary proposed by Du et al. [36][37][38] was adopted, since it not only satisfies the above criterion, but also has the advantages of high accuracy, adequate stability, and convenient application. The three-dimensional sketch of the adopted artificial boundary is shown in Figure 2. adopted, since it not only satisfies the above criterion, but also has the a accuracy, adequate stability, and convenient application. The three-dim the adopted artificial boundary is shown in Figure 2. To simulate the springs and dampers of the artificial bounda element with one end fixed is used in the ANSYS model; then, the F calculating the impedance functions of pile-group foundations is comp

Calculation of the Impedance Funcitons
The motion equation of the cap-pile group-soil system under har expressed as: where m, c, and k are the mass matrix, damping matrix, and respectively; u is the displacement vector; and p is the force vector who Assuming that: The parameters of the springs and dampers that constitute the artificial boundary can be obtained by the following equations: in tangential direction where G, ρ, and λ are the shear modulus, mass density, and Lamé constant of the soil, respectively; c p and c s are the velocity of P-wave and S-wave, respectively; r is the distance between the point load and the boundary; and α and β are coefficients, the recommended empirical values of which are 0.8 and 1.1, respectively. To simulate the springs and dampers of the artificial boundary, the combin14 element with one end fixed is used in the ANSYS model; then, the FE model used for calculating the impedance functions of pile-group foundations is completed.

Calculation of the Impedance Funcitons
The motion equation of the cap-pile group-soil system under harmonic load can be expressed as: m ..
where m, c, and k are the mass matrix, damping matrix, and stiffness matrix, respectively; u is the displacement vector; and p is the force vector whose frequency is ω.
Assuming that: by substituting Equation (9) into Equation (8), it can be obtained that: where: Then, u can be obtained through the decomposition of complex matrix k c . By applying the unit harmonic loads with different frequencies at the center of the cap top in the vertical, transverse, longitudinal, transverse rotation, and longitudinal rotation directions, the dynamic flexibility matrix of the pile foundation can be obtained by calculating the corresponding displacement, which can be expressed as: where the subscripts Z, Y, X, RX, and RY represent the vertical, transverse, longitudinal, longitudinal rotation, and transverse rotation directions, respectively. Finally, the impedance functions of the pile-group foundation can be obtained in virtue of the inverse properties of the dynamic flexibility matrix and dynamic stiffness matrix, as follows:

Verification of the Numerical Algorithm
In this section, the calculation results of the proposed algorithm are compared with those obtained via existing methods and field test results, in order to verify its applicability.

Comparison with Existing Methods
First, a 4 × 4 pile-group foundation was chosen as an example, the parameter values of which are taken from [15] and shown in Table 1. For comparison, both the thin layer method presented in [15] and the proposed algorithm were used to obtain the impedance functions of the single pile as well as the pile-group foundation. The overall size of the FE model in the proposed algorithm is 80 m × 80 m × 80 m, and the size of the soil element is set as 2.5 m. The damping of the piles is neglected in accordance with the benchmark. Figures 3 and 4 show the calculated vertical impedance functions of the single pile and the pile-group foundation, respectively, as obtained by the proposed algorithm and the thin layer method.    The impedance functions are expressed as k+ia0c, where k and c ar and the damping item, respectively. VV s k represents the vertical static s pile in elastic half-space, while a0 is the dimensionless frequency, which as: where ω is the frequency; and d and Vs are the pile diameter and shea the soil, respectively.
It can be seen that the results obtained by the proposed algorit The impedance functions are expressed as k + ia 0 c, where k and c are the stiffness item and the damping item, respectively. k S VV represents the vertical static stiffness of a single pile in elastic half-space, while a 0 is the dimensionless frequency, which can be expressed as: where ω is the frequency; and d and V s are the pile diameter and shear wave velocity of the soil, respectively. It can be seen that the results obtained by the proposed algorithm are consistent with those obtained by the thin layer method. Even though there is some error, which probably occurs due to the relatively large element size or the difference between calculation assumptions, from the perspective of practical application, the results are sufficiently accurate to meet the requirements of engineering.
Second, the results obtained by the proposed algorithm were compared with those obtained through a rigorous method presented in a recent study [39] based on a 4 × 4 pile-group foundation. For the results presented in this work, it was assumed that υ s (Poisson's ratio of the soil) = 0.4, E pile /E soil = 100 (E pile and E soil are the elasticity moduli of the piles and soil, respectively), υ p (Poisson's ratio of the piles) = 0.25, L (length of the piles) = 15d (where d is the pile diameter), S (pile spacing) = 5d, and ρ s /ρ p = 0.7, where ρ s and ρ p are the mass density of the soil and the piles, respectively. Figure 5 shows the calculated vertical impedance functions of the pile-group foundation, as obtained by the proposed algorithm and the existing rigorous method. It should be noted that the symbols in Figure 5 have the same meanings as those in Figure 4. It can be seen that the results obtained by the proposed algorithm are also consistent with those obtained by the rigorous method. The error of the maximum stiffness item is 0.6% and 4% as calculated from Figures 4 and 5, respectively, falling within an acceptable range of values for engineering projects. Therefore, through classic and recent numerical examples, the effectiveness of the proposed algorithm is verified and, compared with other methods that require a higher theoretical basis of users, the proposed algorithm is more friendly and easy to use for most engineers, which has certain advantages. rithms 2021, 14, x FOR PEER REVIEW Figure 5 shows the calculated vertical impedance functions foundation, as obtained by the proposed algorithm and the existing ri should be noted that the symbols in Figure 5 have the same meanings 4. It can be seen that the results obtained by the proposed algorithm with those obtained by the rigorous method. The error of the maximu 0.6% and 4% as calculated from Figures 4 and 5, respectively, acceptable range of values for engineering projects. Therefore, through numerical examples, the effectiveness of the proposed algorithm compared with other methods that require a higher theoretical b proposed algorithm is more friendly and easy to use for most eng certain advantages.

Description of the Experiment
To verify the proposed numerical algorithm, a field experiment o Tianjin was carried out. Figure 6 shows the testing site. It can construction of the pile-group foundation has been completed, while has not yet been constructed.

Description of the Experiment
To verify the proposed numerical algorithm, a field experiment on a real bridge in Tianjin was carried out. Figure 6 shows the testing site. It can be seen that the construction of the pile-group foundation has been completed, while the superstructure has not yet been constructed.
To verify the proposed numerical algorithm, a field experi Tianjin was carried out. Figure 6 shows the testing site. I construction of the pile-group foundation has been completed, has not yet been constructed.  The size of the testing cap was 7.6 m × 5 m × 2 m, and each foundation had 6 piles in total, the layout of which is shown in Figure 7. The length and diameter of each pile were 42.1 m and 1.0 m, respectively. The main kinds of soil in the bridge site were silt, mucky silty clay, silty clay, and clay. The distribution of the soil is shown in Figure 8, and the corresponding characteristics are shown in Table 2. orithms 2021, 14, x FOR PEER REVIEW The size of the testing cap was 7.6 m × 5 m × 2m, and each found total, the layout of which is shown in Figure 7. The length and dia were 42.1 m and 1.0 m, respectively. The main kinds of soil in the b mucky silty clay, silty clay, and clay. The distribution of the soil is and the corresponding characteristics are shown in Table 2.      During the test, a high-elastic cohesive hammer was used to impact the designated position of the pile cap for excitation (as shown in Figure 9), and the responses of the foundation could be acquired by installing a series of vibration sensors ( Figure 10) on the top of the cap. type of sensor selected was 941B, from the Institute of Engineering Mec sensitivity and passband are 0.3 Vs 2 /m and 0.25-80 Hz, respecti microtremors, and their sensitivity can meet the required test precision. frequency was set to 256 Hz and, under each working condition, the knocked about 20 times so that the signal could be collected.

Experimental results
The impedance functions of the pile-group foundation can be obtain the inverse properties of the dynamic flexibility matrix and dynamic stiff other words, by comparing the dynamic flexibility of the pile-grou measured experimentally with that calculated by the proposed a applicability of the latter in calculating the impedance functions foundations can be verified. Figure 11 shows the vertical and horizontal (in both north-south dire There were three working conditions in the test: (1) Impact the top of the cap in a vertical direction and obtain its vertical acceleration response; (2) Impact the north side of the cap and obtain its acceleration response in a north-south direction; (3) Impact the east side of the cap and obtain its acceleration response in an east-west direction.
The type of hammer and the data acquisition instrument were DFC-1 and INV3018A, respectively-both from the China Orient Institute of Noise & Vibration. The type of sensor selected was 941B, from the Institute of Engineering Mechanics, whose sensitivity and passband are 0.3 Vs 2 /m and 0.25-80 Hz, respectively, due to microtremors, and their sensitivity can meet the required test precision. The sampling frequency was set to 256 Hz and, under each working condition, the hammer was knocked about 20 times so that the signal could be collected.

Experimental Results
The impedance functions of the pile-group foundation can be obtained in virtue of the inverse properties of the dynamic flexibility matrix and dynamic stiffness matrix. In other words, by comparing the dynamic flexibility of the pile-group foundation measured experimentally with that calculated by the proposed algorithm, the applicability of the latter in calculating the impedance functions of pile-group foundations can be verified. Figure 11 shows the vertical and horizontal (in both north-south direction and eastwest directions) dynamic flexibility of the pile-group foundation obtained by the proposed algorithm and field experiment. There exists interference from passing trains and largescale mechanical operation around the test site, the low-frequency component of which is large. Meanwhile, in consideration of the limitations of the test equipment, the data within the range of 8-20 Hz are analyzed and compared.
It can be seen that the dynamic flexibility of the pile-group foundation in three directions obtained from the proposed algorithm is in good agreement with the results from the field experiment in the frequency range of 8-20 Hz, which once again confirms the applicability of the proposed algorithm for solving the dynamic impedance of pile-group foundations.

Experimental results
The impedance functions of the pile-group foundation can be obtaine the inverse properties of the dynamic flexibility matrix and dynamic stiffn other words, by comparing the dynamic flexibility of the pile-grou measured experimentally with that calculated by the proposed a applicability of the latter in calculating the impedance functions foundations can be verified. Figure 11 shows the vertical and horizontal (in both north-south direc west directions) dynamic flexibility of the pile-group foundation obt proposed algorithm and field experiment. There exists interference from and large-scale mechanical operation around the test site, the low-frequen of which is large. Meanwhile, in consideration of the limitations of the te the data within the range of 8-20 Hz are analyzed and compared.   It can be seen that the dynamic flexibility of the pile-group fou directions obtained from the proposed algorithm is in good agreement from the field experiment in the frequency range of 8-20 Hz, which onc the applicability of the proposed algorithm for solving the dynam pile-group foundations.

Conclusions
Existing research methods for the dynamic impedance func pile-group foundations are too complex and specialized to be general projects. Therefore, based on the viscous-spring artificial boundary theo element method, a project-oriented numerical algorithm for the impe pile-group foundations is proposed in this paper, through which the dyn function is calculated and obtained according to its definition.
The calculation results are compared with those obtained by the me in existing literature, as well as those from a field experiment. The compa the results obtained by the numerical algorithm are consistent with th existing methods (error of maximum stiffness item is less than 5% experiment, and its accuracy can meet the needs of engineerin

Conclusions
Existing research methods for the dynamic impedance functions of bridge pile-group foundations are too complex and specialized to be generalized in practical projects. Therefore, based on the viscous-spring artificial boundary theory and the finite element method, a project-oriented numerical algorithm for the impedance of bridge pile-group foundations is proposed in this paper, through which the dynamic impedance function is calculated and obtained according to its definition.
The calculation results are compared with those obtained by the methods presented in existing literature, as well as those from a field experiment. The comparisons show that the results obtained by the numerical algorithm are consistent with those obtained by existing methods (error of maximum stiffness item is less than 5%) and the field experiment, and its accuracy can meet the needs of engineering, showing the effectiveness of the proposed numerical algorithm.
The study in this paper may be beneficial for simplifying the modeling of the complex train-bridge-foundation-soil coupled system, and provide an efficient and simple means of further studying the influence of pile-soil interaction on the seismic response of a coupled train-bridge system.