Numerical Simulation and Dynamic Analysis of Single-Hole Cliff- Side Loess Cave Dwelling under Seismic Actions

Loess cave dwellings are the most typical style of regional architecture in northwest China; now, there are still tens of millions of people living in them. The northwest is an earthquake-prone area, and cave dwellings have suffered a lot of damage in previous moderate and strong earthquakes, so their earthquake resistance has attracted people’s attention. At present, the seismic analysis of aboveground building structures is relatively mature, while the seismic analysis of loess cave dwellings is less researched. To study the seismic response of loess cave dwellings, a single-hole cliff-side loess cave dwelling located in Yangjialing revolution former sites of Yan’an City of northwest China was investigated and surveyed; the three-dimensional numerical model was established by MIDAS/GTS NX. Combining the historic earthquake damage investigation, dynamic time-history analysis of the single-hole loess cliff-side cave dwelling subject to four horizontal earthquake actions was conducted to determine the weak positions, failure characteristics, and the corresponding displacement and stress of the loess cave dwelling under earthquake load. The results show that the loess has an amplification effect on the seismic waves, the arch vault is a key factor to the stability of the loess cave dwellings, the cliff-side loess cave dwellings in an 8-degree area cannot be used to continue living, and the entrance to loess cave dwellings is the most dangerous place when the earthquake happens.


Introduction
On the loess plateau of northwest China, as influenced by economic conditions, climatic environment, and geological features, as a kind of vernacular architecture, loess cave dwelling has been sprung up greatly here and even became the main form of rural dwellings in northern Shaanxi. Loess cave dwelling is a category of architecture which gives a priority to the local requirements and make full of local building materials (raw soil) and climatic conditions [1,2]. Because of simple construction, warmth in winter and coolness in summer, and low cost, it ran rather well in the condition of dry and less rainy weather at that time, and thus, it was popular with local residents. However, in recent years, the climatic conditions of Yan'an have been changed; espe-cially in summer and autumn, the precipitation is increased sharply; the continuous rainfall causes the cave damaged or even collapsed, which makes a serious threat to the lives and property safety of the occupants, especially for those loess cave dwellings built on the unfavorable geological conditions, such as loess slope and gully zone. In addition, the loess cave dwelling itself also has the disadvantages of inadequate lighting, poor ventilation, poor moisture-proof and poor waterproof, and so on. With the improvement of people's living conditions and the progress of urbanization, people moved from cave dwellings one after another, and cave dwelling with regional culture characteristics is gradually withdrawing from the historical stage; traditional residential buildings are facing the crisis of gradual disappearance. Anyhow, some existing cave dwellings lying in the historical revolution former sites of Yan'an in northern Shaanxi province have important political value and economic significance; these cave dwellings are revolution-related historical relics and Red Tourism's major attraction,the protection of these loess cave dwellings has been getting increasing attention. In particular, northwest is a multiearthquake area, the Haiyuan earthquake with a magnitude of 8.5 and the Gulang earthquake with a magnitude of 8 occurred in this area, and the loess cave dwellings had suffered a lot of damage in previous moderate and strong earthquakes. And there are still more than 40 million people living in cave dwellings now [3]; strengthening the study on earthquake resistance of loess cave dwellings is helpful in reducing earthquake damage of loess cave dwellings and can provide theoretical reference for improving the safety of the raw soil cave structure and further protection.
There are extensive studies and efforts on the loess cave dwellings up to now, but emphases have been placed on the thermal environment [1,[3][4][5][6][7], architectural form, culture [8][9][10][11], construction techniques, and structural characteristics of cave dwellings [12,13]. However, the research on the seismic performance of the loess cave dwellings is carried out relatively few. Yang and Yang [14] summarized various kinds of earthquake damage characteristics in different regions, evaluated the seismic performance of the cliff-side and separated loess cave dwellings in the moderate strong earthquake, and then proposed some measures to improve the seismic capacity of the loess cave dwellings. Luo [15] summarized cave dwelling buildings and the earthquake damage characteristics, analyzed the seismic stability of the cave and presented the corresponding calculation formula according to the static theory, and then recommended the size of loess cave building in the earthquake zone and seismic resistant measures. Wei et al. [16] conducted the earthquake damage investigation on Shanxi DaYiyang District, summarized and analyzed the earthquake damage characteristics and destroying mechanism of the loess cave dwellings, and then evaluated the seismic performance and application scope of the loess cave dwellings. Chen et al. [17] summarized the historical earthquake damage experience of the cave dwellings on the loess precipice, proposed the estimation method of the critical height of cliff slope in the earthquake, and conducted a preliminary study on the cliff-side cave's seismic performance. Shi et al. [18] elaborated the earthquake damage characteristics and mechanism of the loess cave dwellings, made statistics on the results of earthquake damage under different seismic intensities, and analyzed the earthquake damage prediction results of various cave dwellings under different earthquake intensities. The above researches are mainly carried out based on the results of the earthquake damage investigation, which is not systematic and thorough enough on the force and deformation characteristics of cave dwellings under earthquake action. Gao and Ren [19,20] applied the finite element method, respectively, according to the plane strain problem (using 8-node and other parametric finite elements) and the space 3D problem (using 20-node and other parametric finite elements) and conducted the time-history response analysis on the loess cave dwelling under various kinds of earthquake loads with the self-compiled program; this method is not very practical and difficult to popularize. Guo and Tong [21] and Guo [22] used FLAC software to study the influence of crosswise and longitudinal earthquakes on the stability of the loess cave dwellings, but the analysis results were not mutually verified with seismic damage data, and the conclusions are worthy of discussion. At present, the seismic analysis of aboveground building structures is relatively mature, but from the above literature, the seismic performance of loess cave dwellings is less studied. In addition, the seismic experiments of underground structures are not easy to do; numerical simulation is a good method to evaluate the seismic performance of cave dwellings.
Based on the cliff-side loess cave dwelling located in the Yan'an revolutionary old site of northern Shaanxi, this study utilizes the geotechnical finite element software MIDAS/NX to build a three-dimensional finite element model, applies the 4-node three pyramid unit, uses the viscous element boundary to simulate the boundary of dynamic analysis, then applies the artificial seismic waves, and adopts the linear time-history analysis with direct integration to solve the seismic response. According to the finite element results, the forces and deformations of the single-hole cliff-side loess cave dwellings are analyzed under the action of earthquakes, and the weak positions and damage characteristics are determined under earthquake action, which is of great significance to earthquake prevention and disaster mitigation work and protection and reinforcement.

Calculation Method
2.1. Free Vibration Analysis. When conducting dynamic timehistory analysis, eigenvalue analysis is conducted firstly, that is to conduct free vibration analysis and then dynamic timehistory analysis. Inherent dynamic characteristics of the structure can be obtained through eigenvalue analysis, including vibration mode, self-vibration period, and vibration mode participation coefficient. In order to calculate the vibration mode shape and natural period under the condition of undamped free vibration, the applied characteristic equation is as follows: where ½K is the stiffness matrix of the structure, ½M is the mass matrix of the structure, ω n is the eigenvalue of the n th vibration mode, and φ n is the vibration mode vector of the n th vibration mode.

Dynamic Time-History Analysis.
Time-history analysis refers to the process of calculating the structural response (displacement, internal force, etc.) at any time according to the dynamic characteristics of the structure and the solution of the dynamic equilibrium equation when the structure is subjected to dynamic load. The dynamic equilibrium differential equation adopted in the time-history analysis is where ½M and ½K are the same as the above, ½C is the damping matrix, and pðtÞ represents dynamic load, pðtÞ 2 Geofluids = -½M € u g ðtÞ, € u g ðtÞ is time-history for inputting seismic acceleration and € uðtÞ, _ uðtÞ, and uðtÞ, respectively, are the relative displacement, velocity, and acceleration of the structure at any moment.
There are many methods to represent damping, but modal damping is usually used in numerical analysis. Modal damping is generally divided into proportional damping and nonproportional damping. The proportional damping in MIDAS soft includes mass proportional damping, stiffness proportional damping, and Rayleigh damping; Rayleigh damping was chosen in this paper. To reduce the uncertainty of stiffness proportional damping on a higher mode, the sum of mass proportional damping and stiffness proportional damping is used as the damping matrix, that is, where a 0 and a 1 are the mass and stiffness proportional factors for calculating the damping matrix; they are calculated according to the following equation: where ξ 1 and ξ 2 are, respectively, damping ratios of two different vibration frequencies ω 1 and ω 2 ; this paper takes ξ 1 = ξ 2 = ξ. The mode superposition method and direct integration method are two methods of time-history analysis. The direct integration method that sets time as an integral parameter is a method for finding the solution of the dynamic equilibrium differential equation and is used here for finding the solution. The Newmark method with better convergence is utilized for direct integration in this paper; the basic assumption is From Formula (7), we can obtain In the above formula, δ and α are constant. Under the action of dynamic load, the dynamical equilibrium differential equation of the system at t + Δt moment is Substituting Equations (8) and (6) into (9), it can be simplified as where k 1 = 1/αΔt 2 , k 2 = δ/αΔt, k 3 = 1/2α − 1, k 4 = 1/αΔt, k 5 = ðδ/2α − 1ÞΔt, k 6 = δ/α − 1, α and δ are Newmark integral constants, and Δt is the integral time step. When α = 0:25, δ = 0:5, Δt ≤ T p /10 (T p is the minimum period), the Newmark integral is unconditionally stable, and the result can reach the accuracy requirement. At this time, Equation (8) can be simplified as Let k 7 = Δtð1 − δÞ, k 8 = δΔt, then Equation (6) can be expressed as By solving Equation (10), uðt + ΔtÞ can be obtained; then, putting it in Equations (11) and (12), respectively, € uðt + ΔtÞ and _ uðt + ΔtÞ can be calculated.

Three-Dimensional Finite Element Model and Loess
Material Parameters. Although cave dwellings have various forms, they can be classified into three basic types in terms of layout and structural forms, as shown in Figure 1 [11]. Loess cave dwellings formed by digging holes in the natural cliff are called cliff-side cave dwellings, see Figure 1(a), which are more complicated than the other two in terms of force. Field measurements on cliff-side cave dwellings have been performed in the Yangjialing revolution former site of Yan'an City (see Figure 2); Figure 3 shows the geometrical parameters of the cliff-side cave dwellings. The measured geometrical dimension parameters of a single-hole cliff-side cave dwelling are shown in Table 1. The calculated region which is cut from the semi-infinite space body should be the region with a hole diameter not less than 3-4 times selected along the cavern in all directions [23]. Therefore, the bottom soil in this paper's calculation model is selected 3 times the total height of the cave dwelling, that is, 10 m, up to the earth's surface; more than 3 times the net span of the cave dwelling is selected at the left and right sides; that is, the total width of the left and right sides is 20 m; more than 3 times the depth of the cave dwelling is selected at the front and back sides of the longitudinal soil body of the loess cave dwelling; that is, the length of the soil in front of the cave dwelling is 20 m; the depth of rear soil body is 27 m. The soil body is divided into tetrahedral elements by using the automatic division entity grid method; that is, a triangular pyramid element of four nodes is used. The calculation model is shown in Figure 4. In this paper, loess is defined as an ideal elastic-plastic material, using the Mohr-Coulomb failure criterion, and the determined loess material parameter values [24,25] are shown in Table 2.

Boundary Condition.
Before dynamic time-history analysis, structural eigenvalue analysis is conducted firstly to obtain the inherent frequencies of the first and second modes which are used to compute the damping matrix required for the time-history analysis. When performing eigenvalue analysis, the boundary condition needs to be defined as the elastic boundary; MIDAS/GTS NX utilizes ground curved spring to generate the elastic boundary. To simulate bedrock condition, the bottom of the model is set to fixed condition, which is set by selecting the option of "fixed bottom condition." As conducting dynamic analysis, because the seismic wave will generate a reflection wave on the boundary, the treatment of boundary has a great influence on the analysis result. There are many methods for defining dynamic analysis boundary, mainly including artificial boundary, transfer boundary/superelement, infinite element, boundary element, and other methods. This paper adopts the viscous boundary condition. In order to define the boundary of the model as a viscous boundary, it is necessary to calculate the damping value in all directions of the rock and soil. MIDAS/GTS NX can use ground curved springs to generate viscous boundaries.

Geofluids
The parameters of the viscous boundary elements are automatically calculated according to the properties of the input rock and soil during modeling, and no additional parameters are required.

Application Method of Earthquake
Wave. Due to the lack of earthquake observation records near the local underground structure, this paper adopts the most typical EI-Centro earthquake record of Empire Valley in the United States in 1940 to conduct analysis. Figure 5 is the northsouth component acceleration curve of the EI-Centro earthquake record, the duration time is 53.72 seconds, and the maximum peak acceleration NS component is 0.3569 g, which occurred in 2.14 seconds. Here, the maximum accelerations of the input artificial earthquake motion are, respectively, 0.05 g (6-degree), 0.10 g (7-degree), 0.15 g (7-degree), and 0.2 g (8-degree); the acceleration coordinate peak of the EI-Centro wave is multiplied by the appropriate constant to modulate by the proportional method, so as to fully meet the requirements of maximum acceleration. As for calculation time, intercept the first 30 seconds, input it from the bottom along the horizontal direction, and only consider the earthquake acceleration in the input horizontal direction.

Finite Element Result Analysis and Discussion
4.1. Free Vibration Analysis. Before the dynamic timehistory analysis, the mode shape analysis of the structure is carried out first, and the characteristic period of the top ten modes is obtained. The free vibration period of the first and second modes with the largest contribution was calculated as T 1 = 0:4682564 s and T 2 = 0:3997093 s. The damping ratio of loess ranges from 0.024 to 0.26 [26,27]. Here, the damping ratio of the two modes is selected as 0.05.

Seismic Response Analysis.
In this paper, the linear timehistory analysis of the direct integration method is adopted, and the defined time step is used to perform the analysis according to the Newmark implicit integral method. The accuracy of the result is related to the size of the defined time step. Generally, the correct result can be obtained if the time step is less than 10% of the minimum period. It can be known from eigenvalue analysis that the minimum period of each mode is 0.242 s, take time step Δt = 0:01 s to meet the requirements, the duration time is 30 s, and there is a total of 3000 time step calculation. In order to investigate the seismic response of the cliff-side cave dwelling, there are in total five key points of arch vault 1, left arch camber 2, right arch camber 3, left arch foot 4, and right arch foot 5, which are selected in the position of the cave dwelling entrance, as shown in Figure 6. As calculating the horizontal displacements time-history and horizontal acceleration time-history of the five key points under four earthquake motions, it is found that under four earthquake motions (0.05 g, 0.1 g, 0.15 g, and 0.2 g), the curve variation tendencies of horizontal displacements timehistory diagram and acceleration time-history diagram of these key points are similar (the horizontal displacement time-history diagram and horizontal acceleration timehistory diagram of these key points under 0.05 g earthquake motion are given here, as shown in Figures 7 and 8). Under different horizontal earthquake motions, the horizontal displacement time-history curve and the horizontal acceleration time-history curve of arch vault 1 are the same as the input earthquake acceleration curve in the aspect of variation tendency; that is, the whole shape, the peak occurrence number, and the continuous range are very similar to the earthquake wave; only the peak occurrence moment is slightly delayed. The maximum acceleration of the input ground motion occurs when t = 2:14 s, while the curve peak   However, the variation tendency difference of the horizontal displacement time-history curve in the arch cambers and feet with the input earthquake acceleration curve before 0.1 s is very great; there is a rapid increase effect, but the peak number and shape of the subsequent curve with earthquake wave curve are very close; the displacement peak time at these four places lags behind, which occurs when t = 2:18 s, while the acceleration time-history curve with the input earthquake motion acceleration curve in the cambers and feet is quite different; the accelerations at these four places increase sharply at the initial moment, which, respectively, reaches to peak value when t = 0:10 s and t = 0:11 s. The rea-son may be that the shape of the cave arch cambers and feet change abruptly, occur corner; in addition, the cavern is empty, resulting in the rapid increase of the horizontal displacement and the horizontal acceleration in the arch cambers and feet at the initial moment. Table 3 is the horizontal displacement and vertical displacement of all key points under the gravity function. It can be seen from the number in the table that the vertical displacement of the arch vault is maximal under static load; there is no need to investigate horizontal displacement, while the horizontal displacement of arch cambers and feet must be considered. These conform to the arch vault area of cliff-side cave dwellings where it is easy to occur crack, collapse, and cave wall landslide and generate destruction phenomenon of stress concentration. Figures 9 and 10 are the maximum horizontal displacement and peak horizontal acceleration values of these key points under four earthquake motions. It can be seen that under different earthquake motions, with the increase of the inputting maximum acceleration value, the maximum horizontal displacement and the peak horizontal acceleration of all key points are also increasing continuously. When inputting the maximum accelerations of the earthquake wave which is 0.05 g, which is equivalent to the earthquake motion with the intensity of 6-degree, the horizontal displacement of all key points is equivalent to horizontal displacement under the action of gravity, which indicates that under this earthquake motion, the cave dwelling will not be destroyed. When inputting the maximum accelerations of the earthquake wave which is 0.1 g (7-degree), 0.15 g (7-degree), and 0.2 g (8-degree), the maximum horizontal displacement increases gradually, which indicates that with the increase of earthquake motion, the cave dwelling is changed from the emerging of some cracks to increasing crack width, then generating plastic zone in a wide range, eventually resulting in collapse. According to the statistical earthquake damages of the Haiyuan earthquake disaster, as shown in Table 4, the higher the seismic intensity, the more serious the earthquake damage of the cliff-side cave dwelling. The calculated results of Table 3 and Figures 9 and 10 are consistent with the investigation results of the Haiyuan earthquake disaster. Meanwhile, Figure 9 also shows that the horizontal displacement increases gradually from feet to cave vault, which is in symmetric distribution on the left arch camber and right arch camber. This is because the horizontal displacement response of the bottom of the cave is the smallest due to the constraint effect of the surrounding boundary, and the displacement gradually increases as the position moves upward, while the horizontal displacement of the top free surface is the largest.
As the maximum acceleration of seismic wave is the input of 0.05 g (6-degree), 0.1 g (7-degree), 0.15 g (7-degree), and 0.2 g (8-degree), the maximum horizontal acceleration  times of the input seismic acceleration peak value. The maximum horizontal acceleration on the arch cambers and feet also increases compared with the peak acceleration of the input seismic wave. It can be seen that the loess has an amplification effect on the seismic waves. Table 5 lists the maximum principal stress value and extreme value under the action of self-weight and seismic motion when t = 2:2 s. It can be seen that the maximum principal stress of arch cambers and feet under the action of self-weight is the compressive stress, while under the action of seismic motion, the maximum principal stress of these places occurs the tensile stress. Under the four seismic motions at the same time, there is obvious stress concentration at these concave angle places of the cave dwelling, such as the feet and arch camber places, the maximum principal stress of these parts occurs tensile stress, and the value is larger. In general, when all of these tensile stresses exceed the tensile strength of the loess, it will cause the fracture of the arch camber and the peeling of the soil at the feet.
Under the action of various seismic motions, the distribution regularities of the Von Mises strain nephogram at the cave dwelling entrance area are the same, so Figures 11  and 12 only show the Von Mises strain nephogram of the cave entrance area under the self-weight and the action of 0.05 g seismic motion. It can be seen from Figures 11 and  12 that the plastic zone of the cave is mainly located in two feet and the arch shoulder parts on both sides of the arch camber, and it has a tendency to gradually extend towards the arch vault. The weak part of the cave dwelling is the two feet and the arch camber at the cave entrance area, there is partial stress concentration for these parts, so it is easy to produce upward stretching diagonal cracks from the bottom of the cave entrance, which will cause the upper soil to slide and collapse along the cracks, resulting in the destruction to the front section of the cave entrance. Figure 13 shows the Von Mises strain nephograms of the cave's various cross-sections, under the action of 0.1 g seismic motion. Figure 14 shows the longitudinal Von Mises strain nephogram along the cave dwelling under the action of 0.1 g seismic motion. As can be seen from Figures 13  and 14, the plastic zone area of the cave dwelling entrance area is the largest, and the deeper into the cave, the plastic zone of each cross-section of the cave decreases gradually. It can be seen that as the earthquake happens, the cave entrance section is the most dangerous place. This is because there is an empty free surface in the front of the cave, so its ability to resist earthquakes is poor; as a result, it is very easy for the soil to collapse at the face of the cave dwelling.   Table 4: Earthquake damage characteristics of cliff-side loess cave dwelling in Haiyuan earthquake.

Earthquake intensity
Characteristics of earthquake damage 6-degree area Occasional collapse 7-degree area Usually, only partial areas of the cave dwelling drop earth, and individual cave dwelling faces collapse 8-degree area About 10% of cave dwelling faces collapse or is under bad damage 9-degree area More than 30% (someone thinks it is 50%) of cave dwelling faces collapse 10-degree area Almost all cave dwelling faces collapse; 30%-40% of cave dwelling vaults collapse completely

Conclusions
In this study, we explored the numerical modeling of the single-hole cliff-side loess cave dwelling dynamic timehistory response to seismic motions, determined the horizontal displacement, horizontal acceleration, stress and strain, the weak positions, failure characteristics of the loess cave dwellings under four horizontal earthquake actions with Midas/GTS geotechnical finite element analysis software. The finite element analysis results show the following. (1) Compared with the input seismic wave peak acceleration, the maximum horizontal acceleration of the arch vault, arch cambers and the feet presents an increasing phenomenon. It can be seen that the loess has an amplification effect on seismic waves. (2) Under different seismic motions, as the maximum acceleration of input increases, the horizontal displacement and horizontal acceleration peak of all key points are also increasing, and the horizontal displacement of the arch vault is larger than the feet of loess cave dwelling. This shows that the arch vault is a key factor in the stability of cave dwellings.
(3) When subjected to the seismic motion of 8 degrees, the cave dwelling can gradually generate plastic areas in a larger range, even results in collapse, so cave dwelling in the 8degree area cannot be used to continue living. (4) Furthermore, the plastic area of the cave is mainly located in two feet and the arch shoulder parts on both sides of the arch camber, the weak part of the cave is the two feet and the arch camber of the cave dwelling. As the earthquake happens, the plastic area of the cave dwelling entrance area is the largest, the cave dwelling entrance section is the most dangerous place. (5) The established model dimension is reasonable, the boundary conditions are properly handled, the method of applying seismic wave is correct, and the results of the dynamic time-history analysis are consistent with the investigation results of earthquake damage. This research method can be extended to the analysis of the seismic response of porous cave dwellings.

Data Availability
All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.  12 Geofluids