Abstract: This research has been carried out with the purpose of studying wave propagation in elastic environments as well as its simulation on the basis of numerical and analytical methods. The results obtained have been used as a case study dynamic modeling the Masjed-Soleiman’s cavern. Therefore, first the relationship between the model’s optimum mesh size and the frequency and the input wave amplification was determined so that the results obtained could be used to verify the numerical model developed for analyzing the cavern dynamically. Then, strains obtained from the wave’s propagation of earthquake in Masjed Soleiman were studied under free field boundary conditions. The results obtained show that the power house cavern of Masjed Soleiman was resistant with regard to the record applied.
INTRODUCTION
Due to the recent developments in numerical modeling and computer processing, numerical modeling can evidently be used to study wave propagation. Thus, the effect of seismic wave’s propagation on underground structures can be studied through implementing the above-mentioned methods. However, it should be considered that solving dynamic problems numerically through computer processing is different from solving the static problems numerically. One of such differences concerns the verification of the model created. In static analysis, the model can be verified by making a comparison between the results of the monitoring instruments mounted with those of the model developed. However, it isn’t true concerning dynamic analysis. As a result one cannot ensure that the model propagates the wave truly (Hashash et al., 2001).
Therefore, in this research, in order to overcome this problem, there has been an attempt to develop a computer modeling of seismic wave propagation through developing wave propagation model in one-dimension as well as obtaining the necessary informations. For this purpose, by computer modeling as well as numerical resolution of one-dimensional longitudinal wave propagation along elastic environment, the relationship between the size of model meshes and input wave amplitude has been studied and by using the results obtained, the best size of mesh for modeling of the Masjed Soleiman cavern may be identified.
NUMERICAL MODELING OF ONE-DIMENSIONAL ELASTIC LONGITUDINAL WAVES PROPAGATION
To obtain further information concerning modeling of dynamic wave propagation, there has been an attempt to model the propagation of one-dimensional pressure waves through finite difference software and the relevant results should be compared with the analytical relations of the existing closed form, so that information could be obtained on the optimum (optimized) dimension’s of the model blocks for wave propagation and their relationship with dominant frequency of the input wave as well (Ma and Zhou, 1998).
Description of the created model: To examine the propagation of one-dimensional waves, one can use the propagation of wave along a bar. For modeling a bar, a rectangular cube with elastic behavior model has been considered. The physical and mechanical properties of the given bar are shown in Table 1. Input wave is a sine harmonic wave with frequency of 2.5 Hz, which has been applied to the model as time histories of longitudinal-wave velocity.
| Table 1: | Physical and mechanical properties of given bar |
| Fig. 1: | Created model and its boundary conditions |
The boundary conditions of the model have been assumed as infinite boundary conditions in order to resolve the reflection of wave from the boundaries and the bar would be assumed infinite. To examine model state when the wave passes, velocity time histories for each point (3 point across the bar) has been read and compared to each other. The model which is created and the location of 3 points are shown in Fig. 1.
Defining optimum size of model blocks: In order to identify the optimum size of model blocks and examine modification of wave propagation, the two following criteria can be controlled. The first criterion for controlling the trueness of wave propagation is one in which the plane waves are not subject to distortion after their propagation through stiff materials and shape of wave must be the same both at the beginning and at the end, without having any tangible difference. The second criterion is that the velocity of wave propagation in numerical model must be equal to the velocity calculated from expression 1 for longitudinal wave (Leif, 1975).
(1) |
| E | = | Elastic model |
| ρ | = | Environment density |
Considering these two criteria to obtain optimum mesh size, the selection of mesh size started from 200 m and by gradually making input wave constant, size of blocks decreased to obtain optimum size.
As can be shown from Fig. 2, the size of the block is so long that the numerical model cannot receive input wave and the given figure (history of velocity of 3 points across the bar) does not show true propagation of the wave. This condition continues in Fig. 3 and in the end, with a block size of 80 m (Fig. 4), wave propagation is obtained which is without distortion of the primitive wave. But, the second provision has not been fulfilled and wave velocity from numerical method is different from one obtained from analytical method. Thus, the dimensions of the blocks have still been assumed smaller and eventually, the wave velocity obtained, with mesh size of 60 m (Fig. 5) might be equal to the velocity obtained from expression 1. If the block dimensions continue to be smaller in (Fig. 6, 7) may have no effect on the truth of results and only the calculating time increases. As a result, block dimensions of 60 m made from the materials with the properties mentioned in Table 1, which have the same mechanical properties of rocks encountered around Masjed Soleiman cavern, can be assumed as optimum dimension of block. Displacement time histories of 3 points on the bar obtained from numerical analysis is compared with the displacements obtained from analytical expressions in an attempt to further investigate the trueness of longitudinal wave propagation in the model developed for two models with block sizes of 80 and 60 m.
In analytical method, time histories of displacement may be obtained by using integral calculus of velocity time histories and incorporating boundary conditions:
As can be seen from Fig. 8 and 9, the displacement obtained from the model with block size of 60 m was more congruent with analytical results compared to the one obtained from the model with block size of 80 m. This is another reason for proving that the block size of 60 m is an optimum size for a bar with the properties shown in Table 1.
The effect of input wave’s amplitude on the longitudinal wave propagation: To investigate the effect of input wave amplitude on the longitudinal wave propagation, input wave amplitude was changed. A time history of velocity at 3 points across the bar was calculated for maximum amplitude of 100, 10, 0.01 and the relevant results of which are presented in Fig. 10-13. As it can be seen, the amplitude of the wave has no effect on the propagation and distortion of wave.
| Fig. 2: | Time history of three points at the beginning middle and the end for the bar with elements of 200 m |
| Fig. 3: | Time history of three points at the beginning middle and the end for the bar with elements of 150 m |
| Fig. 4: | Time history of three points at the beginning, middle and the end for the bar with elements of 80 m |
| Fig. 5: | Time history of three points at the beginning, middle and the end for the bar with elements of 60 m |
| Fig. 6: | Time history of three points at the beginning, middle and the end for the bar with elements of 20 m |
| Fig. 7: | Time history of three points at the beginning, middle and the end for the bar with elements of 2 m |
| Fig. 8: | Comparison of displacements obtained from analytical solution and numerical modeling for a model with 80 m mesh size |
| Fig. 9: | Comparison of displacements obtained from analytical solution and numerical modeling for a model with 60 m mesh size |
| Fig. 10: | Effect of input wave altitude on its propagation for sine wave with 100 m amplitude |
| Fig. 11: | Effect of input wave altitude on its propagation for sine wave with 10 m amplitude |
| Fig. 12: | Effect of input wave altitude on its propagation for sine wave with 1 m amplitude |
| Fig. 13: | Effect of input wave altitude on its propagation for sine wave with 0.01 m amplitude |
| Fig. 14: | Effect of input wave frequency on its propagation for sine wave with 200 Hz frequency |
| Fig. 15: | Effect of input wave frequency on its propagation for sine wave with 20 Hz frequency |
| Fig. 16: | Effect of input wave frequency on its propagation for sine wave with 10 Hz frequency |
| Fig. 17: | Effect of input wave frequency on its propagation for sine wave with 5 Hz frequency |
| Fig. 18: | Effect of input wave frequency on its propagation for sine wave with 0.2 Hz frequency |
| Fig. 19: | Effect of input wave frequency on its propagation for sine wave with 0.02 Hz frequency |
The effect of input wave frequency on propagation of wave: In order to study this issue, the frequency of input wave was changed by considering the size of mesh as a constant (60) and the relevant results were taken into consideration.
In a bar with mesh dimensions of 60 m, sine wave with frequencies of 200, 20, 10, 5, 0.2 and 0.02 Hz was propagated. For each wave, velocity time histories were read at the beginning, end and middle of the bar (Fig. 1). The results obtained are presented in Fig. 14-19. It is clear that for a bar with a constant size of mesh, some distortion exist which is due to increase of input wave frequency. Consequently, there is a relationship between model mesh size and frequency of input wave for true propagation of waves. This relationship can be expressed as follows (Kuhlemeyer and Lysmer, 1973):
(2) |
| λ | = | Wave length, input wave |
| V | = | Input wave velocity in numerical model |
| f | = | Frequency of input wave |
It means that the size of model mesh must be smaller than V/10f to overcome the problem of waves distortion in the model. Also in the case of earthquake, it can be examined whether expression 2 in the geometrically model is correct by identifying dominant frequency.
DYNAMIC ANALYSIS OF MASJED SOLEIMAN CAVERN
The Masjed Soleiman dam has two parallel caverns, one is the cavern of power plant and the other is the transformer cavern (with smaller dimensions).
There are many parallel tunnels between these two large spaces which connect these two caverns (Fig. 20). In this analysis, it has been attempted to investigation the effect of maximum design earthquake loads on the Masjed Soleiman cavern. For dynamic analysis of Masjed Soleiman, the following stages have been carried out.
Geometrical model of spaces: In the primitive stage, geometry of the cavern as well as tunnels between those was created by cylindrical and cubic elements. Then, elements of the surrounding environment were developed and eventually, Fig. 21, which is a rectangular cube with 260000 blocks was developed. In selecting dimensions of elements (blocks), by using results of the earlier part, the lengths of elements were selected in a way that the number of equations would be fewer and the length of element be less than 60 m to prevent wave distortion. After generating the geometry model, the geomechanical properties of rock (Table 2) exercised to model (Moshanir et al., 1992, 1994).
Initial and boundary conditions of model and its equilibrium before excavation: One of the important stages in numerical modeling is applying suitable initial and boundary conditions to the model. The boundary conditions should be at a distance from the structure to be congruent with real conditions of environment. Before excavating, the model shown in (Fig. 21) should maintain its equilibrium.
| Fig. 20: | General view of modeled spaces in this study |
| Fig. 21: | Cubic block in which underground spaces of dam would be excavated |
| Table 2: | Geomechanical properties of rock |
In other sense, unbalance forces which have been established during creating the geometry of the model as well as boundary and initial conditions must be set to zero. This is shown in Fig. 22.
Excavation and support installation: After excavating underground spaces, primary supports were installed. The support system consists of a shotcrete layer of 15 cm with a system of anchor-bolt (Moshanir et al., 1992). For modeling shotcrete and anchor-bolts, shell element and cable elements were used respectively. In addition to steel mechanical properties, the properties of grouting have been considered in implementing this element.
Establishing static equilibrium and analysis of results: Examining of model static equilibrium as well as its validity may be performed following excavation and support installation. In order to control the validity of the model developed, the results obtained from numerical modeling have been compared with those obtained from extensometers mounted on the roof of power plant cavern (Moshanir, 2002). For this purpose, three sections at distances of 34.75, 60.75 and 85.75 m from cavern length were selected and at these three places, the rate of displacements in alignment with Z axis was investigated. The results of which have been shown in Fig. 23-25. The comparison of numerical results with those of instrumentation show a relatively proper consistency of numerical results with reality.
Dynamic analysis of model: After establishing the static equilibrium this model may be used for dynamic analysis. In dynamic analysis, the following stages should be performed (Kirzhner and Rosenhouse, 2000; Hashash et al., 2001).
Identification of suitable earthquake record for applying to model: In this research, a record of earthquake was selected and prepared by performing necessary modification, to be used in Masjed Soleiman dam. The selected earthquake record is on the basis of acceleration time histories and concerns Gilroy earthquake. One of reasons for selecting this earthquake is that the maximum altitude of this earthquake is close to that of predicted in Masjed Soleiman area. After selection of suitable earthquake record, the following modifications should be carried out:
| • | Modification of record acceleration maximum altitude: Since the estimated maximum predictable earthquake for Masjed Soleiman dam zone included 0.45 and 0.36 for horizontal and vertical accelerations, respectively (Hajehasani, 2001). By multiplying suitable coefficient, the Gilroy earthquake record altitude reached the maximum predictable earthquake altitude |
| Fig. 22: | Reduction of unbalanced force for primary equilibrium |
| Fig. 23: | Displacement along vertical axis of Z at section of 60.75 in terms of cm |
| Fig. 24: | Displacement along vertical axis of Z at section of 34.75 in terms of cm |
| • | Elimination of frequencies higher than dominant frequency of record: Optimum mesh size changes on the basis of frequency of input wave, thus, some record components which produce frequencies higher than dominant frequency of record, must be eliminated. Because such high frequencies require very small size of meshes in the model. If size of model mesh is identified on the basis of above high frequencies, the cpu time for solving this problem increases considerably. If we ignore these high frequencies, the problem of wave distortion arises, thus, the elimination of these data is inevitable |
| Fig. 25: | Displacement along vertical axis of Z at section of 85.75 in terms of cm |
Defining boundary condition for model for dynamic loads: Wave reflection and diffraction of models boundary is very important in dynamic analysis of underground structure. In fact the mediums around underground structures are infinitive but computer models are finite. In this investigation the boundary conditions of model was considered infinitive boundary which have been shown in Fig. 26. This assumption helps the computer model of Masjed Soleiman underground structure correspond better with real states.
Identification of damping conditions of the model: In this modeling, Rayleigh damping has been employed. In implementing Rayleigh damping, two parameters-optimum frequency of model and the rate of critical damping of environment should be identified. To obtain optimum frequency, which is a mixture of natural frequency of earth and dominant frequency of input wave, natural frequency of earth should be identified at first. Then the input wave frequency should be taken into optimum frequency obtained. The computer model was analyzed dynamically under the effect of its self-weight force and the displacement history was calculated along the vertical axis of one of its points to identify the natural frequency of model. The history mentioned is a cosine function whose frequency is, in fact, the natural frequency of the model. For the given model, the natural frequency of model is 6 Hz. The frequency of input wave is 2.5 Hz. Thus, optimum frequency of model may be assumed to be 4 Hz.
| Fig. 26: | Infinite boundary condition which has been used in modeling |
Critical damping was assumed 4%, considering the rate of this damping ranging from 2-5% for rock environments (Kirzhner and Rosenhouse, 2000).
Applying dynamic load and analysis of results: After identifying parameters regarding dynamic analysis, the records provided in the earlier part may be applied to the model and the results can be studied. The relevant records were applied to the structure in two forms:
| • | Record concerning to earthquake’s horizontal acceleration applied to the floor of structure as shear wave |
| • | Record concerning to earthquake’s vertical acceleration applied to the floor of structure as pressure wave |
After applying earthquake’s records to the model, stresses and strains which are due to dynamic loading can be numerically determined and can identify the Masjed Soleiman cavern stability:
The strain investigation show that the maximum shear strain have been generated in power house cavern roof, near the transformer cavern (Fig. 27), but the amount of them is in stable range because the amount of them is lower than critical shear strain (Table 3) that obtain from Sakurai equation (Sakurai et al., 1994) expression (3,4). Therefore the deformations that have been generated, of maximum design earthquake loads, will be in stable range:
log εc = -0.25 logE-1.59 |
(3) |
γc = (1+v)εc |
(4) |
| Fig. 27: | Maximum shear strain after seismic load |
| Table 3: | Critical shear strain (10) |
| γc | = | Critical shear strain |
| εc | = | Critical strain |
| E | = | Modulus of elasticity [kgf cm-2] |
| υ | = | Poisson’s ratio |
After exercising the seismic loads to model, the generated stresses in support system obtain lower than support system strength, therefore can expect that the Masjed Soleiman cavern support system will resist for maximum design earthquake loads.
This two reasons guide that we have concluded Masjed Soleiman cavern is stabile for seismic loads.
CONCLUSIONS
The results of this investigation can be expressed as follow:
| • | In numerical modeling of underground spaces used to study the effect of dynamic loads on the stability of underground space, the results obtained through decreasing (abasing) the elements dimensions of numerical model are more accurate. However, decreasing the dimensions more than an optimum limit has no effect on better results and it just results in an increased CPU time |
| • | Size of input waves latitude has no effect on truth of wave propagation and optimum length of elements models as well |
| • | The most prominent effective factor in the optimum size of mesh is input wave frequency in addition to physical and mechanical properties of materials. Thus with an increase of input wave frequency, the elements length must decrease. Therefore, by filtering the record applied, the data producing frequency more than the dominant wave frequency should be eliminated from time histories |
| • | The fact that underground structures are resistant against seismic waves is undeniable. But concerning the complicated underground spaces like Masjed Soleiman cavern in which some spaces are located near each other it should be considered more seriously. The results of this research showed that this cavern is stabile for maximum design earthquake but especial dangers such as faults deformations, liquefactions problems and wedges falling must be investigated |