**A simple model for stress-induced anisotropic softening of weak **

**sandstones **

M.C. Weng1_{ L.S. Tsai}2_{ Y.M. Hsieh}3_{ F.S. Jeng}4, *_{ }
1

Assistant professor, Department of Civil and Environmental Engineering, National University of Kaohsiung, Kaohsiung, Taiwan.

2

Engineer, China Engineering Consultants, Inc., Taipei, Taiwan.

3

Assistant professor, Department of Construction Engineering, National Taiwan University of Science and Technology, Taipei, Taiwan.

4

Professor, Department of Civil Engineering, National Taiwan University, Taipei, Taiwan.

*

Corresponding author, fsjeng@ntu.edu.tw; Tel/Fax:+886 2 2364 5734

**Abstract **

Weak sandstones possess deformational behaviors different from hard rocks; these phenomena, including shear dilation and softening of deformational moduli, are much more significant. It was found that: (1) Under hydrostatic loading, the bulk modulus increases as confining pressure arises; and (2) Under shear loading, the weak sandstone may transform from its original isotropy to a stress-induced anisotropy material and the deformational modulus can accordingly be softened as well. These phenomena contribute to the increase of crown settlements during tunnel excavations and accounts for several cases of tunnel squeezing.

Therefore, a model which is capable of simulating major deformational characteristics of weak sandstones is essentially needed for engineering purposes. A simple yet innovative constitutive model is accordingly proposed. This proposed model is characterized during the simulation as having: (1) non-linear volumetric deformation under hydrostatic loading, (2) significant shear dilation prior to the failure state, (3) isotropic stiffening of deformational moduli under hydrostatic loading and (4) anisotropic softening of deformational moduli under shearing condition.

The proposed model was formulated based on the linear elastic model, and it accounts
*for the variations of moduli E and G* through different loading conditions. It was found
**Manuscript**

that the proposed model closely simulates the actual deformational characteristics of weak sandstones. In addition, the proposed model only needs six material parameters and all these parameters can be easily obtained from experiments.

The proposed model was then incorporated into a finite element program and was used to analyze a squeezing tunnel case. Overall, this model can describe the deformation behavior for weak sandstones, especially on the significant shear dilation prior to the failure state. As a result, the proposed model shows the merits in being simple in form, yet being versatile in its applicability.

**Keywords: weak sandstones, constitutive model, shear dilation, deformational modulus, **
softening, stiffening.

**1. Introduction **

Engineering difficulties are often encountered during tunnel constructions in weak rocks, which possess undesirable engineering characteristics such as high porosity, poor cementation, and low strength and stiffness. In western part of Taiwan, weak sedimentary rocks (including sandstone, shale and mudstone), undergone through juvenile rock forming process, are often encountered during tunnel constructions, and these materials exhibit relatively low shear strength and stiffness. In several unsuccessful tunnel constructed in weak rock strata, ground squeezing and tunnel face instability have been reported [1]. In one particular case where a typical weak rock in Taiwan, known as Mushan sandstone, was encountered and several tunnel sections showed a crown settlement ranging from 14 cm to 30 cm, which is much larger than ordinarily a few centimeters.

In order to identify the deformational behavior of weak sandstones, several research
projects were conducted by the authors. Jeng et al. [2] categorized mechanical properties of
many types of sandstones using the petrographic features and categorized the analyzed
*sandstones into Type A and Type B sandstones. The Type B sandstone, compared to Type A, *
showed lower stiffness and larger reduction in both stiffness and strength during shearing.
*These characteristics suggested that the Type B sandstone could be a problematic rock type, *
*and tunneling in the Type B sandstone is prone to tunnel squeezing. *

*For the conventional triaxial compression (CTC) tests, both the confining pressure and *
shear stresses increase; therefore, it is difficult to distinguish the effects of shear and
volumetric stresses on the volumetric deformation. As such, triaxial tests with pure shear
stress path (Fig. 1) were performed on weak sandstones to study the effects of volumetric and
shear stresses on the volumetric deformation [3, 4]. Triaxial tests with pure shear stress path
*(PS) make it possible to separate the effects of shear stresses and volumetric stresses on *

deformations, and have been adopted in recent years for studying rock mechanics in the three-dimensional stress space. Based on the experiment results [3, 4], the deformation characteristics of weak sandstones include: (1) non-linear volumetric deformation during hydrostatic compression and the stiffness increase with the hydrostatic stress (Fig. 2); (2) significant shear dilation and distortion (Figs. 3 & 4) during shearing; and (3) substantial plastic deformation occurs prior to the failure state during shearing.

From micro-mechanical point of view, the plastic behavior is related to complex micro- and meso-crack coupled growth, and accumulation of these cracks and further application of loading will lead to stress-induced anisotropy and softening of rock [5-12]. Sayers et al. [7, 8] showed that initially isotropic rocks would exhibit anisotropic elastic behaviors under deviatoric shearing, and that the degree of anisotropy would be influenced by the micro-crack opening and their orientation. They have also showed the direction of crack propagation is often parallel to the direction of the maximum compressive stress applied. Lockner & Beeler [11] also pointed out that when the specimen was subjected to high differential stresses, new micro-crack-induced damage would be initialized, which in turn induced significant anisotropy. They have also found that increasing confining pressure not only increases the deformation moduli but also tends to suppress anisotropy. Olsson [12] indicated that changes in the incremental shear moduli reflect the development of deformation-induced, anisotropic damage. He measured the incremental shear moduli in two directions, and found that in the initial loading range the moduli were equal, but then the moduli decreased in different rates with increasing loadings.

In addition to understanding deformation behaviors, including shear dilation and stress-induced anisotropy of weak rocks, the assessment of the stress condition and displacements associated with these deformation characteristics is important in practical civil engineering or in rock engineering. Conventionally, rocks are often considered as

linear-elastic, followed by plastic behaviors before and after the stress state reaches the yielding condition of the material, respectively. This model probably works well for hard rocks, in which the shear dilation or the stress-induced anisotropy are not significant, but performs unsatisfactorily for soft rocks in the following regards:

(1) As shear dilation and stress-induced anisotropy exist, they need to be incorporated into a suitable constitutive model, especially for tunnel constructions in the squeezing ground, in which the volumetric deformation is significant; and

(2) The softening of the stiffness modulus occurs even in the so-called elastic stage, and this nonlinear elasticity should be considered in establishing a constitutive model for soft rocks.

Conventional elasto-plastic models may require quite a lot of material parameters on one hand, or may still have limited capability of describing the above-mentioned deformational characteristics of soft rocks on the other hand. Therefore, a simple yet effective constitutive model, which is characterized by only a few material parameters for soft rocks, is needed for engineering design and analysis. Accordingly, an innovative constitutive model accounting for stress induced anisotropy and softening of weak rocks is proposed in this study.

The proposed model was designated to be capable of showing the following characteristics:

1. The model gives reasonable predictions on apparent shear dilation, including the stress-induced, anisotropic softening;

2. The relatively simple model can be readily applicable to engineering practice by requiring a simple and straightforward parameter determination procedure;

3. Relatively fewer material parameters, compared to elasto-plastic models, are needed in the proposed model; and

4. The major deformational behaviors of weak rocks can be well described by the proposed model.

**2. Model **

**formulation **

**2.1 Elastic modeling for stress induced anisotropic softening **

For isotropic linear elastic material, the matrix form of the incremental stress-strain relations in the principal stress space can be expressed as:

⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ∆ ∆ ∆ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − − − − − − − − − − = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ∆ ∆ ∆ 3 2 1 3 2 1 1 2 2 2 2 2 2 1 2 2 2 2 2 2 1 1

### σ

### σ

### σ

### ε

### ε

### ε

*G*

*G*

*E*

*G*

*G*

*E*

*G*

*G*

*E*

*G*

*G*

*E*

*G*

*G*

*E*

*G*

*G*

*E*

*E*(1)

*where E and G are tangent Young’s modulus and shear modulus; * ∆ , σ_{1} ∆ and σ_{2} ∆ are σ_{3}
the stress increments along the three directions of principal stresses; and ∆ , ε_{1} ∆ and ε_{2} ∆ ε_{3}
are the corresponding strain increments.

The laboratorial experiments, as shown in Fig. 2*, found that the bulk modulus K can be *
related to the hydrostatic stress as:

*b*
*I*
*a*
*K* = +
3
1
(2)

*where I*1* is the first stress invariant; a and b are the material parameters. Parameter b is the *
*bulk modulus of the material at zero confinement. Parameter a determines the rate of *
increase of bulk modulus with increasing confining pressure.

A constitutive model for weak sandstones can be developed to describe the softening and
the stress-induced anisotropy by modifying Eqn. 1. The concept of the proposed model is
described as follows: the weak sandstone subjected to hydrostatic stress condition can be
viewed as an isotropic material with deformation modulus that increases or stiffens when
*confining pressure increases. As such, the material parameters a and b in *Eqn. 2 can be
related to the initial shear modulus (i.e. maximum shear modulus) based on the isotropic
elasticity theory as:

1
max
3 (1 2 ) 3 (1 2 )
( )
2 (1 ) 2 (1 ) 3
*K* *I*
*G*

### ν

### ν

*a*

*b*

### ν

### ν

− − = = + + + (3)*where ν is the Poisson’s ratio in *Eqn. 3.

When the material is subjected to shear loadings, the shearing may cause both softening
of deformation moduli and stress-induced anisotropy. Like soil materials, the stress-strain
behavior of weak sandstones is highly nonlinear at all phases of shear loading with a
*maximum shear modulus * *G*_{max} at initial, small strain. Therefore, the modulus degradation
(softening) can be expressed as a function (Eqn. 4) of mobilized stress level, or *G G*/ _{max}
versus shear stress. Similar relationships have also been proposed for other geomaterials.

[13-15].
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
= 2
,
2
2
max 1 ( )
*f*
*J*
*J*
*G*
*G* (4)

where shear strength *J*_{2}_{,}* _{f}* can be calculated by assuming a linear yield criteria which has

the form as: *J*_{2}_{,}* _{f}* =α

*I*

_{1}+

*k*; and

*J is the second deviatoric stress invariant.*

_{2}

Based on observations from experiments [2, 11, 12], it is was found that the application of shear stress will induce anisotropy and two assumptions can accordingly be made for the proposed constitutive model:

(1) The direction of “anisotropic” softening is assumed to be coincident with the axis of the major principal stress based on experimental observations; and

(2) The shear modulus degradations parallel to the major principal stress axis can be described
by Eqn. 4. Furthermore, the shear modulus degradation normal to the major principal
stress direction can be described using a similar relationship as Eqn. 4 by introducing an
*anisotropic parameter e ranging from 0 to 1: *

⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
=
′ 2
,
2
2
max 1 ( )
*f*
*J*
*J*
*e*
*G*
*G* (5)

Therefore, for materials like weak sandstones, the matrix of the stress-strain relation in the principal stress coordinate can be expressedas:

⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ ∆ ∆ ∆ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − − − − − − − − − − − = ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ ∆ ∆ ∆ 3 2 1 3 2 1 1 ' 2 ' 2 2 2 2 ' ' 2 1 2 2 2 2 2 2 1 1

### σ

### σ

### σ

### ε

### ε

### ε

*G*

*G*

*E*

*G*

*G*

*E*

*G*

*G*

*E*

*G*

*G*

*E*

*G*

*G*

*E*

*G*

*G*

*E*

*E*(6)

If *G G′*= , the matrix is identical to the matrix of an isotropic material and the material
is still an isotropic material. As will be shown later, *G′* can differ from *G* under certain

loading conditions and the matrix has a form different from the matrix of isotropic materials. Under such circumstance, the material is referred as a stress-induced anisotropic material.

Remarkably, if *G*′ ≠*G*, the material is anisotropic and the increments (∆ , σ_{1} ∆ , σ_{2} ∆ , σ_{3}
1

ε

∆ , ∆ and ε_{2} ∆ ) in ε_{3} Eqns. 1 and 7, should be associated with the major, intermediate and

minor stress or strain increments, respectively.

*The deformation modulus E can be expressed in terms of * *G* and ν , based on Eqn. 7,
as:
2
2
max
2,
2 (1 )
2 (1 ) 1 ( )
*f*
*E* *G*
*J*
*G*
*J*

### ν

### ν

= + ⎛ ⎞ ⎜ ⎟ = + − ⎜ ⎟ ⎝ ⎠ (7)From Eqns. 4, 5 and 7*, the variations of moduli E , * *G* and *G′* in the stress space are
plotted in Fig 5, which clearly depicts how Young’s modulus and shear modulus varies with
different combinations of hydrostatic stress and shear stress.

Based on Eqn. 6 and the elasticity theory, the volumetric strain increment can be obtained as Eqn. 8, and the volumetric strain can be further separated into two components: volumetric deformation due to hydrostatic stress

## ( )

∆ε

_{v}*and volumetric deformation due to shear stress*

_{p}## (

∆ε

_{v}## )

*as the following expressions:*

_{s}## ( ) ( )

1## (

## )

2 3 2 ' 3 ' ' 2 9 3 1*J*

*G*

*G*

*G*

*G*

*I*

*G*

*E*

*G*

*E*

*E*

*s*

*v*

*p*

*v*

*v*∆ ⋅ − + ∆ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛

_{−}

_{−}= ∆ + ∆ = ∆

### ε

### ε

### ε

(8)## ( )

_{1}' 2 9 3 1

*I*

*G*

*E*

*G*

*E*

*E*

*p*

*v*⎟∆ ⎠ ⎞ ⎜ ⎝ ⎛

_{−}

_{−}= ∆ε (9)

## ( ) (

## )

2 3 2 ' 3 '

_{J}*G*

*G*

*G*

*G*

*s*

*v*

_{⋅}∆ − = ∆ε (10)

Eqn. 10 shows two interesting features of the proposed model: 1) Shearing induces
volumetric deformations, which cannot be described using isotropic linear elastic models; and
2) Shearing can produce both compressive and dilative volumetric strains. If *G* is greater
than *G′*, the volumetric strain increment

## ( )

∆### ε

_{v}*will be compressive; if*

_{s}*G*is smaller than

*G′*, the

### ( )

∆ε

_{v}*will be dilative. As such, this proposed model has the ability of properly describing the shear-induced deformation observed in lab experiments.*

_{s}**2.2 Determination of material parameters **

*The elastic moduli E , * *G and * *G′* of the proposed elastic model can be deduced from
*six material parameters a, b, e, ν, *α* and k. Parameters a and b determine the volumetric *
behavior of the material under isotropic loading, as shown in Eqn. 2*. Parameter e dictates *
the stress-induced anisotropy, or different degrees of shear modulus degradations in different
principal directions, as in Eqn. 4 and 5*. Parameter ν is the Poisson’s ratio used in *Eqn. 3 and
7. Parameters α* and k are parameters related to the shear strength of the material, used to *
define softening behavior, as in Eqn. 4 and 5.

This section demonstrates how these six parameters can be determined during laboratory experiments by using an example of the Mushan sandstone, a weak rock commonly seen in mountain areas in Taiwan.

The Mushan sandstone, which has caused squeezing ground hazards for tunnels in Taiwan, has a porosity of 14.1%, dry density of 2.3 g/cm3, and saturated water content of 5.8%. The average uniaxial compressive strength is 37.1 MPa in dry condition and 28.9 MPa in saturated condition. Based on the petrographic analyses, the percentages of grains,

matrix, and voids are 67.5%, 18.5% and 14.1%. The average particle diameter is 0.72 mm. Mineralogically, Mushan sandstone consists of 90.7% of quartz, 9.0% of rock fragments, and thus is classified as lithic greywacke [2, 3].

To determine the six material parameters required for the proposed constitutive model, two laboratory experiments are needed: 1) uniaxial compression tests, and 2) triaxial compression tests with pure shear stress path.

*Poisson’s ratio ν can be obtained from uniaxial compression test at the initial loading *
stage. From hydrostatic loading stage of triaxial compression tests, as shown in Fig. 6,
*material parameters a and b can be determined by curve fitting *Eqn. 2 with experimental
curves. It can be seen from Fig. 6 that the bulk modulus almost linearly varies with the
confining pressure, and thus the proposed linear relationship (Eqn. 2) is applicable and can be
accordingly determined by fitting the experimental curve shown in Fig. 6.

During shear loading stage of triaxial compression tests, parameters α*, k, and e can be *
determined. In this proposal model, a yield criterion similar to the Drucker-Prager model is
assumed and used to determine the shear strength of the material. These three parameters in
turn will affect the shear modulus degradation during softening. The parameters α *and k are *
the friction angle and the cohesion intercept of the yield criterion, as shown in Fig. 7. The
*parameter e, which ranges from 0 to 1, reflects the degree of shear-induced anisotropy. A *
*decreasing e increases the difference between shear moduli * *G and * *G′* as well as it results
in more shear dilation, which is illustrated in Fig. 8*. The parameter e can be back-calculated *
from the volumetric deformation under shear loading.

**3. Results of simulations using the proposed model **

The performance of the proposed model is demonstrated by comparing the laboratory results with the model simulations. The versatility of the proposed constitutive model is

then tested by using the same set of fitted parameters to predict material responses under different testing conditions or stress paths and these predictions are compared with laboratory experiment data.

**3.1 Simulations of experimental results **

*One pure shear (PS) experiment conducted using the Mushan sandstone in the laboratory *
is simulated by the proposed model with fitted parameters for this experiment listed in Table
1*, and the measured and simulated deformation stress-strain curves obtained from PS *
experiments are shown in Fig. 9. At the hydrostatic loading stage of the experiment, it is
seen, as in Fig. 9a, that the measured nonlinear volumetric behavior can be well simulated by
the propose model. As the experiment continues at the shear loading stage, as shown in Fig.
9b, the material shows significant shear dilation, which is also well simulated by the model
simulation. As a result, the proposed model is capable of modeling deformational behavior
of weak sandstones under both hydrostatic loading and shear loading conditions.

*Deducing from the experimental results, the variations of moduli E , * *G and * *G′* for the
two loading stages are plotted in Fig 9c. At first, all moduli increase linearly, as described
by Eqns. 2 and 3, under hydrostatic loading. Afterwards, as shear stress is applied, all
moduli degrade exponentially according to Eqns. 4 and 5, and *G* starts to differ from *G′*
upon increasing shearing, as shown in Fig. 9c. Therefore, the studied sandstone deforms
isotropically and anisotropically upon hydrostatic and pure shear, respectively. Based on

Eqn. 7, the proposed model is capable of simulating the stress-induced anisotropy of a material, which is initially isotropic and the softening of the deformation can be simulated as well. Overall, the proposed model is capable of simulating the deformation behavior of weak sandstones induced by both the hydrostatic and pure shearing loadings.

**3.2 Versatility of the proposed model for different stress paths **

In order to evaluate the predictive capability of the proposed model, the same set of material parameters listed in Table 1 is used for the proposed model to predict material responses under various confining stresses and shear stress paths, and these predictions are then compared to laboratory measures.

Figure 10* shows the PS results of both model predictions and experimental results under three *
*different hydrostatic stresses, p = 20, 40 and 60 MPa. It can be seen in *Fig. 10 that the proposed
model gives reasonable predictions, which compares well to the measured deformation in the
aspects of volumetric deformation (Fig. 10a) and shear deformation (Fig. 10b). It should be
*emphasized that the parameters are obtained by fitting the experimental results tested under p = 20 *
*MPa, thus the simulations for p = 40 and 60 MPa are actual predictions rather than results obtained *
from curve fittings.

*To further evaluate the versatility of the proposed model, different stress paths, the CTC *
tests, are simulated by the proposed model using parameters in Table 1, which was obtained from
*the PS experiment under p = 20 MPa. As a result, *Figure 11 shows the model prediction compares
*fairly well with the results of CTC test with minor discrepancy. This discrepancy possibly *
originates in the differences in deformation moduli of different specimens; however, as the major
concern of the proposed model is to simulate shear dilation, the observed discrepancy is acceptable.

The predictive capability of the proposed model is also evaluated in a triaxial extension mode.
*The Reduced Triaxial Extension test (RTE), which induces shearing by reducing the axial stress *
after hydrostatic loading stage, was performed in the laboratory and predicted by the proposed
model. Again, a fairly good agreement between the two is seen in Fig. 12, suggesting that the
proposed can predict material responses in a hydrostatic pressure unloading conditions.

Overall, the proposed model can reasonably predict or predicted the deformational behavior of weak sandstones over a wide range of confining pressures and stress paths.

**4. Application to tunnel excavations **

As the proposed model for weak sandstones is capable of modeling material responses realistically with considerations of anisotropic softening, the impact of such behavior to the deformation induced by tunnel excavations can be accordingly analyzed by numerical analyses. The proposed constitutive model was incorporated into a finite element code ABAQUS and is then used to predict the wall deformation of the tunnel.

A well-known tunneling project in Taiwan is selected as a case study. This tunnel, with a height of 16m, a width of 12.4m, and a thickness of overburden ranging from 20m to 170m, was constructed in the range of the Western foothill of Taiwan, in which weak sedimentary rocks prevail. The rock mass along the tunnel was mostly the aforementioned Mushan sandstone. Drill and blast method was used to excavate this tunnel, and the excavation was done in two stages: top heading excavation and invert excavation.

Two constitutive models, linear-elastic-perfectly-plastic with Drucker-Prager yield
criterion and the proposed model, are used to analyze the side wall deformation of the studied
tunnel, and the material parameters and properties are listed in Table 2 based on the studied
rock mass. A section of the tunnel, known to have suffered from ground squeezing, is
shown in Fig. 13 based on the geometry described in Jeng et al. [1]. The tunnel was
analyzed in three steps: geostatic condition with *K _{o}* = 1.0, top heading excavation and invert
excavation. The results of the analysis are presented as the following.

Figure 14 illustrates the variations of the shear modulus *G* in the cross section
throughout the analysis. It could be observed that significant degradations of the shear
moduli *G* at crown, sidewall and invert take place during top heading excavation (stage 2).
These moduli remain very low at the subsequent stage. It is found that softening initiates at
the upper sidewalls first, then at the invert and at last at the crown. Finally, deformational

modulus distributions, including Young’s modulus and shear modulus, are shown in Fig. 15, which reveals that significant softening occurred all around the tunnel, and the tunnel invert has the most degradation.

The stress paths at the crown, the invert and sidewalls are shown Fig. 16. Similar stress
*paths close to RTE path are seen at the crown and at the invert; namely, the material is *
increasingly sheared while the hydrostatic pressure decreases. For sidewalls, however, the
*stress path is similar to CTC path. At last, the stresses at the crown, the invert and the upper *
sidewall of inner section reach the failure envelope.

The inward displacements at selected locations on tunnel walls calculated from finite element calculations based on the two material models used are summarized in Table 3. It shows that the calculations of inward displacements, based on the proposed model, are larger all over the tunnel section, including the crown, sidewalls and the invert, compared to those predicted by elasto-plastic model. The actual crown settlement ranges from 14 ~ 30 cm. As such, without incorporating the plastic model into the proposed model, the proposed model has predicted a crown settlement closer to the reality than the one predicted by the conventional elasto-plastic model.

What accounts for the significant amount of inward ground displacement movements obtained by the proposed model can be further explored by examining the volumetric strain distributions of the two models. Figures 17a and 17b show the distribution of volumetric deformation around the tunnel obtained by the elasto-plastic model and the proposed model, respectively. As shown in Fig. 17a, a dilation zone beneath the invert (the shaded areas) was developed based on the elasto-plastic model. Consequently, the elasto-plastic model predicts greater displacements at the invert than at other locations. However, the dilation zone calculated by the proposed model is much larger than that calculated by the

elasto-plastic model, especially above the crown and beneath the invert, as shown in Fig. 17b. As a result, more dilation pushes the crown more inward, as revealed by the proposed model. The shear dilation prior to yielding accounts for the larger extent of the dilation zone, which cannot be simulated by the conventional elasto-plastic models since they would not allow dilation in the elastic range.

Remarkably, the crown settlement predicted by the proposed constitutive model is about 10.2 cm, which is closer to, but still somewhat underestimated, the measured crown settlement of 14 ~ 30 cm (Jeng et al., [1]). Therefore, other factors such as scale effect and creep may account for the discrepancy and require further study.

**5. Conclusion **

An innovative constitutive model was proposed for weak sandstones to describe the deformational characteristics exhibited by weak sandstones and other geomaterials. This proposed constitutive model is characterized by the following features: (1) the non-linear volumetric deformation under hydrostatic loading; (2) the significant shear dilation prior to the failure state, (3) modulus stiffening under hydrostatic loading and anisotropic softening under shearing condition; and (4) being versatile in various stress paths.

The stress-strain relationship of the model, which accounts for the variations of modulus

*E and G through different loading conditions, was formulated based on the framework of *

linear elasticity. Under hydrostatic stress condition, the weak sandstone is viewed as an isotropic material, and its stiffness increases as confining pressure arises. However, when the weak sandstone is subjected to shearing, the material not only becomes anisotropic but also softens. Such deformational characteristics have been observed through numerous experiments in the literature. As such, the modulus degradation is expressed as a decaying binomial function of the shear stress, and the shear modulus for the plane parallel to the

maximum principal stress varies differently from those at the other two directions, accounting for shear-induced anisotropy. Finally, the proposed model was tested through comparisons between experimental data and model predictions and found to be versatile for various stress paths.

Compared to other sophisticated elasto-plastic models, the proposed model needs only six material parameters, all of which have physical meaning and can be obtained easily. The model has been incorporated into the finite element program ABAQUS and used to analyze a squeezing tunnel case in Taiwan. Larger extent of dilation zones accounting for the greater crown settlement, obtained using the proposed model and this prediction is closer to the actual crown settlement than the prediction obtained by using the conventional elasto-plastic model. Overall, the proposed model can reasonably describe the deformation behavior, with significant shear dilation prior to the failure state for weak sandstones.

**6. Acknowledgements **

The research is supported by the National Science Council of Taiwan, Grant no. NSC-95-2221-E-002-263.

**7. References **

[1] Jeng, F.S. Lin M.L. and Huang T.H. Study of the Geological Barriers of the Tunnels in Northern Taiwan (in Chinese). Taipei: Ministry of Transportation, Research Report, 1996. [2] Jeng, F.S. Weng, M.C. Lin, M. L. and Huang, T.H., Influence of petrographic parameters on

geotechnical properties of Tertiary sandstones from Taiwan. Engineering Geology 2004, 73:71–91.

[3] Weng M.C. Mechanical characteristics and the relations with microstructure factors of foothill sandstones. Taipei: National Taiwan University Civil Engineering, Doctoral Dissertation, 2002. [4] Weng, M.C. Jeng, F.S. Huang, T.H. and Lin, M.L., Characterizing the Deformation Behavior

[5] Wu, B. King, M.S. and Hudson, J.A., Stress-induced ultrasonic wave velocity anisotropy in a sandstone. Int J Rock Mech Min Sci Geomech Abstr, 1991, 28:101-107.

[6] David, C. Menendez, B. and Darot, M., Influence of stress-induced and thermal cracking on physical properties and microstructure of La Peyratte granite. Int J Rock Mech Min Sci, 1999, 36:433-448.

[7] Sayers, C.M., Stress-induced ultrasonic wave velocity anisotropy in fractured rock. Ultrasonics, 1988, 26:311-317.

[8] Sayers, C.M. Van Munster, J.G. King, M.S., Stress-induced ultrasonic anisotropy in Berea sandstone. Int J Rock Mech Min Sci Geomech Abstr, 1990, 27:429-436.

[9] Sayers, C.M., Stress-dependent elastic anisotropy of sandstones. Geophysical Prospecting, 2002, 50:85-95.

[10] Ita, S.L. Cook, N.G.W. Myer, L.R. Nihei, K.T., Effects of stress anisotropy on the static and dynamic properties of Berea sandstone. Int J Rock Mech Min Sci Geomech Abstr, 1993, 30:785-788.

[11] Lockner, D.A. and Beeler N.M., Stress-induced anisotropic poroelasticity response in

sandstone. In: Electronic Proc. 16th ASCE Engin. Mech. Conf., Seattle, WA, 2003. p.16-18. [12] Olsson W., Development of anisotropy in the incremental shear moduli for rock undergoing

inelastic deformation. Mech Mater, 1995; 21:231–42.

[13] Tatsuoka, F. and Shibuya, S., Deformation characteristics of soils and rocks from field and laboratory tests. Rept. of the Inst. of Industrial Science, 1992, 37:136.

[14] LoPresti, D.C.F. Pallara, O. Lancellotta, R. Armandi, M. and Maniscalco, R., Monotonic and cyclic loading behavior of two sands at small strains. ASTM Geotechnical Testing Journal, 1993, 16:409-424.

[15] Tatsuoka, F. Jardine, R.J. LoPresti, D.C.F. DiBenedetto, H. and Kodaka, T., Theme Lecture: Characterizing the pre-failure deformation properties of geomaterials. In: 14th International Conf. on Soil Mechanics & Foundation Engineering, Vol. 4, Hamburg, 1997, p.35.

**8. Nomenclature **

*a * Material parameter for proposed model

*b * Material parameter for proposed model

*CTC Conventional triaxial compression test *
*E * Tangent Young’s modulus

*e * Material parameter for proposed model

*G* Tangent shear modulus

*G′* Tangent shear modulus for the plane parallel to the

maximum principal stress

*o*

*G * Initial tangent shear modulus
max

*G* * Maximum tangent shear modulus *
1

*I* First stress invariant

2

*J* Second deviatoric stress invariant, _{J}* _{s}_{ij}_{s}_{ji}*
2
1

2=
*J2f* Second deviatoric stress invariant at failure

*K * Tangent bulk modulus

*o*

*K * The coefficient of lateral stress at rest

*k * Interception of Drucker-Prager failure criteria

*p * Hydrostatic stress *p* *I* σ* _{kk}*
3
1
=
=

_{1}3 1 (MPa)

*PS * Pure shear test

*RTE Reduced triaxial extension test *
*ij*

*s * Second deviatoric stress tensor

*T * *Interception of failure envelope with I1* axis

*UCS Uniaxial compressive strength (MPa) *

α Slope of Drucker-Prager failure criteria

*i*

ε

*∆ Strain increment in the principal stress i direction. *

*v*

ε

∆ Volume strain increment

## ( )

∆ε*v*

*p*Volumetric strain increment due to hydrostatic stress

## (

∆ε*v*

## )

*Volumetric strain increment due to shear stress*

_{s}γ _{Shear strain }

*ν * Poission’s ratio

*ij*

**List of figure captions **

*Figure 1 Schematic illustration of stress paths. PS = Pure Shear test, CTC = Conventional *
*Triaxial Compression test, RTE = Reduced Triaxial Extension tests. *

Figure 2 Typical deformation curves of volumetric strain induced by volumetric stress (Path OA in Fig. 1) obtained from weak sandstones .

Figure 3 Typical deformation curves of shear strain induced by shear stress (Paths AC in Fig. 1).

Figure 4 Typical deformation curves of volumetric strain induced by shear stress (Path AC in

Fig. 1).

*Figure 5 Schematic illustration of the variation of the deformation modulus E , * *G* and *G′*
in the hydrostatic stress and shear stress space.

*Figure 6 Variation of tangent bulk modulus K during hydrostatic compression state of *
loading. Three test results are shown in the figure. By fitting the experiment results,
the material parameters *a* and *b* are accordingly determined, as listed in Table 1.
Figure 7 Failure envelope of the studied weak sandstone.

*Figure 8 Influences of the material parameter e, which ranges from 0 to 1, on the simulated *
dilation deformation based on the proposed model.

*Figure 9 Comparison of the deformation predicted by the proposed model and by the PS *
experiments. (a) Variation of volumetric strain during hydrostatic loading; (b) Shear
dilation induced by pure shear loading; (c) Variation of deformational moduli deduced
from experiments

Figure 10 Simulation of shear and volumetric strains induced by shear stress under different
*hydrostatic pressures (PS tests). From *Figs. 9 to 11, the material parameters were
*obtained from PS test under p = 40 MPa. *

*Figure 12 Simulation of volumetric strains induced by shear stress for RTE tests. *

Figure 13 Schematic illustration of the geometry of the studied tunnel.

*Figure 14 Variations of the shear modulus G around the cross section during analysis. The *
*shear modulus G at crown, sidewalls and invert initiate degradation at stage 2 (top *
heading).

*Figure 15 Distribution of E and G after tunnel excavation. (a) Young’s modulus E ; (b) *
shear modulus *G*.

Figure 16 Stress paths of the materials located at the crown, invert and sidewalls during tunnel excavation.

Figure 17 Volumetric deformation predicted by different models. Dilation zones are illustrated by shaded areas. (a) Elasto-plastic model; (b) Proposed model.

Table 1 Material coefficients of Mushan Sandstone for the proposed constitutive model

Coefficients Remark Value

*a * _{Material parameter } _{280 }
*b * _{Material parameter } _{621 MPa }
*ν * _{Poission’s ratio } _{0.21 }
*e * _{Material parameter } _{0.8 }

α Slope of failure criteria 0.36
*k * Interception of failure

criteria 9.44 MPa
**Table 01**

Table 2 Corresponding parameters for the analyzed tunnel

Type of model Parameter Remark Value
*a * _{Material parameter } _{20 }
*b * _{Material parameter } _{80 }
*ν * _{Poission’s ratio } _{0.21 }
*e * _{Material parameter } _{0.8 }

α Slope of failure criteria 0.16 Proposed model

*k * Interception of failure _{criteria } 82 kPa

Ε Young’s modulus 200 kPa
*ν * Poission’s ratio 0.21

α Slope of failure criteria 0.16
Elasto-Plastic
model
(Drucker–Prager
model) *
*k * Interception of failure
criteria 82 kPa
* Associated flow rule is assumed.

Table 3 Inward displacements predicted by two models

Displacemen t (cm) model

Crown Sidewall Invert

Proposed model 10.2 8.6 15.8
Elasto-plastic model 6.3 6.5 12.6
**Table 03**

Figure 1
**Figure 01**

**0**
**10**
**20**
**30**
**40**
**50**
**60**
**70**
**0** **2000** **4000** **6000** **8000** **10000** **12000**
**Volumetric strain **

### εεεε

**v****(10 -6 )**

**Hy**

**d**

**ro**

**st**

**at**

**ic**

**st**

**re**

**ss**

**p****(M**

**P**

**a)**p=60 MPa p=40 MPa p=20 MPa Figure 2

**Figure 02**

**0**
**10**
**20**
**30**
**40**
**50**
**60**
**70**
**80**
**0** **4000** **8000** **12000** **16000** **20000**
**Shear strain **

### γγγγ

**(10 -6)**

**She**

**ar**

**s**

**tr**

**es**

**s**

**(J**

**2**

**)**

**0.**

**5****(**

**M**

**P**

**a)**p=60 MPa p=40 MPa p=20 MPa Figure 3

**Figure 03**

**0**
**20**
**40**
**60**
**80**
**100**
**-6000** **-5000** **-4000** **-3000** **-2000** **-1000** **0** **1000** **2000**
**Volumetric strain **

### εεεε

**v****(10 -6 )**

**She**

**ar**

**s**

**tr**

**es**

**s**

**(J****2**

**)**

**0.**

**5****(M**

**P**

**a)**p=60 MPa p=40 MPa p=20 MPa

**Dilation**

**Compression**Figure 4

**Figure 04**

### (a) Variation of

*E*

### (b) Variation of

*G*

*and G′ *

Figure 5

**K = 280 I1/3 + 621**
**R2 = 0.97**
**0**
**2000**
**4000**
**6000**
**8000**
**10000**
**12000**
**14000**
**16000**
0000 10101010 20202020 30303030 40404040 50505050 60606060

**Hydrostatic stress p (MPa)**

**Hydrostatic stress p (MPa)**

**B**

**ulk**

** m**

**od**

**ulu**

**s **

**K**

**K**

** (M**

**P**

**a)**

Figure 6
**Figure 06**

R2 = 0.98
**0**
**20**
**40**
**60**
**80**
**100**
**120**
**0** **50** **100** **150** **200** **250** **300**
**I****1**** (MPa)**
** She**
**ar**
** s**
**tr**
**es**
**s (J****2****)****0.****5**** (M**
**P**
**a)**
44
.
9
36
.
0 1
,
2 = *I* +
*J* *f*
Figure 7
**Figure 07**

**0**
**10**
**20**
**30**
**40**
**50**
**60**
**70**
**80**
**-8000** **-6000** **-4000** **-2000** **0**
**Volumetric strain **

### εεεε

**v****(10-6)**

**S**

**h**

**ea**

**r st**

**re**

**ss**

**(J**

**2**

**)**

**0.**

**5**

**(M**

**P**

**a)**

**e=1**

**e=0.9**

**e=0.8**

**e=0.5**

**e=0**Figure 8

**Figure 08**

**0**
**5**
**10**
**15**
**20**
**25**
**0** **2000** **4000** **6000** **8000**
**Volumetric strain **εεεε**v**** (10-6)**
**H**
**ydrostatic st**
**res**
**s **
**p**** (M**
**P**
**a)**
Measured def.
Simulation

(a) Variation of volumetric strain during hydrostatic loading

**0**
**5**
**10**
**15**
**20**
**25**
**30**
**35**
**-6000** **-5000** **-4000** **-3000** **-2000** **-1000** **0** **1000**
**Volumetric strain **εεεε**v**** (10-6)**
**Shear str**
**ess **
**(J****2****)****0.****5**** (M**
**Pa**
**)**
Measured def.
Simulation

(b) Shear dilation induced by pure shear loading

**0**
**2**
**4**
**6**
**8**
**10**
**12**
**0** **10** **20** **30** **40** **50** **60**

**Axial stress (MPa)**

**M**
**odul**
**us (G**
**Pa**
**)** **E****G'****G**

Hydrostatic loading Pure shear loading

**O** **A** **C**

(c) Variation of deformational moduli deduced during experiments Figure 9

**0**
**10**
**20**
**30**
**40**
**50**
**60**
**70**
**80**
**-6000** **-5000** **-4000** **-3000** **-2000** **-1000** **0** **1000**
**Volumetric strain **

### εεεε

**v****(10 -6 )**

**S**

**h**

**ea**

**r st**

**re**

**ss**

**(J****2**

**)**

**0.**

**5****(**

**M**

**P**

**a)**p=60 MPa p=40 MPa p=20 MPa Simulation

(a) Shear dilation under various hydrostatic pressure conditions

**0**
**10**
**20**
**30**
**40**
**50**
**60**
**70**
**80**
**0** **4000** **8000** **12000** **16000** **20000**
**Shear strain **

### γγγγ

**(10 -6)**

**S**

**hea**

**r s**

**tres**

**s**

**(J****2**

**)**

**0.**

**5****(**

**M**

**P**

**a)**p=60 MPa p=40 MPa p=20 MPa Simulation

(b) Shear deformation under various hydrostatic pressure conditions

Figure 10
**Figure 10**

**0**
**10**
**20**
**30**
**40**
**50**
**60**
**70**
**80**
**-8000** **-6000** **-4000** **-2000** **0** **2000** **4000**
**Volumetric strain **

### εεεε

**v****(10**

**-6**

**)**

**S**

**hea**

**r s**

**tres**

**s**

**(J****2**

**)**

**0.**

**5**

**(M****P**

**a)**CTC, p=18 MPa Simulation Figure 11

**Figure 11**

**0**
**5**
**10**
**15**
**20**
**25**
**30**
**35**
**40**
**-6000** **-5000** **-4000** **-3000** **-2000** **-1000** **0**
**Volumetric strain **

### εεεε

**v****(10 -6 )**

**She**

**ar**

**s**

**tr**

**es**

**s (J**

**2**

**)**

**0.**

**5****(**

**M**

**P**

**a)**RTE, p=40 MPa Simulation Figure 12

**Figure 12**

Figure 13
**Figure 13**

**0**
**10**
**20**
**30**
**40**
**50**
**60**
**70**
**80**
**90**
**Excavation Process**
**Sh**
**ea**
**r mo**
**dul**
**us**
** G**** (M**
**P**
**a)** Crown_{Invert}
Upper sidewall
Lower sidewall

**Top heading excavation** **Invert excavation**

Figure 14
**Figure 14**

(a) Elastic modulus E

(b) Shear modulus G Figure 15

**0.0**
**0.2**
**0.4**
**0.6**
**0.8**
**1.0**
**1.2**
**1.4**
**1.6**
**1.8**
**0.0** **0.5** **1.0** **1.5** **2.0** **2.5** **3.0**

**Hydrostatic stress p (MPa)**

**She**
**ar**
** s**
**tr**
**es**
**s **
**(J****2****)****0.****5**** (M**
**P**
**a)**
Crown
Invert
Upper sidewall
Lower sidewall
**Failure envelope**
Figure 16
**Figure 16**

Figure 17
**Figure 17**