Soil Mechanics – Page 3 – vulcanhammer.net

Academic Issues, Geotechnical Engineering, Soil Mechanics

Deriving and Solving the Equations of Consolidation

Posted on 1 November 20224 November 2022 by Don Warrington

In an earlier post we discussed consolidation settlement. For situations where soil a) decreases its volume due to rearrangement of the particles and b) does so over a relatively long period of time due to difficulties in expelling pore water pressures, we need to know how long it takes to reach maximum settlement, in addition to know what that settlement is and what deflections we might achieve along the way.

This is generally a two part process: a) determining the dissipation of pore water pressure and b) determining the amount of settlement at a given time. This post will focus on the first part of the process; we will deal with the second in a later post. This is a well-worn path in geotechnical engineering but hopefully this derivation is a little simpler than others.

Developing the Differential Equation

Relating Porosity to Strain

Let us begin by considering the system above (from Tsytovich (1976)), a uniformly loaded soil layer which is saturated (always) and clay (usually.) The first thing we need to note is that, while the diagram uses z as the variable of length in the vertical direction, from this post we will use the variable x (sorry for the confusion.) The second thing is that we assume the water and the solids to be incompressible. The third thing is that all the changes that take place do so because of changes in the voids where the water is resident. We first note the definition of porosity as

$n={\frac {V_{{v}}}{V_{{t}}}}$ (1)

where

n = porosity
V_v = volume of voids
V_t = total volume

That being the case, the relationship between the porosity of the soil (due to changes in the volume of the voids) and the flow rate of the water can be expressed as

${\frac {\partial }{\partial x}}q(x,t)=-{\frac {\partial }{\partial t}} {\it n}(x,t)$ (2)

We can envision a differential volume having a height x and an area A. Since the problem is one-dimensional, the areas cancel out and the ratio of the void height to the total height is

$n={\frac {x_{{v}}}{x_{{t}}}}$ (3)

where

x_v = height of the voids
x_t = height of the solids

We want to determine the change in porosity from some state 0 to some state 1, just as we did with void ratio in this post. That change can be expressed as follows:

${\frac {x_{{{\it v0}}}}{x_{{{\it t0}}}}}-{\frac {x_{{{\it v1}}}}{x_{{{\it t0}}}-x_{{{\it v0}}}+x_{{{\it v1}}}}}$ (4)

We can assume that the change in void volume x_v0 – x_v1 << x_t0, in which case Equation (4) can be simplified to

$\Delta{{{\it n}}}={\frac {x_{{{\it v0}}}-x_{{{\it v1}}}}{x_{{{\it t0}}}}} = \frac{\Delta x}{x_t}$ (5)

Now we can say, for the small increments we are dealing with here, that the change in strain is equal to the change in porosity,

$\Delta n = \Delta \epsilon$ (6)

From this and our previous post, we can thus use the change in strain to come to the following:

$\Delta n = m_v \Delta \sigma_x$ (7)

Stating this differentially,

${\frac {\partial }{\partial t}}{\it n}(x,t)=m_{{v}}{\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (8)

In this way we emphasise that both porosity and uniaxial stress are functions of depth and time. We now combine Equations (2) and (8) to yield

${\frac {\partial }{\partial x}}q(x,t)=-m_{{v}}{\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (9)

Including Permeability

Darcy’s Law (or more properly d’Arcy’s Law) states that

$q(x,t) = -k{\frac {\partial }{\partial x}}H(x,t)$ (10)

where k is the coefficient of permeability and H is the hydraulic head. We then combine Equations (9) and (10) to obtain

$k{\frac {\partial ^{2}}{\partial {x}^{2}}}H(x,t)=m_{{v}}{\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (11)

Noting that

$H(x,t) = \frac{u(x,t)}{\gamma_w}$ (12)

where $\gamma_w$ is the unit weight of water and $u(x,t)$ is the excess pore water pressure generated by the decrease in porosity, we substitute this with the result of

${\frac {k{\frac {\partial ^{2}}{\partial {x}^{2}}}u(x,t)}{{\it \gamma_w }}}=m_{{v}}{\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (13)

Putting It All Together

We now define the parameter

$c_v = \frac{k}{m_v \gamma_w}$ (14)

At this point we need to stop and make an observation: this whole process wouldn’t be worth too much if $c_v$ wasn’t constant (or reasonably so.) The reason it is is that the coefficient of permeability $k$ and the coefficient of volume compressibility $m_v$ both decrease as the void ratio/porosity decrease, and do so at roughly the same rate. That being the case, we substitute Equation (14) into Equation (13) to have at last

$c_{{v}}{\frac {\partial ^{2}}{\partial {x}^{2}}}u(x,t)={\frac {\partial }{\partial t}}\sigma_{{x}}(x,t)$ (15)

Tidying up the Physics, Governing Equation, Boundary and Initial Conditions, and the Solution

It should be evident that getting to Equation (15) was a major triumph for geotechnical theory. But it’s also evident that there’s one glaring problem: the dependent variable on the left hand side is not the same as the one on the right. It is here that we need to clarify some assumptions behind our equations.

A basic assumption in consolidation theory is that, when the load at the surface is applied, all of this additional load is initially borne by the pore water pressure. Because of the aforementioned permeability, the water will want to “head for the exits,” i.e., the permeable boundaries of the layer being compressed. When the particles have rearranged themselves and the excess pore water has been squeezed out, the settlement should stop (until secondary compression kicks in.) During this process the load on the pore water is being progressively transferred to the soil particles until consolidation has stopped (which, in theory, it never does, as we will see) and the load is completely handed off to the particles.

That being the case, the consolidation equation (15) should be rewritten as

$c_{{v}}{\frac {\partial ^{2}}{\partial {x}^{2}}}u(x,t)={\frac {\partial }{\partial t}}u(x,t)$ (16)

At this point we should invoke the effective stress equation

$\sigma_x(x,t) = p - u(x,t)$ (17)

where p is the applied pressure, and do the following to determine the time and distance history of the vertical stress $\sigma_x(x,t)$ :

Determine the solution of Equation (16), which will be our governing equation.
Determine the initial and boundary conditions for the problem.
Solve the problem for $\sigma_x(x,t)$ using Equation (17).

There are several ways to solve the governing equation for this problem. Verruijt employed Laplace Transforms to accomplish this; these are very useful, as was demonstrated by Warrington (1997) for the wave equation. For this analysis, we will use separation of variables, a method which is not as fundamental as one might like (this is a demonstration of a method that is) but which is easier to follow than most of the others.

The separation of variables begins by assuming the solution is as follows:

$u(x,t) = X(x)T(t)$ (18)

This means that the solution is a product of two functions, one of time and the other of distance, and that those functions are separate one from another. Substituting this into Equation (16) and rearranging a bit yields

${\frac {{\frac {d^{2}}{d{x}^{2}}}X(x)}{X(x)}}={\frac {{\frac {d}{dt}}T(t)}{{\it c_v}\,T(t)}}$ (19)

Since the left and right hand sides are equal to each other, they are equal to a third variable, which we will designate as $\beta^2$ . In reality they are equivalent to the eigenvalue $\lambda = \beta^2$ , but the utility of the squared term will become evident. In any case,

${\frac {{\frac {d^{2}}{d{x}^{2}}}X(x)}{X(x)}}=-{\beta}^{2}$ (20a)

${\frac {{\frac {d}{dt}}T(t)}{{\it c_v}\,T(t)}}=-{\beta}^{2}$ (20b)

The solutions to these equations are, respectively,

$X(x)={\it C_1}\,\cos(\beta\,x)+{\it C_2}\,\sin(\beta\,x)$ (21a)

$T(t)={e^{-{\it c_v}\,{\beta}^{2}t}}{\it C_3}$ (21b)

We need to pause and consider the boundary conditions. As the problem is shown at the top of the article, the boundary conditions are as follows:

$X(0) = 0$ (22a)

$X'(h) = 0$ (22b)

Equation (22a) represents a Dirichelet boundary condition and Equation (22b) represents a Neumann boundary condition. While the equation for $X(x)$ can certainly be solved for this boundary condition, a simpler way would be to do the following:

Mirror the problem shown above at x = h so that you have two permeable boundaries.
You now have a layer of 2h thickness with Dirichelet boundary conditions on both sides. At the centre of this new layer, there is no flow; the water above it flows upward, and the water below flows downward.
Problems in the field can have either one or two permeable boundaries. The distance h is NOT the thickness of the layer but the longest distance the trapped pore water must travel to escape. Confusing h with the thickness of the layer is a common mistake and should be avoided at all costs.

That said, the second boundary condition is now

$X(2h) = 0$ (22c)

Returning to Equation (21a), because of Equation (22a) $C_1 = 0$ as the cosine is by definition unity at this point. If we then set x = 2h and X(2h) = 0, for real, non-zero values of $C_2$ the boundary condition can be satisfied if and only if

$\beta = 1/2\,{\frac {n\pi }{h}}$ (23)

This is the square root of the eigenvalues. Any integer value of n > 0 is valid for this, and we will have recourse to them all to produce a complete orthogonal set (Fourier Series) to solve the problem.

This change will affect both Equations (21a) and (21b). Combining constants into the coefficient $B_n$ , substituting the results into Equation (18) and making it an infinite sum for the Fourier series yields

$u(x,t) = {\it B_{n}}\,\sin(1/2\,{\frac{n\pi\,x}{h}}){e^{-1/4\,{\frac{{\it c_v}\,{n}^{2}{\pi}^{2}t}{{h}^{2}}}}}$ (24)

At this point we need to consider our initial conditions $f(x)$ . Although many different distributions of initial pore pressure are possible, the simplest one–and the one most commonly used–is a uniform pressure p, which is the same as the surface pressure compressing the layer. For a complete orthogonal set, the coefficients $B_n$ can be computed by (Kreyszig (1988)) as

${\it B_n}=\int _{0}^{2\,h}f(x)\sin(1/2\,{\frac {n\pi \,x}{h}}){dx}{h}^{-1}$ (25)

Substituting $f(x) = p$ and performing the integration,

$B_n = 2\,{\frac{p\left(1-\cos(n\pi)\right)}{n\pi}}$ (26)

Substituting this into Equation (24) and taking the complete sum, we have at last

$u(x,t)=\sum_{n=1}^{\infty}2\,p\left(1-\cos(n\pi)\right)\sin(1/2\,{\frac{n\pi\,x}{h}}){e^{-1/4\,{\frac{{\it c_v}\,{n}^{2}{\pi}^{2}t}{{h}^{2}}}}}{n}^{-1}{\pi}^{-1}$ (27)

More simply we can say that

$u(x,t)=\frac{4p}{\pi}\left(\sin(1/2\,{\frac{\pi\,x}{h}}){e^{-1/4\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}+1/3\,\sin(3/2\,{\frac{\pi\,x}{h}}){e^{-9/4\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}+1/5\,\sin(5/2\,{\frac{\pi\,x}{h}}){e^{-{\frac{25}{4}}\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}\cdots\right)$ (28)

The vertical pressure on the soil skeleton is determined by combining Equations (17) and (28) to yield

$\sigma_{x}(x,t)=p\left(1-\frac{4}{\pi}\left(\sin(1/2\,{\frac{\pi\,x}{h}}){e^{-1/4\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}+1/3\,\sin(3/2\,{\frac{\pi\,x}{h}}){e^{-9/4\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}+1/5\,\sin(5/2\,{\frac{\pi\,x}{h}}){e^{-{\frac{25}{4}}\,{\frac{{\it c_v}\,{\pi}^{2}t}{{h}^{2}}}}}\cdots\right)\right)$ (29)

It is interesting to note that, although Equation (27) includes all positive non-zero values of n, only the odd ones end up in Equations (28) and (29). For even values of n, $B_n = 0$ .

Conclusion

At this point we have derived the equations of pressure dissapation for consolidation settlement. In a future post we will deal with the second part of the problem, namely how much of the total anticipated settlement has taken place at any given time from initial loading.

Reference

Kreyszig, E. (1988) Advanced Engineering Mathematics. Sixth Edition. New York: John Wiley and Sons

Geotechnical Engineering, Soil Mechanics

The “New” NAVFAC DM 7.1 (Soil Mechanics) is Now in Print

Posted on 19 October 2022 by Don Warrington

Recently, NAVFAC revised their half-century old NAVFAC DM 7.1 and have released it to the public. Now we offer this classic (formally designated as UFC 3-220-10) in print format. Click here or on the photo at the right to order.

The chapters are as follows:

IDENTIFICATION AND CLASSIFICATION OF SOIL AND ROCK
FIELD EXPLORATION, TESTING, AND INSTRUMENTATION
LABORATORY TESTING
DISTRIBUTION OF STRESSES
ANALYSIS OF SETTLEMENT AND VOLUME EXPANSION
SEEPAGE AND DRAINAGE
SLOPE STABILITY
CORRELATIONS FOR SOIL AND ROCK

As was the case before, there is a wealth of information here, updated and expanded. Now you can order your copy in print today!

Geotechnical Engineering, Soil Mechanics

The Problem of Size: Gazetas and Stokoe (1991) Revisited

Posted on 28 July 202214 October 2022 by Don Warrington

In a recent exchange Dr. Mark Svinkin, who has contributed several well-read articles to this site, pointed out that he had commented on a paper by Gazetas and Stokoe (1991.) The paper, Dr. Svinkin’s comments and their response can be found here:

Although this research was done a long time ago, it’s worth revisiting because of the issue that Dr. Svinkin brings up: the issue of size, that it’s not a straightforward business to extrapolate the results of model tests in controlled environments to full-scale foundations in actual stratigraphies.

In my fluid mechanics laboratory course, I discuss the issue of dynamic similarity, how one can take an airfoil or other flying object on a small scale and, using things such as the Reynolds Number, extrapolate those results to full-scale aircraft. This has proven very useful in the development of aircraft, especially before (and even long after) the development of simulation using computational fluid dynamics.

With geotechnical engineering, it has not been quite as simple. Attempts to use things such as centrifuge testing have not been as successful as, say, wind tunnel testing has been for aeronautics. Part of the problem is, as I like to say, that geotechinical engineering is not non-linear in the same sense as fluids are. Another problem is that the earth is not as homogeneous as the atmosphere, even when altitude and weather effects are considered (and these influence each other in the course of events.) But underneath all of this there are some fundamental issues that have complicated the issue of foundation size, and Dr. Svinkin points this out. My intent is to amplify on that and remind people that these issues are still relevant.

Dr. Svinkin points out the following figure from Tsytovich (1976.) I’ve referenced this text in several recent posts. Tsytovich looks at many problems in soil mechanics differently from our usual view in this country, and his perspective is frequently insightful. (An excellent example is here.) In this diagram he shows the effect of basic foundation size on the settlement of the foundation, and Tsytovich’s own explanation of this follows:

Relationship between settlement of natural soils and dimensions of loading area (from Tsytovich (1976).) The variable F is the area of the foundation, thus the square root of F is the basic dimension of the foundation and, in the case of square foundations, the exact dimension b or B (see below.)

Thus, Fig. 90 shows a generalized curve of the average results of numerous experiments on studying the settlements of earth bases (at an average degree of compaction) for the same pressure on soil but with different areas of loading. Three different regions may be distinguished on the curve: I —the region of small loading areas (approximately up to 0.25 m²) where soils at average pressures are predominantly in the shear phase, with the settlement being reduced with an increase of area (just opposite to what is predicted by the theory of elasticity for the phase of linear deformations); II — the region of areas from 0.25-0.50 m² to 25-50 m² (for homogeneous soils of medium density, and to higher values for weak soils), where settlements are strictly proportional to $\sqrt{F}$ and at average pressures on soil correspond to the compaction phase, i.e., are very close to the theoretical ones; and III — the region of areas larger than 25-50 m², where settlements are smaller than the theoretical ones, which may be explained by an increase of the soil modulus of elasticity (or a decrease of deform ability) with an increase of depth. For very loose and very dense soils these limits will naturally be somewhat different.

The data given can be used for establishing the limits of applicability of the theoretical solutions obtained for homogeneous massifs to real soils, which is of especial importance in developing rational methods of calculation of foundation settlements.
From Tsytovich (1976)

Although much of the discussion centred on Tsytovich and Barkan (1962,) there is evidence elsewhere to underscore this problem, which Tsytovich sets forth in a very succinct manner.

To begin with, let us consider the basic equation for elastic settlement, which was discussed in this post, and is as follows:

$s = \frac{\omega p b (1-\nu^2)}{E}$ (1)

where

$s =$ settlement of the foundation at the point of interest
$\omega = I =$ influence factor, given in the table below
$p =$ uniform pressure on the foundation
$b = B =$ smaller dimension of rectangle or dimension of square side
$\nu =$ Poisson’s Ratio of the soil
$E =$ Modulus of elasticity of the soil

The values of $\omega$ are shown below.

Similar values can be found in both the Soils and Foundations Manual and NAVFAC DM 7.

It is clear that, once one is past the basic soil properties and the pressure applied on the foundation, the settlement is proportional to the basic dimension of the foundation, which is exactly what is taking place in Region II. This is also why the bulk modulus of the soil is not a basic soil property, as I discuss in this lecture. When we consider plate load tests, we must correct them for the difference between the size of the test plate and the size of the foundation, as this slide presentation shows.

Since we are dealing with foundation dynamics, one item that seems to have fallen out of the whole discussion is that of Lysmer (1965). Lysmer’s Analogue, which reduces the response of a soil under the foundation to a simple spring-damper-mass system, defines the spring constant as follows:

$K = \frac{4Gr}{1-\nu}$ (2)

where $K$ is the spring constant of the soil and $r$ is the foundation’s radius. If we break it down further, as is done in Warrington (1997,) and develop a unit area spring constant under the foundation, we have

$k = \frac{4Gr}{F(1-\nu)}$ (3)

where $k$ is the equivalent unit area spring constant under the soil. Equation (3) in particular shows that, for a given unit load on a foundation, the static portion of the reaction is inversely proportional to the basic size of the foundation. (The unit damping constant is actually independent of the area for round foundations.)

These results show that, while the effect of size may differ from one model to the next, it cannot be overlooked in any attempt to extrapolate physical model tests of any kind to actual use. This effect is further complicated by variations in shear modulus due to either strain softening, layered stratigraphy, effective stresses or other factors. The effect of the stratigraphy is further magnified by the fact that larger foundations have larger “bulbs of influence” into the soil and thus layers that smaller foundations would not interact with become significant with larger ones.

“Sand box” tests have other challenges. While they attempt to simulate a semi-infinite space, reflections from the walls of the box are inevitable, especially with periodic loads such as were present in this test. These challenges were documented in the original study. (An interesting study using another one of these boxes is that of Perry (1963).)

The failure of geotechnical engineering to adequately resolve the size issue, both in terms of design and in terms of using laboratory data to simulate full-scale performance, remains a frustration in geotechinical engineering. Hopefully other types of models will help move things forward, along with advances in our understanding of soil behaviour and our ability to replicate it both experimentally and numerically.

Academic Issues, Soil Mechanics

Manual For The Systematic Study Of The Regime Of Underground Waters – Altovskij, Konopljantsev — Mir Books

Posted on 13 June 2022 by Don Warrington

In this post, we will see the book Manual For The Systematic Study Of The Regime Of Underground Waters edited by M. Altovskij; A. Konopljantsev. About the book The manual, is divided into four chapters. The first chapter-general theoretical problems-treats underground water and its behaviour as a natural historical process reflecting certain peculiarities attendant upon […]
Manual For The Systematic Study Of The Regime Of Underground Waters – Altovskij, Konopljantsev — Mir Books

Geotechnical Engineering, Soil Mechanics

Superposition, and Using Point Loads in Place of Distributed Ones

Posted on 1 June 202211 October 2022 by Don Warrington

Elastic solutions to stresses induced by loads at the surface have their limitations, but they allow the use of the principle of superposition. The principle of superposition states that you can add the effects of different loads on a single point. Superposition requires that the stress state of the point be path independent, which is the case with elastic conditions. No matter how you load and unload a point in a system, if elasticity is maintained the result will be the same for a given set of loads.

This is illustrated for point loads in the graphic below, but it applies to distributed loads (such as this and this) as well.

Action of a Number of Concentrated Forces (from Tsytovich (1976))

The stresses that result from each load affect the total stress at the point of interest. They can be computed and added together. So, since “point” loads are physically impossible, is their computation be useful? The answer to this question is “yes” but it takes some judgement, like so many things in geotechnical engineering.

Let us consider the following case, from NAVFAC DM 7.01.

Separate Column Footings Problem, from NAVFAC DM 7.01

To solve this problem, DM 7 converted the square footings to circular ones and then used the chart shown in the post Going Around in Circles for Rigid and Flexible Foundations. This chart, like others, is hard to read. Is it possible to use point loads as a substitute?

There are three things we need to note here. The first is that the load on each column is 27 tons, and is the same for each column.

The second is that the “r” shown in the table above is not the same as it is for the point loads. The variable “r” for the point loads is the horizontal distance from the load to the point of interest.

The third is that there are three column positions shown in the diagram, with three corresponding values of r:

The column on top of the load, Column B2
The columns in the mid-point of the edges, Columns A2, B1, B3, and C2. These have an value of r of 15′ from the point of interest.
The columns in the corners of the square, Columns A1, A3, C1, and C3. These have a value of r of 21.2′ from the point of interest.

Now, instead of the chart, we apply the formula derived earlier for the influence coefficient, which is

$K=3/2\,{\frac {1}{\pi \,\left (1+(\frac{r}{z})^{2}\right )^{5/2}}}$

from which the stress is computed by the equation

$\sigma_z = K \frac{P}{z^2}$

Using this formula, we can construct the table below for this problem.

(1) Z, ft	(2) r/Z for B2	(3) r/Z for A2, B1, B3, C2	(4) r/z for A1, A3, C1, C3	(5) K for B2	(6) K for A2, B1, B3, C2	(7) K for A1, A3, C1, C3	(8) Stress for B2, tsf	(9) Stress for each of A2, B1, B3, C2, tsf	(10) Stress for each of A1, A3, C1, C3, tst	(11) Total Stress at Z, tsf
2	0	7.500	10.607	0.477	0.000	0.000	3.223	0.000	0.000	3.224
4	0	3.750	5.303	0.477	0.001	0.000	0.806	0.001	0.000	0.810
6	0	2.500	3.536	0.477	0.003	0.001	0.358	0.003	0.001	0.370
10	0	1.500	2.121	0.477	0.025	0.007	0.129	0.007	0.002	0.163
15	0	1.000	1.414	0.477	0.084	0.031	0.057	0.010	0.004	0.113
20	0	0.750	1.061	0.477	0.156	0.073	0.032	0.011	0.005	0.094
25	0	0.600	0.849	0.477	0.221	0.123	0.021	0.010	0.005	0.080

Results from the Separate Column Footings problem in NAVFC DM 7.01, using point loads. Note that the stresses in Columns 9 and 10 are due to a single column, and not all four of them. The total stress in column (11) is the sum of the stress in Column 8 plus 4 times the result in Column 9 plus four times the result in Column 10.

We could have opted to add the influence coefficients and then compute the stresses since, for each elevation Z, both the elevation and the column load were the same. We did not for clarity; it is certainly possible to have columns of different loads.

The results are conservative, they tend to be higher, especially at the lower value elevations. It’s worth noting that the total stress at Z = 2′ is higher than the distributed loads on the footings. One way to make this better is to use the center formulae for stress under circles for Column B2 and point loads for the rest. Given, however, the limitations of the method in general, and the considerably lower effort in obtaining the influence coefficients, the method is a reasonable one to use.

Superposition doesn’t only apply to point loads; the example given in the post Analytical Boussinesq Solutions for Strip, Square and Rectangular Loads uses it for rectangular and square loads. Nevertheless, with appropriate engineering judgement, using point loads in place of distributed loads can be a viable option.