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 abovementioned 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 onedimension as well as
obtaining the necessary informations. For this purpose, by computer modeling
as well as numerical resolution of onedimensional 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 ONEDIMENSIONAL ELASTIC LONGITUDINAL WAVES PROPAGATION
To obtain further information concerning modeling of dynamic wave propagation,
there has been an attempt to model the propagation of onedimensional
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
onedimensional 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 longitudinalwave velocity.

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).
Where
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. 1013. 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. 1419.
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):
Where:
λ 
= 
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 
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 anchorbolt (Moshanir et
al., 1992). For modeling shotcrete and anchorbolts, 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. 2325.
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 parametersoptimum
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 selfweight 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 25% 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 logE1.59 
(3) 

Fig. 27: 
Maximum shear strain after seismic load 
Table 3: 
Critical shear strain (10) 

Where:
γ_{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 