# Constitutive Laws for Geomaterials

E. Papamichos

► To cite this version:

E. Papamichos. Constitutive Laws for Geomaterials. Oil & Gas Science and Technology - Revue d'IFP Energies nouvelles, 1999, 54 (6), pp.759-771. 10.2516/ogst:1999064. hal-02075850

HAL Id: hal-02075850

https://hal.science/hal-02075850

Submitted on 21 Mar 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constitutive Laws for Geomaterials

E. Papamichos¹

1 SINTEF Petroleum Research, N-7465 Trondheim - Norway

e-mail: euripides.papamichos@iker.sintef.no

**Résumé** — Lois constitutives pour les géomatériaux — Il existe un nombre important de lois constitutives permettant d'approcher le comportement des géomatériaux sous diverses contraintes, conditions aux limites et divers chemins de chargement, pressions de pore, températures, etc. L'élaboration de lois constitutives suit deux tendances concurrentes : la tendance à élaborer une loi unifiée décrivant le comportement du matériau sous les conditions les plus générales, et le besoin de lois justifié par une utilisation efficace dans la pratique. La dernière tendance inclut à la fois la possibilité d'une calibration par des essais sur des échantillons disponibles, et la convivialité de la loi. À des fins d'ingénierie, la seconde tendance domine, l'ingénieur doit évaluer le problème et sélectionner le modèle le plus approprié pour décrire le phénomène dominant observé. Dans cet article, un examen des diverses lois constitutives pour les géomatériaux, et en particulier pour leurs applications, est présenté. Divers exemples de problèmes concernant l'ingénierie et les géomatériaux appliqués au pétrole sont utilisés afin de démontrer que la loi doit être aussi simple que possible, mais pas trop non plus.

Mots-clés : loi constitutive, grès, craie, injection d'eau, endommagement.

**Abstract** — Constitutive Laws for Geomaterials — An abundance of constitutive laws exists to approximate the behavior of geomaterials under various stresses, boundary conditions and loading paths, pore pressures, temperatures and so on. The construction of constitutive laws is driven by two competing trends: the tendency for a unifying law describing the material behavior under the most general conditions, and the need for laws that can be used efficiently in practice. The latter incorporates both the possibility of calibration from available specimens and the user-friendliness of the law. For engineering purposes, the second tendency dominates and the engineer scientist has to evaluate the problem and select the most appropriate model to describe the dominant phenomena on hand. A review of various constitutive laws for geomaterials is presented with an emphasis on their application. Various examples of engineering problems and geomaterials with emphasis on petroleum applications are used to demonstrate that the law must be as simple as possible but not simpler.

Keywords: constitutive law, sandstone, chalk, water injection, core damage.

## INTRODUCTION

The constitutive laws that are presented describe the behavior of a geomaterial as a continuum. The perspective of the continuum material description can be appreciated in the following quote by Truesdell and Noll (1965):

> "...Widespread is the misconception that those who formulate continuum theories believe matter really is

> continuous, denying the existence of molecules. This is not so. Continuum physics presumes nothing regarding the structure of the matter. It confines itself to relations among gross phenomena, neglecting the structure of a material on a smaller scale. Whether the continuum approach is justified, in any particular case, is a matter, not for the philosophy or methodology of science, but for the experimental test..."
Continuum models are indeed widely used today for both research and engineering applications and provide useful solutions to technical problems. This of course is not to say that micromechanical and discontinuum models should be disregarded as purely scientific in nature. Without a doubt such models have shown a significant potential and in some cases they constitute the preferred way of modeling.

The literature on constitutive laws for geomaterials is rich in models. The reason is that a wide variety of geomaterials with often very dissimilar behaviors has to be appropriately described, such as the behavior of a sandstone, chalk, sand, limestone, shale, concrete, etc. In addition to the various geomaterials, the behavior of a particular material may depend significantly on the application at hand, that is on the boundary or loading conditions, the pore pressures, the pore fluids, the scale of the problem, the time periods involved, the temperature, etc.

There are two trends in constitutive modeling. The first is to make the model as general as possible in the expense of complexity and the other to make it as problem-specific as possible with the advantage of simplicity. The characteristics of a general constitutive model are:

- Various material types are included, such as elastic, plastic, creep, pore pressure, temperature, water saturation, etc.

- Various problems can be solved. The model can be used in various different problems such as wellbore stability, sand production, hydraulic fracturing, compaction, reservoir modeling, etc.

- The role of the engineer is limited.

- Extensive calibration is needed. The model requires the input of many material parameters, which are often difficult to obtain, as in the case of core material where enough test specimens are not available, or expensive to obtain, as large test series are required.

- Low efficiency. The model is cumbersome to use and often the results are difficult to interpret, as the important physical mechanisms are not always obvious.

On the other hand the characteristics of a problem-specific constitutive model are:

- Limited number of material types are included, such as elastic or elastic-plastic or creep or pore pressure or temperature or water saturation, etc.

- Limited number of problems can be solved, as for example wellbore stability or sand production or hydraulic fracturing or compaction or reservoir modeling, etc.

- The role of the engineer is significant since evaluation of the geomaterial and problem at hand and selection of the appropriate model are required.

- Limited calibration is needed, as only the absolutely necessary material parameters are required.

- Optimum efficiency is attained as non-important effects are excluded.

For engineering problems the problem-specific constitutive models have obvious advantages. Examples of such problems with emphasis on petroleum applications are given to demonstrate how to select or develop a constitutive model based on experimental and/or field data and observations. Section 1 describes a sandstone model for open hole stability and sand production problems, Section 2 a chalk model for reservoir compaction with water saturation effects, and Section 3 a sandstone model for core damage effects. Finally we present the concluding remarks.

# 1 SANDSTONE MODEL FOR OPEN HOLE STABILITY AND SAND PRODUCTION

Open hole stability and sand production problems are often encountered in sandstone reservoirs during hydrocarbon production. The reservoir rocks are usually soft, weak sandstones of considerable porosity in the range of 10-30%. In this problem, the rock is exposed to relatively low mean stresses by high deviatoric stresses due to the stress concentration around the hole.

A nonlinear elastic-plastic model has been developed to describe the near well behavior of these sandstones (Sulem et al., 1999). In brief, the model has an elasticity, which is nonlinearly stress-dependent, and a plasticity based on the flow theory. The plasticity is related to a linear Mohr-Coulomb yield function F with friction hardening and cohesion softening. The plastic dilatancy is described by a linear Mohr-Coulomb plastic potential function, different from the yield surface, which makes the flow rule not associative. This allows for the description of both compactive and dilatant behaviors under shear straining. On the other hand, the appropriate plastic behavior under hydrostatic straining has not been incorporated in the present model as plastic hydrostatic straining is not usually significant around wellbores and perforations. Throughout this section compression is taken negative, as it is more convenient for calculations. Tests results and model predictions are given with compression positive as more often used in rock and soil mechanics.

## 1.1 Model Description

The development of the model is based on behavior of soft, weak sandstones as observed in triaxial compression and extension experiments (Papamichos et al., 1999), which are relative experiments for obtaining the deviatoric behavior of a material. For nonlinear materials, the constitutive relations are formulated in the incremental form which, for time-independent behavior, is equivalent to the rate form employed in the following. A dot over a quantity denotes the rate (or increment) of the quantity.
The rate form of the constitutive stress-strain relationships
can be written as:

$$
\sigma_{ij} = C_{ijkl}^{ep} \epsilon_{kl} \quad (1)
$$

where σᵢⱼ and εᵢⱼ are the stress and strain rates, respec-
tively, and Cᵢⱼₖₗ the elastic-plastic tangent stiffness tensor,
which can be written as:

$$
C_{ijkl}^{ep} = C_{ijkl}^{e} - C_{ijkl}^{p} \quad (2)
$$

Figure 1 shows triaxial-compression test results on humid Red Wildmoor sandstone at three confining stresses. Starting from the elasticity of the material, the load-unload cycles show that the elasticity is nonlinear and depends on the stress. In order to keep the model as simple as possible, it is assumed that the Poisson’s ratio v is constant and thus the elastic shear modulus G must be a function of the stress, that is G = G(σᵢⱼ). The reversibility and path independency of the elastic strain energy requires that the stresses are obtained from the strain energy function W as (Chen and Han, 1988):

$$
\sigma_{ij} = \frac{\partial W}{\partial \varepsilon_{ij}} \qquad (3)
$$

With these assumptions and requirements, the elastic
shear modulus *G* must be a function of the equivalent stress
*σ*<sub>e</sub>, defined in terms of the isotropic (mean) stress *p* and the
shear-stress intensity *τ* as:

$$
\sigma_e = \sqrt{3(1 - 2\nu)p^2 + 2(1 + \nu)\tau^2} \quad (4)
$$

The stress invariants *p* and *τ* are defined as:

$$
p = I_1/3, \qquad \tau = J_2 \tag{5}
$$

where $I_1$ is the first stress invariant and $J_2$, the second
deviatoric stress invariant. The dependency of $G$ on the
equivalent stress $\sigma_e$, is described by the material function:

$$
G = G(\sigma_e) = \begin{cases} G_0 + \frac{(G_2 + G_3\sigma_e)\sigma_e}{1 + G_1\sigma_e} & \text{for } \sigma_e \le \sigma_e^{\text{peak}} \\ G_{\text{peak}} & \text{for } \sigma_e > \sigma_e^{\text{peak}} \end{cases} \quad (6)
$$

where the constants G₀ to G₃, Gₚeak, and σₑₚeak for the various
materials are obtained by fitting this function to the
experimental data. According to Equation (6), the elastic
shear modulus G increases from the initial value G = G₀ at
zero stress to G = Gₚeak at σₑ = σₑₚeak and remains constant
thereafter. The elastic tangent stiffness tensor Cᵢⱼₖ may then
be written as:

$$
C_{ijkl}^{e} = G \left[ \frac{2\nu}{1-2\nu} \delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk} \right. \\
\left. + \delta_{il} \delta_{jk} + 2(1+\nu) \frac{\sigma_{ij} \sigma_{kl}}{\sigma_{e}^{2}} \frac{dG/d\sigma_{e}}{G/\sigma_{e} - dG/d\sigma_{e}} \right] \quad (7)
$$

With respect to plasticity, the standard flow theory is used
where the plastic tangent stiffness tensor C<sub>ijkl</sub><sup>p</sup> is expressed as:

$$
C_{ijkl}^{p} = \frac{\langle 1 \rangle}{h} C_{ijmn}^{e} \frac{\partial Q}{\partial \sigma_{mn}} \frac{\partial F}{\partial \sigma_{st}} C_{stkl}^{e} \quad (8)
$$

with *F* and *Q* being the yield and plastic potential functions,
respectively. The hardening modulus *h* is given by:

$$
h = \frac{\partial F}{\partial \sigma_{ij}} C_{ijkl}^{e} \frac{\partial Q}{\partial \sigma_{kl}} - Q_{s} \left[ \frac{\partial F}{\partial \mu} \frac{d\mu}{dg^{p}} + \frac{\partial F}{\partial q} \frac{dq}{dg^{p}} \right] \quad (9)
$$

<figure><img src="image_2.png" /><figcaption>Figure 1</figcaption></figure>

Triaxial compression test results on humid Red Wildmoor sandstone. a: axial stress, and b: volumetric strain versus axial strain at three confining stresses (Papamichos et al., 1999).
and the following definition of the switch function $\langle 1 \rangle$ expressing the plastic loading criterion:

$$ \langle 1 \rangle = \begin{cases} 1, & \text{for } F = 0 \text{ and } \frac{\partial F}{\partial \sigma_{ij}} C_{ijkl}^e \varepsilon_{kl} > 0 \\ 0, & \text{for } F < 0 \text{ or } F = 0 \text{ and } \frac{\partial F}{\partial \sigma_{ij}} C_{ijkl}^e \varepsilon_{kl} \le 0 \end{cases} \quad (10) $$

Figure 2 shows the peak axial stress during triaxial compression at various confining stresses, which gives the shape of yield surface at peak stress. The yield surface is approximated with a straight line in order not to include additional nonlinearities in the model, such as for example a pressure dependency of the strength parameters. At higher confining stresses the curvature of the yield surface may be large enough, such that the exclusion of pressure-dependent effects on the strength parameters may not be possible. For the selection of the yield surface, reduced triaxial extension test data suggest that the material follows an intermediate behavior between the Mohr-Coulomb and the Drucker-Prager yield surfaces (Papamichos et al., 1999). Other analyses and comparisons between model predictions and test results suggest that the Mohr-Coulomb model is more appropriate for geomaterials. The Mohr-Coulomb is a model widely used in geomechanics and a lot of experience and strength parameters, such as cohesion and friction angle, are available.

Figure 1a shows that the material behaves plastically almost from the start of straining as demonstrated by the residual plastic strains during the unload-reload cycles. During the stress increase, the material hardens with the accumulated plastic shear-strain intensity $g^p$, and after the peak stress is reached, the material softens. The yield F and plastic potential Q surfaces are then functions of the stress

<figure><img src="image_3.png" /><figcaption>Figure 2
Triaxial compression test results on humid Red Wildmoor sandstone. Peak axial stress versus confining stress and straight line fit to the data (Papamichos et al., 1999).</figcaption></figure>

$\sigma_{ij}$ and the accumulated plastic shear-strain intensity $g^p$. The use of nonassociativity where the yield and plastic potential functions do not coincide is necessary in soft weak sandstones. The experimental results in Figure 1b show that the material initially compacts and then dilates as the peak stress is approached. With associative behavior, dilation is always predicted. For the Mohr-Coulomb criterion, F and Q can be written in terms of the stress invariants $p, \tau$, and $\theta$, as:

$$ \begin{aligned} F &= F(\sigma_{ij}, g^p) = \tau \left[ \cos \theta - \frac{\mu \sin \theta}{\sqrt{3}} \right] - \mu(q-p) = 0 \\ Q &= Q(\sigma_{ij}, g^p) = \tau \left[ \cos \theta - \frac{\beta \sin \theta}{\sqrt{3}} \right] + \beta p \end{aligned} \quad (11) $$

where the stress invariant $\theta$ is called the Lode angle and describes the relative magnitude of the intermediate principal stress. It is defined through the second $J_2$ and third $J_3$ deviatoric stress invariants as:

$$ \theta = \frac{1}{3} a \sin \left( - \frac{3\sqrt{3}J_3}{2J_2^{3/2}} \right) \quad (12) $$

The parameter $q = q(g^p)$ is the tension cut-off, and $\mu = \mu(g^p)$ and $\beta = \beta(g^p)$ are the mobilized friction and dilatancy coefficients, respectively. These are related to the mobilized friction angle $\phi_m$, dilation angle $\psi_m$, and cohesion $c$ through the expressions:

$$ \mu = \sin \phi_m, \quad \beta = \sin \psi_m, \quad q = c \cos \phi_m \quad (13) $$

The parameter $Q_s$ in Equation (9) is related to the plastic potential surface Q and it is given as a function of the stress invariants $J_2, J_3$ as:

$$ Q_s = 2 \sqrt{\frac{J_2 B_1^{Q^2} + 3 J_3 B_1^Q B_2^Q + J_2^2 B_2^{Q^2}}{3}} \quad (14) $$

where:

$$ B_1^Q = \frac{\partial Q}{\partial J_2} = \begin{cases} \frac{\cos \theta}{2\tau} \left[ \frac{1 + \tan \theta \tan 3\theta + \beta(\tan 3\theta - \tan \theta)}{\sqrt{3}} \right] & \text{if } \theta = -\frac{\pi}{6} \\ \frac{3 + \beta}{4\sqrt{3}\tau} & \text{if } \theta = -\frac{\pi}{6} \\ \frac{3 - \beta}{4\sqrt{3}\tau} & \text{if } \theta = \frac{\pi}{6} \end{cases} \quad (15) $$

$$ B_2^Q = \frac{\partial Q}{\partial J_3} = \begin{cases} \frac{\sqrt{3} \sin \theta + \beta \cos \theta}{2\tau^2 \cos 3\theta} & \text{if } \theta = \pm \frac{\pi}{6} \\ 0 & \text{if } \theta = \pm \frac{\pi}{6} \end{cases} $$

Finally, the plastic strain-rates $\varepsilon_{ij}^p$ are given by the flow rules:

$$ \varepsilon_{ij}^p = \psi \frac{\partial Q}{\partial \sigma_{ij}}, \quad \psi = \frac{\langle l \rangle}{h} \frac{\partial F}{\partial \sigma_{ij}} C_{ijkl}^e \varepsilon_{kl} \quad (16) $$
which result into the following expression for the plastic shear-strain intensity rate $g^p$:

$$g^p = \psi Q_s \quad (17)$$

For the elastic strain rates $\epsilon_{ij}^p$ it holds that:

$$\epsilon_{ij}^e = \epsilon_{ij} - \epsilon_{ij}^p \quad (18)$$

The model requires the calibration of the material functions $q = q(g^p)$, $\mu = \mu(g^p)$ and $\beta = \beta(g^p)$. For these sandstones, a model with both hardening and softening with plastic shear strain was used, as illustrated in Figure 3. Hardening was modeled by an increasing friction coefficient $\mu$, which reaches its maximum value at plastic shear strain $g^p = g_{peak}^p$ and remains constant thereafter. Softening was modeled by a tension cut-off $q$, which starts to decrease from its initial value $q_0$ after the friction coefficient has reached its peak value, and thus after the plastic shear strain $g^p$ has reached the value $g_{peak}^p$. In particular, the following functions were calibrated for each material:

$$q = \begin{cases} q_0 & \text{for } g^p \le g_{peak}^p \\ q_0 - q_1(g^p - g_{peak}^p) & \text{for } g^p > g_{peak}^p \end{cases}$$
$$\mu = \frac{3f}{2\sqrt{3+f}}, f = \begin{cases} f_0 + \frac{(f_2 + f_3 g^p)g^p}{1 + f_1 g^p} & \text{for } g^p \le g_{peak}^p \\ f_{peak} & \text{for } g^p > g_{peak}^p \end{cases} \quad (19)$$
$$\beta = \frac{3d}{2\sqrt{3+d}}, d = \begin{cases} \frac{d_2 + d_3 g^p(2+d_1g^p)}{(1+d_1g^p)} & \text{for } g^p \le g_{peak}^p \\ d_{peak} & \text{for } g^p > g_{peak}^p \end{cases}$$

where the constants $q_0$, $q_1$, $g_{peak}^p$, $f_0$ to $f_3$, $f_{peak}$, and $d_1$ to $d_3$, $d_{peak}$ are obtained by fitting the functions (19) on experimental data.

<figure><img src="image_4.png" /><figcaption>Figure 3

Motion of the Mohr-Coulomb yield surface in stress-space: (0-1) friction-hardening phase, (1-2) tension-softening phase (Sulem et al., 1999).</figcaption></figure>

The model can be calibrated through a minimum of two triaxial compression tests with load-unload cycles at various stress levels. Table 1 lists values of the various elastic and plastic material constants for humid Red Wildmoor sandstone.

<caption>TABLE 1</caption>

Humid Red Wildmoor sandstone: elastic and plastic material constants

<table><tbody><tr><td>v [-]</td><td>0.2</td><td>f<sub>0</sub> [-]</td><td>0.067148</td><td>q<sub>0</sub> [MPa]</td><td>3.82</td></tr><tr><td>G<sub>0</sub> [MPa]</td><td>349.17</td><td>f<sub>1</sub> [-]</td><td>109.88</td><td>q<sub>1</sub> [MPa]</td><td>25.0</td></tr><tr><td>G<sub>1</sub> [MPa<sup>-1</sup>]</td><td>0.042907</td><td>f<sub>2</sub> [-]</td><td>273.75</td><td></td><td></td></tr><tr><td>G<sub>2</sub> [-]</td><td>250.32</td><td>f<sub>3</sub> [-]</td><td>-8946.6</td><td>d<sub>1</sub> [-]</td><td>5.834</td></tr><tr><td>G<sub>3</sub> [MPa<sup>-1</sup>]</td><td>0.0</td><td>f<sub>peak</sub> [-]</td><td>0.94523</td><td>d<sub>2</sub> [-]</td><td>-0.78636</td></tr><tr><td>G<sub>peak</sub> [MPa]</td><td>6183.2</td><td>g<sub>peak</sub><sup>p</sup> [-]</td><td>0.0099069</td><td>d<sub>3</sub> [-]</td><td>90.453</td></tr><tr><td>&sigma;<sub>e</sub><sup>peak</sup> [MPa]</td><td>&infin;</td><td></td><td></td><td>d<sub>peak</sub> [-]</td><td>0.94523</td></tr></tbody></table>

## 1.2 Calculation Examples

The back analysis of the triaxial compression and reduced extension test results is shown in Figure 4. The comparison with the experimental results for both test types is very satisfactory despite the complexity of the soft weak sandstone behavior. The model was also used to predict the response of Red Wildmoor sandstone hollow cylinders loaded isotropically (Papamichos et al., 1996). In these tests a thick-wall hollow cylinder was loaded externally, while the internal hole remained unsupported. The external radial stress was equal to the axial stress. The internal hole and external deformations of the cylinder were measured. In Figure 5 experimental results for the external radial stress versus the internal and external tangential strains are compared with the model predictions. The nonlinear elastic-plastic model predicts very well the hollow cylinder response up to the point of failure of the inner hole. This is expected since after this point the continuum description of the material ceases.

The point of initial failure is predicted through the bifurcation theory and continua with microstructure. In addition to the nonlinear elastic-plastic model, predictions were made using a simpler model where the nonlinear elasticity was suppressed and an average value of elastic shear modulus was used. It can be seen in Figure 5 that the linear elastic-plastic model predictions are rather poor, which illustrates the need for a nonlinear elasticity model in this case.

For weakly cemented, unconsolidated sandstones the constitutive law may need to be enriched with additional features to model adequately the near well behavior of the rock. In the presented model, a linear Mohr-Coulomb model was adopted for relatively low mean stresses. In
<figure><img src="image_5.png" /><figcaption>Figure 4
Back analysis of a: triaxial compression, and b: reduced triaxial extension test results (Sulem et al., 1999).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5
Hollow cylinder stability. Nonlinear elastic-plastic versus linear elastic-plastic model simulations and comparison with experimental data.</figcaption></figure>

unconsolidated sandstones, however, a curved yield surface may be required even at low mean stresses. This may be modeled with a pressure-sensitive friction angle. The addition of a pressure cap may also be necessary. Another observation in these sandstones is their volumetric response upon shearing. Figure 6 shows the volumetric response of an unconsolidated reservoir sandstone tested in triaxial compression at three confining stresses. The experimental results clearly indicate a pressure-sensitive dilatancy, which may be necessary to be included in the constitutive law.

## 2 CHALK MODEL FOR RESERVOIR COMPACTION DUE TO WATER INJECTION

The mechanical behavior of weak, high-porosity chalks is important in many aspects of oil production such as reservoir compaction and solids production. Among several factors affecting the mechanical behavior of such rocks is the relative fluid saturation (gas/oil/water), which has an impact on the capillary forces between the grains. Capillarity arises from a fluid property known as surface tension, a phenomenon
<figure><img src="image_7.png" /><figcaption>Figure 6
Triaxial compression test results on an unconsolidated reservoir sandstone. Volumetric strain versus axial strain at three confining stresses, showing pressure-sensitive dilatancy.</figcaption></figure>

that occurs at the interface between different materials. For rocks and soils it occurs between surfaces of water (normally the wetting phase), mineral grains, and air, gas or oil (normally the nonwetting phase). At low degrees of saturation, water is bound in small pores and in the narrow spaces between grains forming capillary menisci. Low degrees of water saturation are often representative of the initial conditions in oil or gas reservoirs. It is noted that for some rocks, bridges of connate water are known to be stable for temperatures up to 300°C, posing problems in defining a dry reference state. The surface tension at the boundary between the water and the air in the adjoining voids pulls the grains together with a force known as contact pressure. The frictional resistance produced by the contact pressure has the same effect as if the grains were held with a certain amount of cohesion. At higher degrees of saturation, the voids between the grains are filled with water, the surface tension strongly decreases and the contact pressure vanishes. The magnitude of the capillary suction can reach the order of several MPa, depending on the radius of the capillary menisci, which is a function of the pore size. Another characteristic of chalks is their collapsible behavior. Chalks have a high porosity because they comprise of coccoliths, that is, skeletons of marine animals. In some chalks, it is found that capillary menisci formed by connate water are the only cohesive force that holds together the coccoliths. Upon elimination of the cohesive forces between the coccoliths, the chalk structure collapses with a significant reduction in porosity. Similar collapsible materials are the bulking structures observed in unconsolidated sands and sandstones.

An elastic-plastic model with water effects was developed to model the behavior of partially water-saturated chalks

(Papamichos et al., 1997). The model is based on experimental evidence from tests on Pietra Leccese chalk and incorporates Bishop's effective stress principle for the capillary pressure, a water-saturation-dependent elastic stiffness, and a Mohr-Coulomb yield surface with a pressure cap. The cohesion and pressure cap parameters in the yield surface are water-saturation-dependent to model the capillary effect on the strength parameters. Compression is taken negative in this section.

## 2.1 Experimental Evidence

The model is based on experimental evidence from triaxial and hydrostatic compression tests on Pietra Leccese chalk at various degrees of water saturation (Papamichos et al., 1997). The Pietra Leccese chalk is a high-porosity outcrop bio-calcarene from the south of Italy. Figure 7 shows uniaxial compression test results for specimens at various water saturations. The results show that the uniaxial compressive strength (UCS) decreases substantially with increasing water saturation. Load-unload cycles in these tests showed that for this material, the elasticity does not change significantly with the applied axial or confining stress. However, a significant decrease in the Young's modulus with increasing water saturation was observed. The Poisson's ratio on the other hand did not appear to be influenced significantly by the saturation. From the test data, an average elastic Young's modulus *E* for each saturation, and a constant average Poisson's ratio for all saturations were calculated.

In Figure 8, the peak axial stress obtained at uniaxial and triaxial compression tests is plotted versus the confining stress and saturation. From this plot and assuming a linear Mohr-Coulomb yield criterion for the peak axial stress, the

<figure><img src="image_8.png" /><figcaption>Figure 7
Uniaxial compression test results on Pietra Leccese chalk. Axial stress versus axial strain for various water saturations *S*.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8
Uniaxial and triaxial compression test results on Pietra Leccese chalk. Peak axial stress versus confining stress and saturation $S$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9
Hydrostatic compression tests on Pietra Leccese chalk showing the influence of water saturation on the pore collapse stress.</figcaption></figure>

peak friction angle $\phi$ and cohesion $c$ can be calculated for the two saturations where triaxial test data are available. The calculated values of peak friction angle and cohesion are tabulated in Table 2. From these limited test results, it appears that changes in the water saturation affect primarily the cohesion of the material, which reduces substantially in fully water-saturated specimens. This is attributed to the destruction, with increasing water saturation, of the capillary menisci, which act as cohesive bonds between the grains. On the other hand, the internal friction angle, which is primarily influenced by such factors as porosity, particle shape and size, and grain size distribution, appears not to be affected significantly.

<caption>TABLE 2</caption>

Pietra Leccese chalk: peak friction angle and cohesion at various water saturations calculated from triaxial compression test data

<table><thead><tr><td>Saturation</td><td>Peak friction angle (°)</td><td>Peak cohesion (MPa)</td></tr></thead><tbody><tr><td>0.047</td><td>29.1</td><td>10.34</td></tr><tr><td>1</td><td>27.5</td><td>6.46</td></tr></tbody></table>

Figure 9 shows the results from four drained hydrostatic compression tests performed on specimens at different water saturations, ranging from 0.013 to 1. During plastic yield under hydrostatic compression, the specimens undergo a characteristic structural change for collapsible materials which is referred to as pore collapse. This gradual structural change is characterized with a significant increase in compressibility. The results show that the elastic bulk modulus and the plastic yield stress decrease significantly with increasing water saturation. The tangent bulk modulus does not change substantially until plastic yield and thus the elasticity of Pietra Leccese can be considered at a first

approximation as stress-independent. The post-pore-collapse tangent bulk modulus is approximately the same in all tested specimens, and thus it can be considered as saturation-independent.

Three water-injection tests were also performed. The water-injection tests comprise of three phases. In the first phase, a specimen with initial saturation 0.013 or 0.019 is hydrostatically compressed under drained conditions until a prescribed hydrostatic stress level is reached. In the second phase, the hydrostatic stress is maintained constant, and water injection is performed by flowing few pore volumes of equilibrium water through the specimen at a constant flow rate. The flow rate is kept low to ensure no significant pore pressure built-up in the specimen. In the third phase, fluid flow is halted and the specimen is hydrostatically unloaded under drained conditions to zero stress. Figure 10 shows plots of the experimental results of water-injection tests at stresses above (test Plw17) and below (tests Plw19 and Plw19a), respectively, the hydrostatic plastic yield stress of the fully saturated material. The load-unload cycles for tests Plw19 and Plw19a show indeed an elastic behavior both before and after water injection. In the figure, the results from the hydrostatic tests without water injection on specimens with saturation 0.013, 0.019 and 1 are superimposed for comparison. During water injection the specimens experience initially a small dilation, followed by compaction. The total final compaction approaches the compaction obtained in the fully saturated specimen at the corresponding stress level. In the specimen water-injected at a stress level above the plastic yield stress of the fully saturated material, the additional water-injection-related compaction is significant, due to pore collapse during water injection. The initial dilation may be related to the elimination of the capillary suction during water injection.
<figure><img src="image_11.png" /><figcaption>Figure 10</figcaption></figure>

Water-injection tests on Pietra Leccese chalk at hydrostatic stress 10, 20 and 40 MPa. The results from the hydrostatic compression tests without water injection are included for comparison.

## 2.2 Model Description

In partially saturated media, the generalized effective stress $\sigma'_{ij}$ is given as:

$$ \sigma'_{ij} = \sigma_{ij} + p_{\text{gas}}\delta_{ij} - S(p_{\text{gas}} - p_{\text{w}})\delta_{ij} \quad (20) $$

where $\sigma_{ij}$ is the total stress, $S$ is the water saturation, $p_{\text{gas}}$ and $p_{\text{w}}$ the pore gas and pore water pressures, respectively, and $\delta_{ij}$ the Kronecker delta. By referring all pressures to the atmospheric pressure, $p_{\text{gas}}$ is zero when equal to the atmospheric, and Equation (20) can be written as:

$$ \sigma'_{ij} = \sigma_{ij} - S p_{cp} \delta_{ij} \quad (21) $$

where the pore-water capillary suction $p_{cp}$, defined as $p_{cp} = p_{\text{gas}} - p_{\text{w}}$, is a function of the saturation $S$, that is $p_{cp} = p_{cp}(S)$. Plastic yield is associated with an isotropic yield surface $F$, the size of which depends on the accumulated plastic strain parameter $\psi$ and the saturation $S$, that is $F(\sigma'_{ij}, \psi, S) = 0$. The total strain rate $\mathcal{E}_{ij}$ is decomposed into two parts:

$$ \mathcal{E}_{ij} = \mathcal{E}_{ij}^{e} + \mathcal{E}_{ij}^{p} \quad (22) $$

The elastic strain $\mathcal{E}_{ij}$ is related to the effective stress $\sigma'_{ij}$ by a linear elastic law as:

$$ \sigma'_{ij} = C_{ijkl}^{e} \mathcal{E}_{kl}^{e} \quad (23) $$

where the elastic stiffness tensor $C_{ijkl}^{e}$ is assumed to be a function of the saturation $S$, that is $C_{ijkl}^{e} = C_{ijkl}^{e}(S)$. Time

differentiation of Equation (23) yields a relation between the effective stress $\mathcal{E}'_{ij}$ rate and the elastic strain rate $\mathcal{E}_{ij}^{e}$:

$$ \mathcal{E}'_{ij} = C_{ijkl}^{e} \mathcal{E}_{kl}^{e} + D_{ij}^{e} \mathcal{S}, \quad D_{ij}^{e} = \frac{dC_{ijkl}^{e}}{dS} \mathcal{E}_{kl}^{e} \quad (24) $$

The plastic strain rate $\mathcal{E}_{ij}^{p}$ is obtained by a nonassociated flow rule as:

$$ \mathcal{E}_{ij}^{p} = \psi \frac{\partial Q}{\partial \sigma'_{ij}} \quad (25) $$

where $Q = Q(\sigma'_{ij}, \psi, S)$ is the plastic potential function, and $\psi$ is a scalar function. From the consistency condition:

$$ \mathcal{F} = \frac{\partial F}{\partial \sigma'_{ij}} \mathcal{E}'_{ij} + \frac{\partial F}{\partial \psi} \psi + \frac{\partial F}{\partial S} \mathcal{S} = 0 \quad (26) $$

and using Equations (22), (24) and (25), $\psi$ can be solved as:

$$ \begin{aligned} \psi &= \frac{1}{h} \frac{\partial F}{\partial \sigma'_{ij}} C_{ijkl}^{e} \mathcal{E}_{kl} + \frac{1}{h} \left( \frac{\partial F}{\partial S} + \frac{\partial F}{\partial \sigma'_{ij}} D_{ij}^{e} \right) \mathcal{S} \\ h &= \frac{\partial F}{\partial \sigma'_{ij}} C_{ijkl}^{e} \frac{\partial Q}{\partial \sigma'_{kl}} - \frac{\partial F}{\partial \psi} \end{aligned} \quad (27) $$

For plastic loading, $\psi$ is a nonnegative factor and $h\psi$ is always positive. It then follows from Equation (27) that the plastic loading criterion can be expressed in terms of the switch function $\langle 1 \rangle$ defined as:

$$ \langle 1 \rangle = \begin{cases} 1 & \text{for } F=0 \text{ and} \\ \frac{\partial F}{\partial \sigma'_{ij}} C_{ijkl}^{e} \mathcal{E}_{kl} + \left( \frac{\partial F}{\partial S} + \frac{\partial F}{\partial \sigma'_{ij}} D_{ij}^{e} \right) \mathcal{S} > 0 \\ 0 & \text{for } F < 0 \text{ or } F=0 \text{ and} \\ \frac{\partial F}{\partial \sigma'_{ij}} C_{ijkl}^{e} \mathcal{E} + \left( \frac{\partial F}{\partial S} + \frac{\partial F}{\partial \sigma'_{ij}} D_{ij}^{e} \right) \mathcal{S} \le 0 \end{cases} \quad (28) $$

The plastic strain rate can be calculated from the flow rule Equation (25), and Equation (27) as:

$$ \begin{aligned} \mathcal{E}_{ij}^{p} = \frac{\langle 1 \rangle}{h} \frac{\partial Q}{\partial \sigma'_{ij}} \frac{\partial F}{\partial \sigma'_{st}} C_{ijkl}^{e} \mathcal{E}_{kl} \\ + \frac{\langle 1 \rangle}{h} \frac{\partial Q}{\partial \sigma'_{ij}} \left( \frac{\partial F}{\partial S} + \frac{\partial F}{\partial \sigma'_{kl}} D_{kl}^{e} \right) \mathcal{S} \end{aligned} \quad (29) $$

The effective stress rate $\sigma'_{ij}$ can be determined from Equation (24), using Equations (22) and (29) as:

$$ \mathcal{E}'_{ij} = C_{ijkl}^{ep} \mathcal{E}_{kl} + D_{ij}^{ep} \mathcal{S} \quad (30) $$
where the tangent elasto-plastic stiffness $C_{ijkl}^{ep}$ and saturation $D_{ij}^{ep}$ tensors are given as:

$$C_{ijkl}^{ep} = C_{ijkl}^{e} - \frac{\langle 1 \rangle}{h} C_{ijmn}^{e} \frac{\partial Q}{\partial \sigma'_{mn}} \frac{\partial F}{\partial \sigma'_{st}} C_{stkl}^{e} \quad (31)$$

$$D_{ij}^{ep} = D_{ij}^{e} - \frac{\langle 1 \rangle}{h} C_{ijkl}^{e} \frac{\partial Q}{\partial \sigma'_{kl}} \left( \frac{\partial F}{\partial S} + \frac{\partial F}{\partial \sigma'_{st}} D_{st}^{e} \right)$$

Finally, the total stress rate $\boldsymbol{\sigma}_{ij}$ can be obtained by time differentiation of Equation (21) and substitution of $\boldsymbol{\sigma}'_{ij}$ from Equation (30), as:

$$\boldsymbol{\sigma}_{ij} = C_{ijkl}^{ep} \boldsymbol{\varepsilon}_{kl} + P_{ij} \boldsymbol{S}, \quad P_{ij} = \frac{d(S p_{cp})}{dS} \delta_{ij} + D_{ij}^{ep} \quad (32)$$

The employed elasticity is isotropic and linear with stress, in which case the elastic stiffness $C_{ijkl}^{ep}$ and saturation $D_{ij}^{e}$ tensors are given by the expressions:

$$C_{ijkl}^{e} = G(S) \left[ \frac{2\nu}{1 - 2\nu} \delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk} \right] \quad (33)$$

$$D_{ij}^{e} = \frac{1}{G} \frac{dG}{dS} \boldsymbol{\sigma}'_{ij}$$

where G is the elastic shear modulus and v the Poisson's ratio. Based on the experimental data for the Pietra Leccese chalk, the dependency of the elasticity on the saturation S is reflected only in the shear modulus G, which is assumed as a function of the saturation.

The employed yield surface is a modified Mohr-Coulomb yield function to include a pressure cap. It is described by the expression:

$$F(\boldsymbol{\sigma}'_{ij}, \mu, q, p_c, S) = \tau \left[ \cos \theta - \frac{\mu}{\sqrt{3}} \sin \theta \right] - \mu(q - p) \sqrt{\frac{p + p_c}{q + p_c}} = 0 \quad (34)$$

where three invariants, termed the isotropic (mean) effective stress p, the shear stress intensity $\tau$, and the Lode angle $\theta$, have been introduced for the effective stress $\boldsymbol{\sigma}'_{ij}$:

$$p = \frac{1}{3} I_1, \quad \tau = \sqrt{J_2}, \quad \theta = \frac{1}{3} a \sin \left( - \frac{3\sqrt{3}}{2} \frac{J_3}{J_2^{3/2}} \right) \quad (35)$$

with $I_1$ the first invariant of $\boldsymbol{\sigma}'_{ij}$, and $J_2, J_3$ the second and third invariants of the deviatoric stress $S_{ij}$, respectively. In Equation (34) the hardening/softening parameters q, $p_c$ and $\mu$ represent the size of the yield surface. The parameters q and $p_c$ are identified as the intercepts of the yield surface F with the positive (tension cut-off) and negative (pressure cap) p-axis, respectively. A comparison of Equation (34) with the Mohr-Coulomb yield function shows that the parameter $\mu$

conceptually replaces the sinus of the friction angle $\phi$, and that the modified yield surface has been augmented with a square root term for the mean pressure p and a new parameter $p_c$ related to the pressure cap. Figure 11 compares the Mohr-Coulomb and the modified Mohr-Coulomb yield surfaces in the $\tau-p$ plane for $\mu = 0.6, q = 5$ MPa, $p_c = 10$ MPa, and conditions of conventional triaxial compression where $\theta = \pi/6$.

<figure><img src="image_12.png" /><figcaption>Figure 11
Illustration in the $\tau-p$ plane of the original Mohr-Coulomb and the Mohr-Coulomb with a pressure cap yield surfaces.</figcaption></figure>

The corresponding plastic potential function to Equation (34) is:

$$Q(\boldsymbol{\sigma}'_{ij}, \beta, q, p_c, S) = \tau \left[ \cos \theta - \frac{\mu}{\sqrt{3}} \sin \theta \right] - \beta(q - p) \sqrt{\frac{p + p_c}{q + p_c}} \quad (36)$$

where the parameter $\beta$ replaces the sinus of the dilatancy angle in the Mohr-Coulomb model. The hardening parameters $\mu$, $\beta$ are expressed as functions of the accumulated plastic shear strain $g^p$, that is $\mu = \mu(g^p)$, $\beta = \beta(g^p)$. The parameter q is expressed as a function of the saturation S, while the parameter $p_c$ as a function of the saturation and the accumulated plastic volumetric strain $v^p$, that is $q = q(S)$, $p_c = p_c(S, v^p)$. The plastic shear $g^p$ and volumetric $v^p$ strains are defined as:

$$g^p = \int g^p dt, \quad g^p = \sqrt{2\varepsilon_{ij}^p \varepsilon_{ij}^p} \quad (37) \\ v^p = \int v^p dt, \quad v^p = \varepsilon_{kk}^p$$

where $\varepsilon_{ij}^p$ is the plastic deviatoric strain.
## 2.3 Calculation Examples

The calibration procedure and the model parameters for the Pietra Leccese chalk can be found in Papamichos et al. (1997). The hydrostatic compression tests on the specimens with saturation $S = 0.013, 0.019, 0.047$, and 1 were back-analyzed with the developed model and the results are presented and compared with the experimental data in Figure 12. The simulations compare well with the experimental data, showing that the model can simulate properly the saturation-dependent elastic stiffness during the initial part of the test, the saturation-dependent hydrostatic stress at plastic yield, and reasonably well the pore collapse behavior of the Pietra Leccese chalk. The back analysis of the water-injection test at 40 MPa and the comparison with the experimental data are presented in Figure 13. The specimen enters the plastic regime of deformation and experiences pore collapse during the water-injection phase. The simulation compares well with the experimental data, and the model can simulate properly the saturation-dependent elastic stiffness during the initial loading part of the test, the additional compaction during the water-injection phase at both the elastic and plastic regimes, and the elastic stiffness during the unloading part of the test. In the model, the specimen straining during water injection comprises of three parts:

- a dilatant part due to the elimination of the capillary suction, which essentially results in a reduction of the effective compressive stress;

- a compactive elastic part due to the reduction of the elastic stiffness with increasing water saturation;

- a compactive plastic part due to plastification from the softening of the pressure cap with increasing water saturation.

<figure><img src="image_13.png" /><figcaption>Figure 12</figcaption></figure>

Theoretical simulations and experimental results of hydrostatic compression tests on specimens with various water saturations $S$.

<figure><img src="image_14.png" /><figcaption>Figure 13</figcaption></figure>

Theoretical simulations and experimental results of hydrostatic compression tests with water injection at hydrostatic compressive stress 40 MPa.

## 3 SANDSTONE MODEL FOR CORE DAMAGE

Core damage is a permanent alteration of the rock properties as a result of drill-out and retrieval of the core from the *in situ* environment to the surface. For reservoir rocks, stress release during coring is thought to be a major core damage mechanism. Holt *et al.* (1998) have investigated core damage in reservoir sandstones using a synthetic sandstone, analogous to the reservoir rock with respect to rock mechanical and petrophysical properties and texture. The synthetic sandstone is created under *in situ* stress conditions to simulate the diagenesis process. Comparative tests are then performed to simulate virgin compaction behavior, and the behavior of an unloaded and reloaded simulated core.

The stress paths followed in the virgin compaction tests and the simulated core compaction tests are shown in Figure 14. In the virgin compaction tests, sand was loaded to the *in situ* stress conditions where cementation took place. After cementation, the specimen was compacted under oedometric $K_0$ conditions. In the simulated core compaction tests, after cementation at *in situ* stress conditions, the specimen was unloaded to simulate the coring process. Subsequently the specimen was loaded to the *in situ* stress conditions and an oedometric $K_0$ compaction test was performed to simulate a test on a cored specimen. Figure 15 shows a comparison between the virgin compaction and the simulated core compaction for this sandstone. Figure 15a shows only the $K_0$ part of the tests. The virgin compaction curves are nonlinear with a higher initial stiffness. The stiffness during virgin compaction matches the stiffness of the core compaction above 45 MPa. Thus, if core data are used directly for reservoir rock compaction estimation, the
<figure><img src="image_15.png" /><figcaption>Figure 14</figcaption></figure>

Stress paths in compaction tests. a: virgin $K_0$ compaction, and b: coring simulation, reloading to *in situ* conditions and $K_0$ compaction (Holt *et al.*, 1998).

<figure><img src="image_16.png" /><figcaption>Figure 15</figcaption></figure>

Compaction behavior of virgin versus simulated core of a synthetic reservoir sandstone (Holt *et al.*, 1998). a: $K_0$ phase of the tests, and b: all tests.

initial compaction will be overestimated. In Figure 15b, the unloading during the coring simulation is also included. The two specimens behave exactly the same above an axial stress of 45 MPa.

The test results show that during coring, unrecoverable damage takes place with large plastic strains. This can be modeled with an elastoplastic model with a pressure cap to model correctly the behavior under uniaxial compaction. A way to describe the core damage process is a model with two yield surfaces. Both yield surfaces have a pressure cap. Such models are referred to as overlay or mechanical sublayer models. The one yield surface is for the sand matrix and the other for the cement. The sand matrix at *in situ* stresses is loaded and thus, during coring, it simply unloads. The yield surface for sand is centered at zero stresses and extends such

that the *in situ* stress conditions lie on it. During virgin $K_0$ loading, the sand matrix is further loaded and the yield surface extends further such that the current stress conditions lie always on the yield surface. This is illustrated schematically in Figure 16.

The cement at *in situ* stresses is not loaded and thus during either coring or virgin $K_0$ loading, it loads, and after an elastic region it yields. The yield surface for the cement is centered at the *in situ* stress conditions and extends such that it gives an initial elastic region upon reduction of the stresses for coring simulation or increase of the stresses for virgin $K_0$ loading. Yielding leads to the destruction of the cement cohesion, such that the cement yield surface reduces in size during stress changes, which can be either due to coring or virgin $K_0$ loading. Upon reloading to *in situ* stresses,
<figure><img src="image_17.png" /><figcaption>Figure 16
Yield surfaces *F* for the sand matrix and the cement and stress paths during virgin and core compaction.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17
Simulation of virgin $K_0$ compaction and coring $K_0$ compaction.</figcaption></figure>

nonlinear but elastic behavior is observed. After that the material continues to yield since the cement has been destructed. Essentially destruction of the cement occurs either by loading or by unloading from *in situ* conditions, such that after the cement destruction the two materials behave exactly the same.

The proposed two-material model can also be viewed as a model with internal, often referred to also as residual,

stresses. The internal stresses result from the fact that the two materials, the sand matrix and the cement, act in opposite ways in the case of unloading, that is coring. Internally viewed the compression in the sand matrix is counter balanced by tension in the cement, resulting in zero external stress.

Figure 17 shows a simulation of virgin $K_0$ compaction and coring $K_0$ compaction. The model results compare favorable with the experimental results in Figure 15a showing the potential of the model to predict core damage in reservoir sandstones.

## CONCLUSIONS

The presented modeling examples show that the behavior of geomaterials is often too complex to generalize the models. Enough testing material is often not available for a calibration of a general constitutive model at a reasonable cost. General models lead also to cumbersome computations, which often obscure the physics of the problem at hand. The role of the engineer becomes thus important, as he is required to identify and model first-order effects.

## REFERENCES

Chen, W.F. and Han, D.J. (1988) *Plasticity for Structural Engineers*, Springer-Verlag, 173.

Holt, R.M., Brignoli, M., Kenter, C.J., Meij, R. and Schutjens, P.M.T.M (1998) From Core Compaction to Reservoir Compaction: Correction for Core Damage Effects. *SPE/ISRM 47263*. *Proc. Eurock '98*, 1, 311-320.

Papamichos, E., Brignoli, M. and Santarelli, F.J. (1997) An Experimental and Theoretical Study of a Partially-Saturated Collapsible Rock. *Mech. Cohes.-Frict. Mater.*, 2, 3, 251-278.

Papamichos, E., Tronvoll, J., Skjærstein, A., Unander, T.E., Vardoulakis, I. and Sulem, J. (1996) The Effect of Plane-Strain and Isotropic Loading in Hollow-Cylinder Strength. *Proc. Eurock '96*, Torino, Barla, G. (ed.), 1, Balkema, Rotterdam, 1997-204.

Papamichos, E., Tronvoll, J., Vardoulakis, I., Labuz, J.F., Skjærstein, A., Unander, T.E. and Sulem, J. (1999) Constitutive Testing of Red Wildmoor Sandstone. *Mech. Cohes.-Frict. Mater.*, 4.

Sulem, J., Vardoulakis, I., Papamichos, E., Oulahna, A. and Tronvoll, J. (1999) Elasto-Plastic Behaviour of Red Wildmoor Sandstone. *Mech. Cohes.-Frict. Mater.*, 4.

Truesdell, C. and Noll, W. (1965) The Non-Linear Field Theories of Mechanics. In *Encyclopedia of Physics*, Flügge, S. (ed.), III/3, Sect. 3, Springer-Verlag.