Skip to content

ge.design — API reference

Effective stress

total_stress

Python
total_stress(gamma: Union[float, Sequence[float]], depth: Union[float, Sequence[float]]) -> float

Total vertical stress sigma_v = sum(gamma_i * H_i).

PARAMETER DESCRIPTION
gamma

Unit weight(s) of each layer (kN/m^3).

TYPE: float or sequence

depth

Either a single depth (with scalar gamma) -- returns gamma*depth -- or layer thicknesses matching gamma -- returns sum(gamma_i * H_i).

TYPE: float or sequence

RETURNS DESCRIPTION
sigma

Total vertical stress (kPa).

TYPE: float

Reference

Das (2010) Eq. 5.1; Terzaghi (1943).

Source code in geoeq/design/stress.py
Python
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
def total_stress(
    gamma: Union[float, Sequence[float]],
    depth: Union[float, Sequence[float]],
) -> float:
    """Total vertical stress sigma_v = sum(gamma_i * H_i).

    Parameters
    ----------
    gamma : float or sequence
        Unit weight(s) of each layer (kN/m^3).
    depth : float or sequence
        Either a single depth (with scalar gamma) -- returns gamma*depth --
        or layer thicknesses matching ``gamma`` -- returns sum(gamma_i * H_i).

    Returns
    -------
    sigma : float
        Total vertical stress (kPa).

    Reference
    ---------
    Das (2010) Eq. 5.1; Terzaghi (1943).
    """
    g = np.atleast_1d(np.asarray(gamma, dtype=float))
    d = np.atleast_1d(np.asarray(depth, dtype=float))
    check_positive(g, "gamma")
    check_non_negative(d, "depth")
    if g.size == 1 and d.size == 1:
        return float(g[0] * d[0])
    if g.size != d.size:
        raise ValueError(
            "gamma and depth must have equal length (one entry per layer).")
    return float(np.sum(g * d))

pore_pressure

Python
pore_pressure(z: Union[float, Iterable[float]], z_w: float = 0.0, gamma_w: float = GAMMA_WATER) -> Union[float, np.ndarray]

Hydrostatic pore pressure u = gamma_w * (z - z_w), clipped at 0.

PARAMETER DESCRIPTION
z

Depth(s) of interest (m, positive downward).

TYPE: float or array

z_w

Depth of the water table (m). Default 0 (water table at surface).

TYPE: float DEFAULT: 0.0

gamma_w

Unit weight of water (kN/m^3). Default 9.81.

TYPE: float DEFAULT: GAMMA_WATER

RETURNS DESCRIPTION
u

Pore pressure (kPa). Zero above the water table.

TYPE: float or array

Reference

Das (2010) Eq. 5.2.

Source code in geoeq/design/stress.py
Python
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
def pore_pressure(
    z: Union[float, Iterable[float]],
    z_w: float = 0.0,
    gamma_w: float = GAMMA_WATER,
) -> Union[float, np.ndarray]:
    """Hydrostatic pore pressure u = gamma_w * (z - z_w), clipped at 0.

    Parameters
    ----------
    z : float or array
        Depth(s) of interest (m, positive downward).
    z_w : float
        Depth of the water table (m). Default 0 (water table at surface).
    gamma_w : float
        Unit weight of water (kN/m^3). Default 9.81.

    Returns
    -------
    u : float or array
        Pore pressure (kPa). Zero above the water table.

    Reference
    ---------
    Das (2010) Eq. 5.2.
    """
    check_positive(gamma_w, "gamma_w")
    z_arr = np.asarray(z, dtype=float)
    u = gamma_w * np.maximum(0.0, z_arr - z_w)
    return float(u) if u.ndim == 0 else u

effective_stress

Python
effective_stress(sigma: Union[float, ndarray], u: Union[float, ndarray]) -> Union[float, np.ndarray]

Terzaghi effective stress sigma' = sigma - u (kPa).

Reference

Terzaghi (1943); Das (2010) Eq. 5.5.

Source code in geoeq/design/stress.py
Python
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def effective_stress(
    sigma: Union[float, np.ndarray],
    u: Union[float, np.ndarray],
) -> Union[float, np.ndarray]:
    """Terzaghi effective stress sigma' = sigma - u (kPa).

    Reference
    ---------
    Terzaghi (1943); Das (2010) Eq. 5.5.
    """
    sigma = np.asarray(sigma, dtype=float)
    u = np.asarray(u, dtype=float)
    res = sigma - u
    return float(res) if res.ndim == 0 else res

capillary_rise

Python
capillary_rise(D10: float, e: float, C: float = 0.2) -> float

Capillary rise height in soil (Hazen-type estimate).

h_c [cm] = C / (e * D10 [cm])

PARAMETER DESCRIPTION
D10

Effective grain size, the 10% passing diameter (mm).

TYPE: float

e

Void ratio (-).

TYPE: float

C

Empirical constant; 0.1 < C < 0.5 typical. Default 0.2.

TYPE: float DEFAULT: 0.2

RETURNS DESCRIPTION
h_c

Capillary rise (m).

TYPE: float

Reference

Hazen (1930); Das (2010) Eq. 5.8.

Source code in geoeq/design/stress.py
Python
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def capillary_rise(D10: float, e: float, C: float = 0.2) -> float:
    """Capillary rise height in soil (Hazen-type estimate).

    h_c [cm] = C / (e * D10 [cm])

    Parameters
    ----------
    D10 : float
        Effective grain size, the 10% passing diameter (mm).
    e : float
        Void ratio (-).
    C : float
        Empirical constant; 0.1 < C < 0.5 typical. Default 0.2.

    Returns
    -------
    h_c : float
        Capillary rise (m).

    Reference
    ---------
    Hazen (1930); Das (2010) Eq. 5.8.
    """
    check_positive(D10, "D10")
    check_positive(e, "e")
    D10_cm = D10 / 10.0  # mm to cm
    h_c_cm = C / (e * D10_cm)
    return h_c_cm / 100.0  # to m

stress_plot

Python
stress_plot(profile, dz: float = 0.1, **kwargs)

Plot total, pore, and effective stress vs depth.

Delegates to profile.plot(). See SoilProfile.plot.

Source code in geoeq/design/stress.py
Python
133
134
135
136
137
138
def stress_plot(profile, dz: float = 0.1, **kwargs):
    """Plot total, pore, and effective stress vs depth.

    Delegates to ``profile.plot()``. See ``SoilProfile.plot``.
    """
    return profile.plot(dz=dz, **kwargs)

Seepage

darcy_flow

Python
darcy_flow(k: float, i: float, A: float) -> float

Darcy's law: Q = k * i * A (m^3/s).

PARAMETER DESCRIPTION
k

Hydraulic conductivity (m/s).

TYPE: float

i

Hydraulic gradient (-, dimensionless).

TYPE: float

A

Cross-sectional area (m^2).

TYPE: float

Reference

Darcy (1856); Das (2010) Eq. 5.11.

Source code in geoeq/design/seepage.py
Python
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
def darcy_flow(k: float, i: float, A: float) -> float:
    """Darcy's law: Q = k * i * A (m^3/s).

    Parameters
    ----------
    k : float
        Hydraulic conductivity (m/s).
    i : float
        Hydraulic gradient (-, dimensionless).
    A : float
        Cross-sectional area (m^2).

    Reference
    ---------
    Darcy (1856); Das (2010) Eq. 5.11.
    """
    check_positive(k, "k")
    check_positive(A, "A")
    return k * i * A

hydraulic_gradient

Python
hydraulic_gradient(dh: float, L: float) -> float

Hydraulic gradient i = delta h / L.

PARAMETER DESCRIPTION
dh

Head loss (m).

TYPE: float

L

Flow-path length (m).

TYPE: float

Reference

Das (2010) Eq. 5.10.

Source code in geoeq/design/seepage.py
Python
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def hydraulic_gradient(dh: float, L: float) -> float:
    """Hydraulic gradient i = delta h / L.

    Parameters
    ----------
    dh : float
        Head loss (m).
    L : float
        Flow-path length (m).

    Reference
    ---------
    Das (2010) Eq. 5.10.
    """
    check_positive(L, "L")
    return dh / L

critical_gradient

Python
critical_gradient(Gs: float, e: float) -> float

Critical (boiling) hydraulic gradient i_cr = (Gs - 1)/(1 + e).

When i >= i_cr in upward flow, effective stress vanishes and quicksand (heave) develops.

Reference

Das (2010) Eq. 5.18; Terzaghi & Peck (1948).

Source code in geoeq/design/seepage.py
Python
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def critical_gradient(Gs: float, e: float) -> float:
    """Critical (boiling) hydraulic gradient i_cr = (Gs - 1)/(1 + e).

    When i >= i_cr in upward flow, effective stress vanishes and quicksand
    (heave) develops.

    Reference
    ---------
    Das (2010) Eq. 5.18; Terzaghi & Peck (1948).
    """
    check_positive(Gs, "Gs")
    check_positive(e, "e")
    if Gs <= 1:
        raise ValueError("Gs must be > 1.")
    return (Gs - 1) / (1 + e)

equivalent_k

Python
equivalent_k(k_layers: Sequence[float], H_layers: Sequence[float], direction: str = 'horizontal') -> float

Equivalent permeability of a layered system.

horizontal: k_eq = sum(k_i * H_i) / sum(H_i) (parallel flow) vertical: k_eq = sum(H_i) / sum(H_i / k_i) (series flow)

Reference

Das (2010) Eq. 5.14 (horizontal), 5.15 (vertical).

Source code in geoeq/design/seepage.py
Python
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
def equivalent_k(
    k_layers: Sequence[float],
    H_layers: Sequence[float],
    direction: str = "horizontal",
) -> float:
    """Equivalent permeability of a layered system.

    horizontal: k_eq = sum(k_i * H_i) / sum(H_i)        (parallel flow)
    vertical:   k_eq = sum(H_i) / sum(H_i / k_i)        (series flow)

    Reference
    ---------
    Das (2010) Eq. 5.14 (horizontal), 5.15 (vertical).
    """
    k = np.asarray(k_layers, dtype=float)
    H = np.asarray(H_layers, dtype=float)
    if k.shape != H.shape:
        raise ValueError("k_layers and H_layers must have the same length.")
    check_positive(k, "k_layers")
    check_positive(H, "H_layers")
    direction = direction.lower()
    if direction in ("horizontal", "h", "parallel"):
        return float(np.sum(k * H) / np.sum(H))
    elif direction in ("vertical", "v", "series"):
        return float(np.sum(H) / np.sum(H / k))
    raise ValueError("direction must be 'horizontal' or 'vertical'.")

flow_net

Python
flow_net(Nf: float, Nd: float, k: float, dh: float, L: float = 1.0) -> float

Seepage discharge from a flow net.

Q = k * dh * (Nf / Nd) * L (m^3/s per metre length normal to plane)

PARAMETER DESCRIPTION
Nf

Number of flow channels.

TYPE: float

Nd

Number of equipotential drops.

TYPE: float

k

Hydraulic conductivity (m/s).

TYPE: float

dh

Total head loss across the system (m).

TYPE: float

L

Length perpendicular to the flow plane (m). Default 1 m.

TYPE: float DEFAULT: 1.0

Reference

Das (2010) Eq. 5.20.

Source code in geoeq/design/seepage.py
Python
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def flow_net(Nf: float, Nd: float, k: float, dh: float, L: float = 1.0) -> float:
    """Seepage discharge from a flow net.

    Q = k * dh * (Nf / Nd) * L   (m^3/s per metre length normal to plane)

    Parameters
    ----------
    Nf : float
        Number of flow channels.
    Nd : float
        Number of equipotential drops.
    k : float
        Hydraulic conductivity (m/s).
    dh : float
        Total head loss across the system (m).
    L : float
        Length perpendicular to the flow plane (m). Default 1 m.

    Reference
    ---------
    Das (2010) Eq. 5.20.
    """
    check_positive(Nf, "Nf")
    check_positive(Nd, "Nd")
    check_positive(k, "k")
    return k * dh * (Nf / Nd) * L

Stress distribution

boussinesq_point

Python
boussinesq_point(P: float, z: float, r: float = 0.0) -> float

Vertical stress under a point load on a semi-infinite elastic mass.

Text Only
delta_sigma_z = (3 P / (2 pi z^2)) * (1 / (1 + (r/z)^2))^(5/2)
PARAMETER DESCRIPTION
P

Point load (kN).

TYPE: float

z

Depth below ground surface (m).

TYPE: float

r

Radial distance from the line of action (m).

TYPE: float DEFAULT: 0.0

RETURNS DESCRIPTION
delta_sigma

Vertical stress increment (kPa).

TYPE: float

Reference

Boussinesq (1885); Das (2010) Eq. 6.4.

Source code in geoeq/design/boussinesq.py
Python
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def boussinesq_point(P: float, z: float, r: float = 0.0) -> float:
    """Vertical stress under a point load on a semi-infinite elastic mass.

        delta_sigma_z = (3 P / (2 pi z^2)) * (1 / (1 + (r/z)^2))^(5/2)

    Parameters
    ----------
    P : float
        Point load (kN).
    z : float
        Depth below ground surface (m).
    r : float
        Radial distance from the line of action (m).

    Returns
    -------
    delta_sigma : float
        Vertical stress increment (kPa).

    Reference
    ---------
    Boussinesq (1885); Das (2010) Eq. 6.4.
    """
    check_positive(z, "z")
    check_non_negative(r, "r")
    ratio = 1.0 / (1.0 + (r / z) ** 2)
    return float(3 * P / (2 * np.pi * z ** 2) * ratio ** 2.5)

boussinesq_line

Python
boussinesq_line(q: float, z: float, x: float = 0.0) -> float

Vertical stress under an infinite line load q (kN/m).

Text Only
delta_sigma_z = (2 q / pi z) * 1 / (1 + (x/z)^2)^2
Reference

Das (2010) Eq. 6.7.

Source code in geoeq/design/boussinesq.py
Python
66
67
68
69
70
71
72
73
74
75
76
def boussinesq_line(q: float, z: float, x: float = 0.0) -> float:
    """Vertical stress under an infinite line load q (kN/m).

        delta_sigma_z = (2 q / pi z) * 1 / (1 + (x/z)^2)^2

    Reference
    ---------
    Das (2010) Eq. 6.7.
    """
    check_positive(z, "z")
    return float(2 * q / (np.pi * z) * 1 / (1 + (x / z) ** 2) ** 2)

boussinesq_strip

Python
boussinesq_strip(q: float, B: float, z: float, x: float = 0.0) -> float

Vertical stress under a uniformly loaded infinite strip of width B.

Text Only
delta_sigma_z = (q / pi) * [ alpha + sin(alpha) cos(alpha + 2 beta) ]

where alpha and beta are the angles subtended at the point by the edges of the strip (Das Eq. 6.10).

PARAMETER DESCRIPTION
q

Surface pressure (kPa).

TYPE: float

B

Strip width (m).

TYPE: float

z

Depth (m).

TYPE: float

x

Horizontal distance from the centre of the strip (m).

TYPE: float DEFAULT: 0.0

Reference

Das (2010) Eq. 6.10.

Source code in geoeq/design/boussinesq.py
Python
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def boussinesq_strip(q: float, B: float, z: float, x: float = 0.0) -> float:
    """Vertical stress under a uniformly loaded infinite strip of width B.

        delta_sigma_z = (q / pi) * [ alpha + sin(alpha) cos(alpha + 2 beta) ]

    where alpha and beta are the angles subtended at the point by the
    edges of the strip (Das Eq. 6.10).

    Parameters
    ----------
    q : float
        Surface pressure (kPa).
    B : float
        Strip width (m).
    z : float
        Depth (m).
    x : float
        Horizontal distance from the *centre* of the strip (m).

    Reference
    ---------
    Das (2010) Eq. 6.10.
    """
    check_positive(z, "z")
    check_positive(B, "B")
    # Geometry: edges at x - B/2 and x + B/2, both at depth z.
    x1 = x - B / 2.0
    x2 = x + B / 2.0
    alpha = np.arctan(x2 / z) - np.arctan(x1 / z)
    # beta is the angle to the closer edge from the vertical (Das defn).
    beta = np.arctan(x1 / z)
    return float(q / np.pi * (alpha + np.sin(alpha) * np.cos(alpha + 2 * beta)))

boussinesq_circular

Python
boussinesq_circular(q: float, R: float, z: float) -> float

Vertical stress on the centreline beneath a uniformly loaded circle.

Text Only
delta_sigma_z = q * [ 1 - 1 / (1 + (R/z)^2)^(3/2) ]
Reference

Das (2010) Eq. 6.15.

Source code in geoeq/design/boussinesq.py
Python
119
120
121
122
123
124
125
126
127
128
129
130
def boussinesq_circular(q: float, R: float, z: float) -> float:
    """Vertical stress on the centreline beneath a uniformly loaded circle.

        delta_sigma_z = q * [ 1 - 1 / (1 + (R/z)^2)^(3/2) ]

    Reference
    ---------
    Das (2010) Eq. 6.15.
    """
    check_positive(z, "z")
    check_positive(R, "R")
    return float(q * (1 - 1 / (1 + (R / z) ** 2) ** 1.5))

boussinesq_rect

Python
boussinesq_rect(q: float, B: float, L: float, z: float, position: str = 'corner') -> float

Vertical stress under a uniformly loaded rectangle B x L.

Uses Fadum (1948) corner influence value, with superposition for 'centre' (4x mB/2 x L/2 sub-rectangles) and 'edge_mid' positions.

PARAMETER DESCRIPTION
q

Surface pressure (kPa).

TYPE: float

B

Rectangle dimensions (m). Convention: L >= B.

TYPE: float

L

Rectangle dimensions (m). Convention: L >= B.

TYPE: float

z

Depth below surface (m).

TYPE: float

position

'corner' -- under one corner of the rectangle 'centre' -- under the centre (sum of 4 sub-rectangles) 'edge' -- under the midpoint of a long edge (2 sub-rectangles)

TYPE: str DEFAULT: 'corner'

Reference

Fadum (1948); Das (2010) Eq. 6.29-6.31.

Source code in geoeq/design/boussinesq.py
Python
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
def boussinesq_rect(
    q: float, B: float, L: float, z: float, position: str = "corner",
) -> float:
    """Vertical stress under a uniformly loaded rectangle B x L.

    Uses Fadum (1948) corner influence value, with superposition for
    'centre' (4x mB/2 x L/2 sub-rectangles) and 'edge_mid' positions.

    Parameters
    ----------
    q : float
        Surface pressure (kPa).
    B, L : float
        Rectangle dimensions (m). Convention: L >= B.
    z : float
        Depth below surface (m).
    position : str
        'corner'  -- under one corner of the rectangle
        'centre'  -- under the centre (sum of 4 sub-rectangles)
        'edge'    -- under the midpoint of a long edge (2 sub-rectangles)

    Reference
    ---------
    Fadum (1948); Das (2010) Eq. 6.29-6.31.
    """
    check_positive(z, "z")
    check_positive(B, "B")
    check_positive(L, "L")
    position = position.lower()

    if position == "corner":
        I = newmark_influence(B / z, L / z)
        return float(q * I)
    elif position in ("centre", "center"):
        # Four rectangles each B/2 x L/2 sharing a corner at the centre.
        I = newmark_influence((B / 2) / z, (L / 2) / z)
        return float(q * 4 * I)
    elif position == "edge":
        # Two rectangles each B/2 x L sharing a corner at the edge mid.
        I = newmark_influence((B / 2) / z, L / z)
        return float(q * 2 * I)
    raise ValueError("position must be 'corner', 'centre', or 'edge'.")

westergaard_point

Python
westergaard_point(P: float, z: float, r: float = 0.0, mu: float = 0.0) -> float

Westergaard vertical stress for a point load on a layered medium.

Text Only
delta_sigma_z = (P / (pi z^2)) * eta / (eta^2 + (r/z)^2)^(3/2)

where eta = sqrt( (1 - 2mu) / (2 - 2mu) ).

For thinly stratified soils (alternating stiff/soft) the Westergaard solution gives lower stresses than Boussinesq.

Reference

Westergaard (1938); Das (2010) Eq. 6.13.

Source code in geoeq/design/boussinesq.py
Python
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
def westergaard_point(P: float, z: float, r: float = 0.0, mu: float = 0.0) -> float:
    """Westergaard vertical stress for a point load on a layered medium.

        delta_sigma_z = (P / (pi z^2)) * eta / (eta^2 + (r/z)^2)^(3/2)

    where eta = sqrt( (1 - 2mu) / (2 - 2mu) ).

    For thinly stratified soils (alternating stiff/soft) the Westergaard
    solution gives lower stresses than Boussinesq.

    Reference
    ---------
    Westergaard (1938); Das (2010) Eq. 6.13.
    """
    check_positive(z, "z")
    check_non_negative(r, "r")
    if not 0 <= mu < 0.5:
        raise ValueError("mu must be in [0, 0.5).")
    eta_sq = (1 - 2 * mu) / (2 - 2 * mu) if mu < 0.5 else 0.0
    eta = np.sqrt(eta_sq) if eta_sq > 0 else 1e-9
    return float(P / (np.pi * z ** 2) * eta / (eta_sq + (r / z) ** 2) ** 1.5)

newmark_influence

Python
newmark_influence(m: float, n: float) -> float

Fadum / Newmark influence value I for a rectangle B x L at depth z.

Text Only
m = B / z,   n = L / z
RETURNS DESCRIPTION
I

Dimensionless influence factor (Das Eq. 6.30). The stress at the corner of a uniformly loaded B x L rectangle is q * I.

TYPE: float

Reference

Fadum (1948); Newmark (1935); Das Eq. 6.30, Table 6.5.

Source code in geoeq/design/boussinesq.py
Python
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
def newmark_influence(m: float, n: float) -> float:
    """Fadum / Newmark influence value I for a rectangle B x L at depth z.

        m = B / z,   n = L / z

    Returns
    -------
    I : float
        Dimensionless influence factor (Das Eq. 6.30). The stress at the
        corner of a uniformly loaded B x L rectangle is q * I.

    Reference
    ---------
    Fadum (1948); Newmark (1935); Das Eq. 6.30, Table 6.5.
    """
    check_positive(m, "m")
    check_positive(n, "n")
    # Fadum's closed-form.
    s = m * m + n * n + 1.0
    term1 = (2 * m * n * np.sqrt(s)) / (s + m * m * n * n) * (s + 1) / s
    # Newmark formula (Fadum 1948):
    # I = (1/4 pi) * { 2 m n sqrt(m^2+n^2+1) / (m^2+n^2+1 + m^2 n^2)
    #                  * (m^2+n^2+2)/(m^2+n^2+1)
    #                 + atan( 2 m n sqrt(m^2+n^2+1) / (m^2+n^2+1 - m^2 n^2) ) }
    a = 2 * m * n * np.sqrt(s)
    b = s + m * m * n * n
    c = (s + 1) / s
    d = s - m * m * n * n
    atan_arg = a / d
    # When d < 0, the angle is in the second quadrant and we add pi.
    if d > 0:
        atan_term = np.arctan(atan_arg)
    else:
        atan_term = np.arctan(atan_arg) + np.pi
    I = (1.0 / (4 * np.pi)) * (a / b * c + atan_term)
    return float(I)

stress_2to1

Python
stress_2to1(q: float, B: float, L: float, z: float) -> float

2:1 vertical-spread approximation for a B x L footing.

Text Only
delta_sigma_z = q * B * L / ((B + z) * (L + z))

Crude but engineering-useful for spread footings up to z ~ 2B.

Reference

Das (2010) Eq. 6.34.

Source code in geoeq/design/boussinesq.py
Python
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
def stress_2to1(q: float, B: float, L: float, z: float) -> float:
    """2:1 vertical-spread approximation for a B x L footing.

        delta_sigma_z = q * B * L / ((B + z) * (L + z))

    Crude but engineering-useful for spread footings up to z ~ 2B.

    Reference
    ---------
    Das (2010) Eq. 6.34.
    """
    check_positive(z, "z")
    check_positive(B, "B")
    check_positive(L, "L")
    return float(q * B * L / ((B + z) * (L + z)))

stress_bulb

Python
stress_bulb(P: float, z_range=None, r_range=None, n: int = 50)

Compute a stress-bulb mesh under a point load.

RETURNS DESCRIPTION
dict with keys 'z', 'r', 'sigma' -- 2D grids suitable for contourf.
Reference

Das (2010) Fig. 6.6 (qualitative isobar pattern).

Source code in geoeq/design/boussinesq.py
Python
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
def stress_bulb(P: float, z_range=None, r_range=None, n: int = 50):
    """Compute a stress-bulb mesh under a point load.

    Returns
    -------
    dict with keys 'z', 'r', 'sigma' -- 2D grids suitable for contourf.

    Reference
    ---------
    Das (2010) Fig. 6.6 (qualitative isobar pattern).
    """
    if z_range is None:
        z_range = (0.1, 10.0)
    if r_range is None:
        r_range = (-5.0, 5.0)
    z = np.linspace(z_range[0], z_range[1], n)
    r = np.linspace(r_range[0], r_range[1], n)
    R, Z = np.meshgrid(r, z)
    ratio = 1.0 / (1.0 + (R / Z) ** 2)
    sigma = 3 * P / (2 * np.pi * Z ** 2) * ratio ** 2.5
    return {"z": Z, "r": R, "sigma": sigma}

stress_isobar_plot

Python
stress_isobar_plot(P: float = 100.0, z_max: float = 8.0, r_max: float = 4.0, levels=None, ax=None, save_as: str = None)

Stress isobars (pressure bulb) beneath a point load.

PARAMETER DESCRIPTION
P

Point load (kN).

TYPE: float DEFAULT: 100.0

z_max

Plot extents (m).

TYPE: float DEFAULT: 8.0

r_max

Plot extents (m).

TYPE: float DEFAULT: 8.0

levels

Iso-stress contour levels (kPa).

TYPE: sequence DEFAULT: None

Source code in geoeq/design/plots.py
Python
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
def stress_isobar_plot(P: float = 100.0, z_max: float = 8.0,
                       r_max: float = 4.0, levels=None,
                       ax=None, save_as: str = None):
    """Stress isobars (pressure bulb) beneath a point load.

    Parameters
    ----------
    P : float
        Point load (kN).
    z_max, r_max : float
        Plot extents (m).
    levels : sequence
        Iso-stress contour levels (kPa).
    """
    import matplotlib.pyplot as plt
    res = stress_bulb(P=P, z_range=(0.05, z_max),
                       r_range=(-r_max, r_max), n=120)
    if levels is None:
        # Geometric spacing from max stress in the grid.
        smax = np.nanmax(res["sigma"])
        levels = np.geomspace(max(smax * 0.01, 0.1), smax * 0.9, 8)
    if ax is None:
        fig, ax = plt.subplots(figsize=(7, 7))
    else:
        fig = ax.figure
    cs = ax.contourf(res["r"], res["z"], res["sigma"],
                     levels=levels, cmap="RdYlBu_r", alpha=0.8)
    ax.contour(res["r"], res["z"], res["sigma"],
               levels=levels, colors="black", linewidths=0.5)
    cbar = fig.colorbar(cs, ax=ax, label=r"$\Delta\sigma_z$ (kPa)")
    # Mark the load application point.
    ax.plot([0], [0], "v", color="black", markersize=14)
    ax.text(0, -0.2, f"  $P = {P:.0f}$ kN", ha="center", va="bottom",
            fontsize=10)
    ax.invert_yaxis()
    ax.set_xlabel(r"Radial distance $r$ (m)")
    ax.set_ylabel(r"Depth $z$ (m)")
    ax.set_title("Boussinesq pressure bulb beneath a point load")
    ax.set_aspect("equal")
    if save_as:
        fig.savefig(save_as, dpi=300, bbox_inches="tight")
    return fig

Bearing capacity

bearing_factors

Python
bearing_factors(phi: float, method: str = 'meyerhof') -> dict

Bearing-capacity factors Nc, Nq, N_gamma for a given friction angle.

PARAMETER DESCRIPTION
phi

Effective friction angle (degrees).

TYPE: float

method

'terzaghi' | 'meyerhof' | 'hansen' | 'vesic'.

TYPE: str DEFAULT: 'meyerhof'

RETURNS DESCRIPTION
dict

{'Nc': ..., 'Nq': ..., 'Ngamma': ...}.
Reference

Das (2014) Tables 4.1-4.3.

Source code in geoeq/design/bearing.py
Python
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
def bearing_factors(phi: float, method: str = "meyerhof") -> dict:
    """Bearing-capacity factors Nc, Nq, N_gamma for a given friction angle.

    Parameters
    ----------
    phi : float
        Effective friction angle (degrees).
    method : str
        ``'terzaghi'`` | ``'meyerhof'`` | ``'hansen'`` | ``'vesic'``.

    Returns
    -------
    dict
        ``{'Nc': ..., 'Nq': ..., 'Ngamma': ...}``.

    Reference
    ---------
    Das (2014) Tables 4.1-4.3.
    """
    check_non_negative(phi, "phi")
    check_range(phi, "phi", 0, 50)
    method = method.lower()
    phi_rad = np.radians(phi)

    # Nq is identical across Meyerhof/Hansen/Vesic (Prandtl-Reissner).
    if phi == 0:
        Nq = 1.0
    else:
        Nq = np.exp(np.pi * np.tan(phi_rad)) * \
            np.tan(np.pi / 4 + phi_rad / 2) ** 2

    # Nc via classical relation Nc = (Nq - 1) / tan(phi); Nc = 5.14 at phi=0.
    if phi == 0:
        Nc = 5.14  # = pi + 2
    else:
        Nc = (Nq - 1) / np.tan(phi_rad)

    if method == "terzaghi":
        # Terzaghi's original Nq, Nc differ slightly (he used local shear);
        # use the Bowles-tabulated forms for general shear.
        Nq = np.exp((1.5 * np.pi - phi_rad) * np.tan(phi_rad)) / \
            (2 * np.cos(np.pi / 4 + phi_rad / 2) ** 2)
        if phi == 0:
            Nc = 5.7  # Terzaghi original
        else:
            Nc = (Nq - 1) / np.tan(phi_rad)
        # Terzaghi N_gamma (empirical fit to chart, Bowles 1996):
        Ngamma = (Nq - 1) * np.tan(1.4 * phi_rad) if phi > 0 else 0.0
    elif method == "meyerhof":
        Ngamma = (Nq - 1) * np.tan(1.4 * phi_rad)
    elif method == "hansen":
        Ngamma = 1.5 * (Nq - 1) * np.tan(phi_rad)
    elif method == "vesic":
        Ngamma = 2 * (Nq + 1) * np.tan(phi_rad)
    else:
        raise ValueError(
            "method must be 'terzaghi', 'meyerhof', 'hansen', or 'vesic'.")
    return {"Nc": float(Nc), "Nq": float(Nq), "Ngamma": float(Ngamma)}

bearing_capacity

Python
bearing_capacity(c: float, gamma: float, Df: float, B: float, phi: float, L: float = None, method: str = 'meyerhof', shape: bool = True, depth: bool = True, inclination: float = 0.0, gamma_above: float = None, water_table: str = 'deep') -> dict

Ultimate bearing capacity q_u of a shallow footing.

PARAMETER DESCRIPTION
c

Effective cohesion (kPa). Use Su for undrained (phi = 0) analysis.

TYPE: float

gamma

Effective unit weight of soil below the footing (kN/m^3). Use gamma' = gamma_sat - gamma_w if the water table is at the base of the footing.

TYPE: float

Df

Depth of embedment (m).

TYPE: float

B

Footing width (m).

TYPE: float

phi

Effective friction angle (degrees).

TYPE: float

L

Footing length (m). None -> strip footing (B/L -> 0).

TYPE: float DEFAULT: None

method

'terzaghi' | 'meyerhof' | 'hansen' | 'vesic'.

TYPE: str DEFAULT: 'meyerhof'

shape

Whether to apply shape/depth corrections. Terzaghi ignores both.

TYPE: bool DEFAULT: True

depth

Whether to apply shape/depth corrections. Terzaghi ignores both.

TYPE: bool DEFAULT: True

inclination

Load inclination from vertical (deg). Default 0.

TYPE: float DEFAULT: 0.0

gamma_above

Effective unit weight above the footing (kN/m^3). Default = gamma.

TYPE: float DEFAULT: None

water_table

'deep' (no correction) | 'at_base' (gamma below = gamma') | 'above_base' (full submergence -- conservative).

TYPE: str DEFAULT: 'deep'

RETURNS DESCRIPTION
dict with 'q_u' (kPa), 'Nc', 'Nq', 'Ngamma', and the applied factors.
Reference

Das (2014), Ch. 4.

Source code in geoeq/design/bearing.py
Python
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
def bearing_capacity(
    c: float, gamma: float, Df: float, B: float, phi: float,
    L: float = None, method: str = "meyerhof",
    shape: bool = True, depth: bool = True,
    inclination: float = 0.0, gamma_above: float = None,
    water_table: str = "deep",
) -> dict:
    """Ultimate bearing capacity q_u of a shallow footing.

    Parameters
    ----------
    c : float
        Effective cohesion (kPa). Use Su for undrained (phi = 0) analysis.
    gamma : float
        Effective unit weight of soil *below* the footing (kN/m^3).
        Use gamma' = gamma_sat - gamma_w if the water table is at the
        base of the footing.
    Df : float
        Depth of embedment (m).
    B : float
        Footing width (m).
    phi : float
        Effective friction angle (degrees).
    L : float, optional
        Footing length (m). ``None`` -> strip footing (B/L -> 0).
    method : str
        'terzaghi' | 'meyerhof' | 'hansen' | 'vesic'.
    shape, depth : bool
        Whether to apply shape/depth corrections. Terzaghi ignores both.
    inclination : float
        Load inclination from vertical (deg). Default 0.
    gamma_above : float
        Effective unit weight above the footing (kN/m^3). Default = gamma.
    water_table : str
        'deep' (no correction) | 'at_base' (gamma below = gamma') |
        'above_base' (full submergence -- conservative).

    Returns
    -------
    dict with 'q_u' (kPa), 'Nc', 'Nq', 'Ngamma', and the applied factors.

    Reference
    ---------
    Das (2014), Ch. 4.
    """
    check_positive(B, "B")
    check_non_negative(Df, "Df")
    check_non_negative(c, "c")
    check_positive(gamma, "gamma")
    if gamma_above is None:
        gamma_above = gamma

    factors = bearing_factors(phi, method=method)
    Nc, Nq, Ng = factors["Nc"], factors["Nq"], factors["Ngamma"]

    if method == "terzaghi":
        # Apply Terzaghi's shape multipliers for square/circle (Das Eq. 4.6-4.7).
        if L is not None and np.isclose(L, B):
            # Square
            sc_, sg_ = 1.3, 0.8
        elif L is None or L > 10 * B:
            sc_, sg_ = 1.0, 1.0
        else:
            sc_, sg_ = 1.3, 0.8  # treat near-square as square; conservative
        q = gamma_above * Df
        q_u = c * Nc * sc_ + q * Nq + 0.5 * gamma * B * Ng * sg_
        return {
            "q_u": float(q_u), "Nc": Nc, "Nq": Nq, "Ngamma": Ng,
            "method": method,
        }

    # Strip default
    L_eff = L if L is not None else 1e6 * B

    sc = sq = sg = 1.0
    dc = dq = dg = 1.0
    ic = iq = ig = 1.0
    if shape:
        sf = bearing_shape_factors(B, L_eff, phi, method=method)
        sc, sq, sg = sf["sc"], sf["sq"], sf["s_gamma"]
    if depth:
        df_ = bearing_depth_factors(Df, B, phi, method=method)
        dc, dq, dg = df_["dc"], df_["dq"], df_["d_gamma"]
    if inclination > 0:
        i_ = bearing_inclination_factors(inclination, phi, c, method=method)
        ic, iq, ig = i_["ic"], i_["iq"], i_["i_gamma"]

    q = gamma_above * Df
    q_u = (c * Nc * sc * dc * ic
           + q * Nq * sq * dq * iq
           + 0.5 * gamma * B * Ng * sg * dg * ig)
    return {
        "q_u": float(q_u),
        "Nc": Nc, "Nq": Nq, "Ngamma": Ng,
        "sc": sc, "sq": sq, "s_gamma": sg,
        "dc": dc, "dq": dq, "d_gamma": dg,
        "method": method,
    }

bearing_allowable

Python
bearing_allowable(q_u: float, FS: float = 3.0) -> float

Allowable bearing pressure q_all = q_u / FS.

Standard FS=3 for shallow foundations under static loads (Das 2014, p. 224).

Source code in geoeq/design/bearing.py
Python
295
296
297
298
299
300
301
302
def bearing_allowable(q_u: float, FS: float = 3.0) -> float:
    """Allowable bearing pressure q_all = q_u / FS.

    Standard FS=3 for shallow foundations under static loads
    (Das 2014, p. 224).
    """
    check_positive(FS, "FS")
    return q_u / FS

bearing_shape_factors

Python
bearing_shape_factors(B: float, L: float, phi: float, method: str = 'meyerhof') -> dict

Shape factors sc, sq, s_gamma for B x L footing.

Reference

Das (2014) Table 4.3.

Source code in geoeq/design/bearing.py
Python
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def bearing_shape_factors(B: float, L: float, phi: float,
                          method: str = "meyerhof") -> dict:
    """Shape factors sc, sq, s_gamma for B x L footing.

    Reference
    ---------
    Das (2014) Table 4.3.
    """
    check_positive(B, "B")
    check_positive(L, "L")
    if B > L:
        B, L = L, B  # ensure B is the shorter side
    phi_rad = np.radians(phi)
    method = method.lower()
    Kp = np.tan(np.pi / 4 + phi_rad / 2) ** 2
    if method == "meyerhof":
        sc = 1 + 0.2 * Kp * B / L
        if phi >= 10:
            sq = s_gamma = 1 + 0.1 * Kp * B / L
        else:
            sq = s_gamma = 1.0
    elif method in ("hansen", "vesic"):
        Nq = np.exp(np.pi * np.tan(phi_rad)) * \
            np.tan(np.pi / 4 + phi_rad / 2) ** 2 if phi > 0 else 1.0
        Nc = (Nq - 1) / np.tan(phi_rad) if phi > 0 else 5.14
        sc = 1 + (Nq / Nc) * (B / L)
        sq = 1 + (B / L) * np.tan(phi_rad)
        s_gamma = max(1 - 0.4 * (B / L), 0.6)
    else:
        raise ValueError("method must be 'meyerhof', 'hansen', or 'vesic'.")
    return {"sc": float(sc), "sq": float(sq), "s_gamma": float(s_gamma)}

bearing_depth_factors

Python
bearing_depth_factors(Df: float, B: float, phi: float, method: str = 'meyerhof') -> dict

Depth factors dc, dq, d_gamma.

Reference

Das (2014) Table 4.3.

Source code in geoeq/design/bearing.py
Python
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
def bearing_depth_factors(Df: float, B: float, phi: float,
                          method: str = "meyerhof") -> dict:
    """Depth factors dc, dq, d_gamma.

    Reference
    ---------
    Das (2014) Table 4.3.
    """
    check_positive(B, "B")
    check_non_negative(Df, "Df")
    phi_rad = np.radians(phi)
    method = method.lower()
    Kp = np.tan(np.pi / 4 + phi_rad / 2) ** 2
    if method == "meyerhof":
        dc = 1 + 0.2 * np.sqrt(Kp) * Df / B
        if phi >= 10:
            dq = d_gamma = 1 + 0.1 * np.sqrt(Kp) * Df / B
        else:
            dq = d_gamma = 1.0
    elif method in ("hansen", "vesic"):
        # k factor (Hansen): Df/B if <= 1; atan(Df/B) [rad] if > 1.
        k = Df / B if Df / B <= 1 else np.arctan(Df / B)
        if phi == 0:
            dc = 1 + 0.4 * k
        else:
            dc = 1 + 0.4 * k
        dq = 1 + 2 * np.tan(phi_rad) * (1 - np.sin(phi_rad)) ** 2 * k
        d_gamma = 1.0
    else:
        raise ValueError("method must be 'meyerhof', 'hansen', or 'vesic'.")
    return {"dc": float(dc), "dq": float(dq), "d_gamma": float(d_gamma)}

bearing_inclination_factors

Python
bearing_inclination_factors(beta: float, phi: float, c: float = 0.0, method: str = 'meyerhof') -> dict

Inclination factors ic, iq, i_gamma for load inclined beta degrees from vertical.

Reference

Das (2014) Table 4.3; Hansen (1970).

Source code in geoeq/design/bearing.py
Python
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
def bearing_inclination_factors(beta: float, phi: float, c: float = 0.0,
                                method: str = "meyerhof") -> dict:
    """Inclination factors ic, iq, i_gamma for load inclined beta degrees
    from vertical.

    Reference
    ---------
    Das (2014) Table 4.3; Hansen (1970).
    """
    check_non_negative(beta, "beta")
    check_range(beta, "beta", 0, 90)
    method = method.lower()
    if method == "meyerhof":
        ic = iq = (1 - beta / 90) ** 2
        i_gamma = (1 - beta / phi) ** 2 if phi > 0 else 1.0
    else:
        ic = iq = (1 - beta / 90) ** 2
        i_gamma = (1 - beta / phi) ** 2 if phi > 0 else 1.0
    return {"ic": float(ic), "iq": float(iq), "i_gamma": float(i_gamma)}

bearing_capacity_plot

Python
bearing_capacity_plot(c: float, gamma: float, Df: float, phi: float, B_range: Sequence[float] = None, methods: Sequence[str] = None, L_over_B: float = 1.0, ax=None, save_as: str = None)

Plot ultimate bearing capacity vs. footing width B.

Useful design-chart helper: shows how q_u scales with B and how the four classical methods diverge for granular soils.

PARAMETER DESCRIPTION
c

Soil parameters (kPa, kN/m^3, m, degrees).

TYPE: float

gamma

Soil parameters (kPa, kN/m^3, m, degrees).

TYPE: float

Df

Soil parameters (kPa, kN/m^3, m, degrees).

TYPE: float

phi

Soil parameters (kPa, kN/m^3, m, degrees).

TYPE: float

B_range

Footing widths to evaluate (m). Default 0.5 .. 5 m.

TYPE: sequence DEFAULT: None

methods

Subset of {"terzaghi", "meyerhof", "hansen", "vesic"}. Default: all four.

TYPE: sequence of str DEFAULT: None

L_over_B

Footing aspect ratio (1.0 = square). Use 1e6 for strip.

TYPE: float DEFAULT: 1.0

Source code in geoeq/design/plots.py
Python
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def bearing_capacity_plot(
    c: float, gamma: float, Df: float, phi: float,
    B_range: Sequence[float] = None, methods: Sequence[str] = None,
    L_over_B: float = 1.0, ax=None, save_as: str = None,
):
    """Plot ultimate bearing capacity vs. footing width B.

    Useful design-chart helper: shows how q_u scales with B and how
    the four classical methods diverge for granular soils.

    Parameters
    ----------
    c, gamma, Df, phi : float
        Soil parameters (kPa, kN/m^3, m, degrees).
    B_range : sequence
        Footing widths to evaluate (m). Default 0.5 .. 5 m.
    methods : sequence of str
        Subset of ``{"terzaghi", "meyerhof", "hansen", "vesic"}``.
        Default: all four.
    L_over_B : float
        Footing aspect ratio (1.0 = square). Use 1e6 for strip.
    """
    import matplotlib.pyplot as plt
    if B_range is None:
        B_range = np.linspace(0.5, 5.0, 30)
    if methods is None:
        methods = ["terzaghi", "meyerhof", "hansen", "vesic"]
    if ax is None:
        fig, ax = plt.subplots(figsize=(7, 5))
    else:
        fig = ax.figure
    colours = {"terzaghi": "#888888", "meyerhof": "#1f3a93",
               "hansen": "#27ae60", "vesic": "#c0392b"}
    for m in methods:
        q = []
        for B in B_range:
            L = B * L_over_B if np.isfinite(L_over_B) else None
            res = bearing_capacity(c, gamma, Df, B, phi, L=L, method=m)
            q.append(res["q_u"])
        ax.plot(B_range, q, label=m.title(),
                color=colours.get(m, "k"), linewidth=1.8)
    ax.set_xlabel("Footing width $B$ (m)")
    ax.set_ylabel("Ultimate bearing capacity $q_u$ (kPa)")
    ax.set_title(
        f"$q_u$ vs $B$ — $c={c}$, $\\phi={phi}°$, $D_f={Df}$ m, $\\gamma={gamma}$ kN/m³")
    ax.grid(alpha=0.3)
    ax.legend()
    if save_as:
        fig.savefig(save_as, dpi=300, bbox_inches="tight")
    return fig

Settlement

settlement_immediate

Python
settlement_immediate(q: float, B: float, Es: float, mu: float = 0.3, Ip: float = None, shape: str = 'rigid_square') -> float

Elastic (immediate) settlement under a uniformly loaded footing.

Text Only
S_i = q * B * (1 - mu^2) / Es * Ip
PARAMETER DESCRIPTION
q

Surface pressure (kPa).

TYPE: float

B

Footing width (m) -- shorter dimension for rectangles.

TYPE: float

Es

Soil elastic modulus (kPa).

TYPE: float

mu

Poisson's ratio (-). Default 0.3.

TYPE: float DEFAULT: 0.3

Ip

Influence factor. If None, picked from shape.

TYPE: float DEFAULT: None

shape

'rigid_square' -> Ip = 0.88 'rigid_circle' -> Ip = 0.79 'flexible_centre' -> Ip = 1.12 'flexible_corner' -> Ip = 0.56 'rigid_strip' -> Ip = 2.0 (approx, depth-dependent)

TYPE: str DEFAULT: 'rigid_square'

Reference

Janbu et al. (1956); Das (2014) Eq. 5.20, Table 5.1.

Source code in geoeq/design/settlement.py
Python
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
def settlement_immediate(
    q: float, B: float, Es: float, mu: float = 0.3,
    Ip: float = None, shape: str = "rigid_square",
) -> float:
    """Elastic (immediate) settlement under a uniformly loaded footing.

        S_i = q * B * (1 - mu^2) / Es * Ip

    Parameters
    ----------
    q : float
        Surface pressure (kPa).
    B : float
        Footing width (m) -- shorter dimension for rectangles.
    Es : float
        Soil elastic modulus (kPa).
    mu : float
        Poisson's ratio (-). Default 0.3.
    Ip : float, optional
        Influence factor. If None, picked from ``shape``.
    shape : str
        'rigid_square'  -> Ip = 0.88
        'rigid_circle'  -> Ip = 0.79
        'flexible_centre' -> Ip = 1.12
        'flexible_corner' -> Ip = 0.56
        'rigid_strip'   -> Ip = 2.0 (approx, depth-dependent)

    Reference
    ---------
    Janbu et al. (1956); Das (2014) Eq. 5.20, Table 5.1.
    """
    check_positive(B, "B")
    check_positive(Es, "Es")
    if Ip is None:
        table = {
            "rigid_square": 0.88, "rigid_circle": 0.79,
            "flexible_centre": 1.12, "flexible_center": 1.12,
            "flexible_corner": 0.56, "rigid_strip": 2.0,
        }
        if shape not in table:
            raise ValueError(
                f"shape must be one of {list(table.keys())}, got {shape!r}.")
        Ip = table[shape]
    return float(q * B * (1 - mu ** 2) / Es * Ip)

settlement_primary

Python
settlement_primary(Cc: float, e0: float, H: float, sigma0: float, delta_sigma: float, Cr: float = None, sigma_pc: float = None, kind: str = 'auto') -> float

1-D primary consolidation settlement.

Three regimes:

  1. NC (normally consolidated): S = (Cc * H / (1 + e0)) * log10((sigma0 + delta_sigma) / sigma0)

  2. OC, final stress below preconsolidation: S = (Cr * H / (1 + e0)) * log10((sigma0 + delta_sigma) / sigma0)

  3. OC, final stress crosses preconsolidation: S = (Cr*H/(1+e0)) * log10(sigma_pc / sigma0)

    • (Cc*H/(1+e0)) * log10((sigma0+delta_sigma)/sigma_pc)
PARAMETER DESCRIPTION
Cc

Compression index (-).

TYPE: float

e0

Initial void ratio (-).

TYPE: float

H

Drainage path length (m) -- full layer thickness for one-way drainage or half-thickness for two-way drainage; pass actual layer thickness.

TYPE: float

sigma0

Initial effective vertical stress (kPa) at mid-depth.

TYPE: float

delta_sigma

Stress increase at mid-depth (kPa).

TYPE: float

Cr

Recompression (swell) index (-). Required for OC cases.

TYPE: float DEFAULT: None

sigma_pc

Preconsolidation pressure (kPa). Required for OC cases.

TYPE: float DEFAULT: None

kind

'NC', 'OC', or 'auto' (decide from sigma_pc).

TYPE: str DEFAULT: 'auto'

RETURNS DESCRIPTION
S

Primary settlement (m).

TYPE: float

Reference

Terzaghi (1925); Das (2014) Eq. 5.31-5.34.

Source code in geoeq/design/settlement.py
Python
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
def settlement_primary(
    Cc: float, e0: float, H: float, sigma0: float, delta_sigma: float,
    Cr: float = None, sigma_pc: float = None, kind: str = "auto",
) -> float:
    """1-D primary consolidation settlement.

    Three regimes:

    1. **NC** (normally consolidated):
       S = (Cc * H / (1 + e0)) * log10((sigma0 + delta_sigma) / sigma0)

    2. **OC**, final stress below preconsolidation:
       S = (Cr * H / (1 + e0)) * log10((sigma0 + delta_sigma) / sigma0)

    3. **OC**, final stress crosses preconsolidation:
       S = (Cr*H/(1+e0)) * log10(sigma_pc / sigma0)
         + (Cc*H/(1+e0)) * log10((sigma0+delta_sigma)/sigma_pc)

    Parameters
    ----------
    Cc : float
        Compression index (-).
    e0 : float
        Initial void ratio (-).
    H : float
        Drainage path length (m) -- full layer thickness for one-way drainage
        or half-thickness for two-way drainage; pass actual layer thickness.
    sigma0 : float
        Initial effective vertical stress (kPa) at mid-depth.
    delta_sigma : float
        Stress increase at mid-depth (kPa).
    Cr : float
        Recompression (swell) index (-). Required for OC cases.
    sigma_pc : float
        Preconsolidation pressure (kPa). Required for OC cases.
    kind : str
        'NC', 'OC', or 'auto' (decide from sigma_pc).

    Returns
    -------
    S : float
        Primary settlement (m).

    Reference
    ---------
    Terzaghi (1925); Das (2014) Eq. 5.31-5.34.
    """
    check_positive(H, "H")
    check_positive(sigma0, "sigma0")
    check_positive(delta_sigma, "delta_sigma")
    if e0 <= 0:
        raise ValueError("e0 must be positive.")

    sigma_f = sigma0 + delta_sigma

    if kind.lower() == "auto":
        if sigma_pc is None or sigma_pc <= sigma0:
            kind = "NC"
        elif sigma_f <= sigma_pc:
            kind = "OC_below"
        else:
            kind = "OC_cross"
    else:
        kind = kind.upper()

    if kind == "NC":
        return float(Cc * H / (1 + e0) * np.log10(sigma_f / sigma0))

    if Cr is None:
        raise ValueError("Cr is required for an OC analysis.")
    if sigma_pc is None:
        raise ValueError("sigma_pc is required for an OC analysis.")

    if kind in ("OC_below", "OC"):
        return float(Cr * H / (1 + e0) * np.log10(sigma_f / sigma0))

    if kind == "OC_cross":
        S1 = Cr * H / (1 + e0) * np.log10(sigma_pc / sigma0)
        S2 = Cc * H / (1 + e0) * np.log10(sigma_f / sigma_pc)
        return float(S1 + S2)

    raise ValueError(
        f"kind must be 'NC', 'OC', 'OC_below', 'OC_cross', or 'auto', "
        f"got {kind!r}.")

settlement_secondary

Python
settlement_secondary(C_alpha: float, H: float, t1: float, t2: float, e_p: float = None) -> float

Secondary compression S_s = C_alpha_eps * H * log10(t2/t1).

PARAMETER DESCRIPTION
C_alpha

Secondary compression index. If e_p is given this is the e-vs-log-t form C_alpha (slope of e vs log t after EOP); converted to strain form via C_alpha_eps = C_alpha / (1 + e_p). Otherwise treated directly as the strain form.

TYPE: float

H

Layer thickness (m).

TYPE: float

t1

Reference (end of primary) and target times (any consistent unit).

TYPE: float

t2

Reference (end of primary) and target times (any consistent unit).

TYPE: float

e_p

Void ratio at end-of-primary, for converting C_alpha -> C_alpha_eps.

TYPE: float DEFAULT: None

Reference

Mesri & Godlewski (1977); Das (2014) Eq. 5.51.

Source code in geoeq/design/settlement.py
Python
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
def settlement_secondary(
    C_alpha: float, H: float, t1: float, t2: float, e_p: float = None,
) -> float:
    """Secondary compression S_s = C_alpha_eps * H * log10(t2/t1).

    Parameters
    ----------
    C_alpha : float
        Secondary compression index. If ``e_p`` is given this is the
        e-vs-log-t form C_alpha (slope of e vs log t after EOP);
        converted to strain form via C_alpha_eps = C_alpha / (1 + e_p).
        Otherwise treated directly as the strain form.
    H : float
        Layer thickness (m).
    t1, t2 : float
        Reference (end of primary) and target times (any consistent unit).
    e_p : float, optional
        Void ratio at end-of-primary, for converting C_alpha -> C_alpha_eps.

    Reference
    ---------
    Mesri & Godlewski (1977); Das (2014) Eq. 5.51.
    """
    check_positive(H, "H")
    check_positive(t1, "t1")
    check_positive(t2, "t2")
    if t2 <= t1:
        raise ValueError("t2 must be > t1.")
    if e_p is not None:
        C_alpha_eps = C_alpha / (1 + e_p)
    else:
        C_alpha_eps = C_alpha
    return float(C_alpha_eps * H * np.log10(t2 / t1))

settlement_schmertmann

Python
settlement_schmertmann(q_net: float, B: float, Es: Union[float, Sequence[float]], layers: Sequence[tuple] = None, Df: float = 0.0, sigma_v_at_base: float = None, gamma: float = 18.0, t_years: float = 1.0, shape: str = 'square') -> dict

Schmertmann's strain-influence settlement for sandy soils.

Text Only
S = C1 * C2 * q_net * sum_i ( Iz_i / Es_i * dz_i )

where: C1 = 1 - 0.5 * (sigma_v' / q_net) -- embedment correction C2 = 1 + 0.2 * log10(t_years / 0.1) -- creep correction

The strain-influence diagram Iz is taken from Schmertmann et al. (1978):

  • Square / circular (axisymmetric): Iz peaks at 0.5*B (Iz_peak), = 0.1 at z=0, 0 at z=2B.
  • Strip (plane strain): Iz peaks at 1.0*B, = 0.2 at z=0, 0 at z=4B.

Iz_peak = 0.5 + 0.1 * sqrt(q_net / sigma_v_at_peak).

PARAMETER DESCRIPTION
q_net

Net surface pressure (kPa) above sigma_v at footing base.

TYPE: float

B

Footing width (m).

TYPE: float

Es

Elastic modulus of subsoil (kPa). If a single number, applied to one layer of depth 0 .. (2B or 4B). If sequence, must match layers.

TYPE: float or sequence

layers

Optional. If None and Es is scalar, the function builds a single-layer model spanning the full influence depth.

TYPE: sequence of (top, bot) in metres below the footing. DEFAULT: None

Df

Footing depth (m).

TYPE: float DEFAULT: 0.0

sigma_v_at_base

Effective vertical stress at the base of the footing (kPa). If None, estimated as gamma * Df.

TYPE: float DEFAULT: None

gamma

Soil unit weight, only used if sigma_v_at_base is None.

TYPE: float DEFAULT: 18.0

t_years

Time after construction for the creep correction (yrs). Default 1.

TYPE: float DEFAULT: 1.0

shape

'square' / 'circle' (axisymmetric) or 'strip' (plane strain).

TYPE: str DEFAULT: 'square'

RETURNS DESCRIPTION
dict

{'S': settlement_metres, 'C1', 'C2', 'Iz_peak', 'layers': [...]}
Reference

Schmertmann (1970); Schmertmann et al. (1978); Das (2014) Ch. 5.10.

Source code in geoeq/design/settlement.py
Python
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
def settlement_schmertmann(
    q_net: float, B: float, Es: Union[float, Sequence[float]],
    layers: Sequence[tuple] = None, Df: float = 0.0,
    sigma_v_at_base: float = None, gamma: float = 18.0,
    t_years: float = 1.0, shape: str = "square",
) -> dict:
    """Schmertmann's strain-influence settlement for sandy soils.

        S = C1 * C2 * q_net * sum_i ( Iz_i / Es_i * dz_i )

    where:
        C1 = 1 - 0.5 * (sigma_v' / q_net)            -- embedment correction
        C2 = 1 + 0.2 * log10(t_years / 0.1)          -- creep correction

    The strain-influence diagram Iz is taken from Schmertmann et al. (1978):

    * Square / circular (axisymmetric):  Iz peaks at 0.5*B (Iz_peak),
      = 0.1 at z=0, 0 at z=2B.
    * Strip (plane strain):             Iz peaks at 1.0*B,
      = 0.2 at z=0, 0 at z=4B.

    Iz_peak = 0.5 + 0.1 * sqrt(q_net / sigma_v_at_peak).

    Parameters
    ----------
    q_net : float
        Net surface pressure (kPa) above sigma_v at footing base.
    B : float
        Footing width (m).
    Es : float or sequence
        Elastic modulus of subsoil (kPa). If a single number, applied
        to one layer of depth 0 .. (2B or 4B). If sequence, must match ``layers``.
    layers : sequence of (top, bot) in metres below the footing.
        Optional. If None and ``Es`` is scalar, the function builds a
        single-layer model spanning the full influence depth.
    Df : float
        Footing depth (m).
    sigma_v_at_base : float, optional
        Effective vertical stress at the base of the footing (kPa).
        If None, estimated as ``gamma * Df``.
    gamma : float
        Soil unit weight, only used if sigma_v_at_base is None.
    t_years : float
        Time after construction for the creep correction (yrs). Default 1.
    shape : str
        'square' / 'circle' (axisymmetric) or 'strip' (plane strain).

    Returns
    -------
    dict
        {'S': settlement_metres, 'C1', 'C2', 'Iz_peak', 'layers': [...]}

    Reference
    ---------
    Schmertmann (1970); Schmertmann et al. (1978); Das (2014) Ch. 5.10.
    """
    check_positive(q_net, "q_net")
    check_positive(B, "B")
    shape = shape.lower()
    if shape in ("square", "circle", "circular", "axisymmetric"):
        z_peak = 0.5 * B
        z_max = 2.0 * B
        Iz_z0 = 0.1
    elif shape == "strip":
        z_peak = 1.0 * B
        z_max = 4.0 * B
        Iz_z0 = 0.2
    else:
        raise ValueError("shape must be 'square', 'circle', or 'strip'.")

    if sigma_v_at_base is None:
        sigma_v_at_base = gamma * Df

    # sigma_v at z_peak (initial effective stress).
    sigma_v_peak = sigma_v_at_base + gamma * z_peak
    Iz_peak = 0.5 + 0.1 * np.sqrt(q_net / max(sigma_v_peak, 1e-6))

    # Build the Iz(z) function: linear from (0, Iz_z0) to (z_peak, Iz_peak)
    # then linear to (z_max, 0).
    def Iz(z):
        if z < 0 or z > z_max:
            return 0.0
        if z <= z_peak:
            return Iz_z0 + (Iz_peak - Iz_z0) * (z / z_peak)
        return Iz_peak * (1 - (z - z_peak) / (z_max - z_peak))

    # Build layers if not provided.
    if layers is None:
        if np.isscalar(Es):
            layers = [(0.0, z_max)]
            Es = [float(Es)]
        else:
            raise ValueError(
                "If Es is a sequence, layers must be provided.")
    Es_arr = np.atleast_1d(np.asarray(Es, dtype=float))
    if len(Es_arr) != len(layers):
        raise ValueError("Es and layers must have the same length.")

    # Integrate sum(Iz / Es * dz) using midpoint of each layer.
    integ = 0.0
    layer_info = []
    for (top, bot), Esi in zip(layers, Es_arr):
        dz = bot - top
        z_mid = 0.5 * (top + bot)
        contribution = Iz(z_mid) / Esi * dz
        integ += contribution
        layer_info.append({"top": top, "bot": bot, "Es": Esi,
                           "Iz_mid": Iz(z_mid), "dS": contribution})

    C1 = max(0.5, 1 - 0.5 * sigma_v_at_base / q_net)
    C2 = max(1.0, 1 + 0.2 * np.log10(max(t_years, 0.1) / 0.1))
    S = C1 * C2 * q_net * integ

    return {
        "S": float(S), "C1": float(C1), "C2": float(C2),
        "Iz_peak": float(Iz_peak), "z_peak": float(z_peak),
        "z_max": float(z_max), "layers": layer_info,
    }

time_factor

Python
time_factor(U: float) -> float

Time factor Tv from average degree of consolidation U (0..1).

Approximate Terzaghi solution (Das Eq. 5.46): Tv = (pi/4) * U^2 for U < 0.6 Tv = 1.781 - 0.933 * log10(100*(1-U)) for U >= 0.6

Reference

Terzaghi (1925); Das (2014) Eq. 5.46.

Source code in geoeq/design/settlement.py
Python
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
def time_factor(U: float) -> float:
    """Time factor Tv from average degree of consolidation U (0..1).

    Approximate Terzaghi solution (Das Eq. 5.46):
        Tv = (pi/4) * U^2                   for U < 0.6
        Tv = 1.781 - 0.933 * log10(100*(1-U))   for U >= 0.6

    Reference
    ---------
    Terzaghi (1925); Das (2014) Eq. 5.46.
    """
    if not 0 <= U < 1:
        raise ValueError("U must be in [0, 1).")
    if U < 0.6:
        return (np.pi / 4) * U ** 2
    return 1.781 - 0.933 * np.log10(100 * (1 - U))

consolidation_degree

Python
consolidation_degree(Tv: float) -> float

Average degree of consolidation U (0..1) from time factor Tv.

Inverse of time_factor (Das Eq. 5.46).

Source code in geoeq/design/settlement.py
Python
356
357
358
359
360
361
362
363
364
365
def consolidation_degree(Tv: float) -> float:
    """Average degree of consolidation U (0..1) from time factor Tv.

    Inverse of ``time_factor`` (Das Eq. 5.46).
    """
    if Tv < 0:
        raise ValueError("Tv must be non-negative.")
    if Tv <= 0.286:  # boundary at U = 0.6
        return float(np.sqrt(4 * Tv / np.pi))
    return float(1 - 10 ** ((1.781 - Tv) / 0.933) / 100)

consolidation_time

Python
consolidation_time(U: float, Hdr: float, cv: float) -> float

Time to reach degree of consolidation U.

Text Only
t = Tv * Hdr^2 / cv
PARAMETER DESCRIPTION
U

Target degree of consolidation (0..1).

TYPE: float

Hdr

Drainage path length (m). Half the layer thickness for double drainage.

TYPE: float

cv

Coefficient of consolidation (m^2/s, or any consistent unit).

TYPE: float

RETURNS DESCRIPTION
t

Time (same time units as cv).

TYPE: float

Reference

Das (2014) Eq. 5.43.

Source code in geoeq/design/settlement.py
Python
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
def consolidation_time(U: float, Hdr: float, cv: float) -> float:
    """Time to reach degree of consolidation U.

        t = Tv * Hdr^2 / cv

    Parameters
    ----------
    U : float
        Target degree of consolidation (0..1).
    Hdr : float
        Drainage path length (m). Half the layer thickness for double drainage.
    cv : float
        Coefficient of consolidation (m^2/s, or any consistent unit).

    Returns
    -------
    t : float
        Time (same time units as cv).

    Reference
    ---------
    Das (2014) Eq. 5.43.
    """
    Tv = time_factor(U)
    check_positive(Hdr, "Hdr")
    check_positive(cv, "cv")
    return float(Tv * Hdr ** 2 / cv)

settlement_time_plot

Python
settlement_time_plot(ax=None, save_as: str = None)

The classical Terzaghi consolidation curve: U vs Tv (log scale).

Plots the approximate piecewise formula used by time_factor: Tv = (pi/4) U^2 for U < 0.6, Tv = 1.781 - 0.933 log10(100(1-U)) for U >= 0.6.

Source code in geoeq/design/plots.py
Python
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
def settlement_time_plot(ax=None, save_as: str = None):
    """The classical Terzaghi consolidation curve: U vs Tv (log scale).

    Plots the approximate piecewise formula used by ``time_factor``:
    Tv = (pi/4) U^2 for U < 0.6, Tv = 1.781 - 0.933 log10(100(1-U))
    for U >= 0.6.
    """
    import matplotlib.pyplot as plt
    U_arr = np.linspace(0.001, 0.99, 200)
    Tv_arr = np.array([time_factor(u) for u in U_arr])
    if ax is None:
        fig, ax = plt.subplots(figsize=(7, 5))
    else:
        fig = ax.figure
    ax.plot(Tv_arr, U_arr * 100, color="#1f3a93", linewidth=2.0)
    # Annotated reference points.
    for U in (0.5, 0.9):
        Tv = time_factor(U)
        ax.plot([Tv], [U * 100], "o", color="#c0392b", markersize=6)
        ax.annotate(f"$U={int(U*100)}\\%$\n$T_v={Tv:.3f}$",
                    xy=(Tv, U * 100), xytext=(10, -10),
                    textcoords="offset points", fontsize=9,
                    color="#c0392b")
    ax.set_xscale("log")
    ax.set_xlim(1e-3, 3)
    ax.set_ylim(0, 100)
    ax.set_xlabel(r"Time factor $T_v$")
    ax.set_ylabel(r"Average degree of consolidation $U$ (%)")
    ax.set_title("Terzaghi 1-D consolidation: $U$ vs $T_v$")
    ax.grid(True, which="both", alpha=0.3)
    if save_as:
        fig.savefig(save_as, dpi=300, bbox_inches="tight")
    return fig

Earth pressure

K0

Python
K0(phi: float, OCR: float = 1.0, method: str = 'jaky') -> float

At-rest earth pressure coefficient.

PARAMETER DESCRIPTION
phi

Effective friction angle (degrees).

TYPE: float

OCR

Over-consolidation ratio (1 for NC soils).

TYPE: float DEFAULT: 1.0

method

'jaky' -> K0 = 1 - sin(phi) (NC; Jaky 1944). 'mayne' -> K0 = (1 - sin phi) * OCR^sin(phi) (Mayne & Kulhawy 1982). 'alpan' -> Alpan (1967) for OC clays.

TYPE: str DEFAULT: 'jaky'

Reference

Jaky (1944); Mayne & Kulhawy (1982).

Source code in geoeq/design/earth_pressure.py
Python
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def K0(phi: float, OCR: float = 1.0, method: str = "jaky") -> float:
    """At-rest earth pressure coefficient.

    Parameters
    ----------
    phi : float
        Effective friction angle (degrees).
    OCR : float
        Over-consolidation ratio (1 for NC soils).
    method : str
        'jaky'  -> K0 = 1 - sin(phi) (NC; Jaky 1944).
        'mayne' -> K0 = (1 - sin phi) * OCR^sin(phi)  (Mayne & Kulhawy 1982).
        'alpan' -> Alpan (1967) for OC clays.

    Reference
    ---------
    Jaky (1944); Mayne & Kulhawy (1982).
    """
    check_range(phi, "phi", 0, 50)
    check_positive(OCR, "OCR")
    phi_rad = np.radians(phi)
    method = method.lower()
    if method == "jaky":
        return float(1 - np.sin(phi_rad))
    if method == "mayne":
        return float((1 - np.sin(phi_rad)) * OCR ** np.sin(phi_rad))
    if method == "alpan":
        return float(0.19 + 0.233 * np.log10(max(phi, 1)))  # rough
    raise ValueError("method must be 'jaky', 'mayne', or 'alpan'.")

Ka

Python
Ka(phi: float, delta: float = 0.0, alpha: float = 0.0, beta: float = 0.0, method: str = 'rankine') -> float

Active earth pressure coefficient.

PARAMETER DESCRIPTION
phi

Effective friction angle of backfill (degrees).

TYPE: float

delta

Wall-soil friction angle (degrees). Default 0 (smooth wall).

TYPE: float DEFAULT: 0.0

alpha

Wall back-face inclination from vertical (degrees, positive = wall leans toward backfill). Default 0.

TYPE: float DEFAULT: 0.0

beta

Backfill slope angle (degrees, positive = ascending). Default 0.

TYPE: float DEFAULT: 0.0

method

'rankine' -- requires delta=alpha=0; uses Rankine (1857) formula. 'coulomb' -- general (Coulomb 1776).

TYPE: str DEFAULT: 'rankine'

Reference

Rankine (1857); Coulomb (1776); Das (2014) Ch. 7-8.

Source code in geoeq/design/earth_pressure.py
Python
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def Ka(phi: float, delta: float = 0.0, alpha: float = 0.0,
       beta: float = 0.0, method: str = "rankine") -> float:
    """Active earth pressure coefficient.

    Parameters
    ----------
    phi : float
        Effective friction angle of backfill (degrees).
    delta : float
        Wall-soil friction angle (degrees). Default 0 (smooth wall).
    alpha : float
        Wall back-face inclination from vertical (degrees, positive = wall
        leans toward backfill). Default 0.
    beta : float
        Backfill slope angle (degrees, positive = ascending). Default 0.
    method : str
        'rankine' -- requires delta=alpha=0; uses Rankine (1857) formula.
        'coulomb' -- general (Coulomb 1776).

    Reference
    ---------
    Rankine (1857); Coulomb (1776); Das (2014) Ch. 7-8.
    """
    check_range(phi, "phi", 0, 50)
    phi_r = np.radians(phi)
    delta_r = np.radians(delta)
    alpha_r = np.radians(alpha)
    beta_r = np.radians(beta)
    method = method.lower()

    if method == "rankine":
        if abs(beta) < 1e-6:
            return float(np.tan(np.pi / 4 - phi_r / 2) ** 2)
        # Sloping backfill, vertical smooth wall (Das Eq. 7.21).
        cos_b = np.cos(beta_r)
        cos_p = np.cos(phi_r)
        num = cos_b - np.sqrt(cos_b ** 2 - cos_p ** 2)
        den = cos_b + np.sqrt(cos_b ** 2 - cos_p ** 2)
        return float(cos_b * num / den)

    if method == "coulomb":
        # Coulomb's general Ka (Das Eq. 8.5).
        num = np.cos(phi_r - alpha_r) ** 2
        sin_term = (np.sqrt(np.sin(phi_r + delta_r)
                            * np.sin(phi_r - beta_r)
                            / (np.cos(alpha_r - beta_r)
                               * np.cos(alpha_r + delta_r))))
        den = (np.cos(alpha_r) ** 2 * np.cos(alpha_r + delta_r)
               * (1 + sin_term) ** 2)
        return float(num / den)

    raise ValueError("method must be 'rankine' or 'coulomb'.")

Kp

Python
Kp(phi: float, delta: float = 0.0, alpha: float = 0.0, beta: float = 0.0, method: str = 'rankine') -> float

Passive earth pressure coefficient.

See Ka for parameters. Same conventions; passive coefficient is the sign-reversed version of active.

Reference

Rankine (1857); Coulomb (1776); Das (2014) Ch. 7-8.

Source code in geoeq/design/earth_pressure.py
Python
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
def Kp(phi: float, delta: float = 0.0, alpha: float = 0.0,
       beta: float = 0.0, method: str = "rankine") -> float:
    """Passive earth pressure coefficient.

    See ``Ka`` for parameters. Same conventions; passive coefficient is the
    sign-reversed version of active.

    Reference
    ---------
    Rankine (1857); Coulomb (1776); Das (2014) Ch. 7-8.
    """
    check_range(phi, "phi", 0, 50)
    phi_r = np.radians(phi)
    delta_r = np.radians(delta)
    alpha_r = np.radians(alpha)
    beta_r = np.radians(beta)
    method = method.lower()
    if method == "rankine":
        if abs(beta) < 1e-6:
            return float(np.tan(np.pi / 4 + phi_r / 2) ** 2)
        cos_b = np.cos(beta_r)
        cos_p = np.cos(phi_r)
        num = cos_b + np.sqrt(cos_b ** 2 - cos_p ** 2)
        den = cos_b - np.sqrt(cos_b ** 2 - cos_p ** 2)
        return float(cos_b * num / den)
    if method == "coulomb":
        num = np.cos(phi_r + alpha_r) ** 2
        sin_term = (np.sqrt(np.sin(phi_r + delta_r)
                            * np.sin(phi_r + beta_r)
                            / (np.cos(alpha_r - beta_r)
                               * np.cos(alpha_r - delta_r))))
        den = (np.cos(alpha_r) ** 2 * np.cos(alpha_r - delta_r)
               * (1 - sin_term) ** 2)
        return float(num / den)
    raise ValueError("method must be 'rankine' or 'coulomb'.")

earth_pressure

Lateral earth pressure.

K0 (at-rest), Ka (active), Kp (passive) -- Rankine and Coulomb theories, plus distributions including water, surcharge, and cohesion (tension cracks).

References
  • Coulomb, C. A. (1776). "Essai sur une application des regles de maximis et minimis a quelques problemes de statique relatifs a l'architecture." Memoires de Mathematique et de Physique, Paris.
  • Rankine, W. J. M. (1857). "On the stability of loose earth." Phil. Trans. Royal Society, 147, 9-27.
  • Jaky, J. (1944). "The coefficient of earth pressure at rest." J. Society Hungarian Architects and Engineers, 22, 355-358.
  • Mayne, P. W., Kulhawy, F. H. (1982). "K0-OCR relationships in soil." J. Geotech. Eng., ASCE, 108(GT6), 851-872.
  • Das (2014), Ch. 7-8.

K0

Python
K0(phi: float, OCR: float = 1.0, method: str = 'jaky') -> float

At-rest earth pressure coefficient.

PARAMETER DESCRIPTION
phi

Effective friction angle (degrees).

TYPE: float

OCR

Over-consolidation ratio (1 for NC soils).

TYPE: float DEFAULT: 1.0

method

'jaky' -> K0 = 1 - sin(phi) (NC; Jaky 1944). 'mayne' -> K0 = (1 - sin phi) * OCR^sin(phi) (Mayne & Kulhawy 1982). 'alpan' -> Alpan (1967) for OC clays.

TYPE: str DEFAULT: 'jaky'

Reference

Jaky (1944); Mayne & Kulhawy (1982).

Source code in geoeq/design/earth_pressure.py
Python
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def K0(phi: float, OCR: float = 1.0, method: str = "jaky") -> float:
    """At-rest earth pressure coefficient.

    Parameters
    ----------
    phi : float
        Effective friction angle (degrees).
    OCR : float
        Over-consolidation ratio (1 for NC soils).
    method : str
        'jaky'  -> K0 = 1 - sin(phi) (NC; Jaky 1944).
        'mayne' -> K0 = (1 - sin phi) * OCR^sin(phi)  (Mayne & Kulhawy 1982).
        'alpan' -> Alpan (1967) for OC clays.

    Reference
    ---------
    Jaky (1944); Mayne & Kulhawy (1982).
    """
    check_range(phi, "phi", 0, 50)
    check_positive(OCR, "OCR")
    phi_rad = np.radians(phi)
    method = method.lower()
    if method == "jaky":
        return float(1 - np.sin(phi_rad))
    if method == "mayne":
        return float((1 - np.sin(phi_rad)) * OCR ** np.sin(phi_rad))
    if method == "alpan":
        return float(0.19 + 0.233 * np.log10(max(phi, 1)))  # rough
    raise ValueError("method must be 'jaky', 'mayne', or 'alpan'.")

Ka

Python
Ka(phi: float, delta: float = 0.0, alpha: float = 0.0, beta: float = 0.0, method: str = 'rankine') -> float

Active earth pressure coefficient.

PARAMETER DESCRIPTION
phi

Effective friction angle of backfill (degrees).

TYPE: float

delta

Wall-soil friction angle (degrees). Default 0 (smooth wall).

TYPE: float DEFAULT: 0.0

alpha

Wall back-face inclination from vertical (degrees, positive = wall leans toward backfill). Default 0.

TYPE: float DEFAULT: 0.0

beta

Backfill slope angle (degrees, positive = ascending). Default 0.

TYPE: float DEFAULT: 0.0

method

'rankine' -- requires delta=alpha=0; uses Rankine (1857) formula. 'coulomb' -- general (Coulomb 1776).

TYPE: str DEFAULT: 'rankine'

Reference

Rankine (1857); Coulomb (1776); Das (2014) Ch. 7-8.

Source code in geoeq/design/earth_pressure.py
Python
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def Ka(phi: float, delta: float = 0.0, alpha: float = 0.0,
       beta: float = 0.0, method: str = "rankine") -> float:
    """Active earth pressure coefficient.

    Parameters
    ----------
    phi : float
        Effective friction angle of backfill (degrees).
    delta : float
        Wall-soil friction angle (degrees). Default 0 (smooth wall).
    alpha : float
        Wall back-face inclination from vertical (degrees, positive = wall
        leans toward backfill). Default 0.
    beta : float
        Backfill slope angle (degrees, positive = ascending). Default 0.
    method : str
        'rankine' -- requires delta=alpha=0; uses Rankine (1857) formula.
        'coulomb' -- general (Coulomb 1776).

    Reference
    ---------
    Rankine (1857); Coulomb (1776); Das (2014) Ch. 7-8.
    """
    check_range(phi, "phi", 0, 50)
    phi_r = np.radians(phi)
    delta_r = np.radians(delta)
    alpha_r = np.radians(alpha)
    beta_r = np.radians(beta)
    method = method.lower()

    if method == "rankine":
        if abs(beta) < 1e-6:
            return float(np.tan(np.pi / 4 - phi_r / 2) ** 2)
        # Sloping backfill, vertical smooth wall (Das Eq. 7.21).
        cos_b = np.cos(beta_r)
        cos_p = np.cos(phi_r)
        num = cos_b - np.sqrt(cos_b ** 2 - cos_p ** 2)
        den = cos_b + np.sqrt(cos_b ** 2 - cos_p ** 2)
        return float(cos_b * num / den)

    if method == "coulomb":
        # Coulomb's general Ka (Das Eq. 8.5).
        num = np.cos(phi_r - alpha_r) ** 2
        sin_term = (np.sqrt(np.sin(phi_r + delta_r)
                            * np.sin(phi_r - beta_r)
                            / (np.cos(alpha_r - beta_r)
                               * np.cos(alpha_r + delta_r))))
        den = (np.cos(alpha_r) ** 2 * np.cos(alpha_r + delta_r)
               * (1 + sin_term) ** 2)
        return float(num / den)

    raise ValueError("method must be 'rankine' or 'coulomb'.")

Kp

Python
Kp(phi: float, delta: float = 0.0, alpha: float = 0.0, beta: float = 0.0, method: str = 'rankine') -> float

Passive earth pressure coefficient.

See Ka for parameters. Same conventions; passive coefficient is the sign-reversed version of active.

Reference

Rankine (1857); Coulomb (1776); Das (2014) Ch. 7-8.

Source code in geoeq/design/earth_pressure.py
Python
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
def Kp(phi: float, delta: float = 0.0, alpha: float = 0.0,
       beta: float = 0.0, method: str = "rankine") -> float:
    """Passive earth pressure coefficient.

    See ``Ka`` for parameters. Same conventions; passive coefficient is the
    sign-reversed version of active.

    Reference
    ---------
    Rankine (1857); Coulomb (1776); Das (2014) Ch. 7-8.
    """
    check_range(phi, "phi", 0, 50)
    phi_r = np.radians(phi)
    delta_r = np.radians(delta)
    alpha_r = np.radians(alpha)
    beta_r = np.radians(beta)
    method = method.lower()
    if method == "rankine":
        if abs(beta) < 1e-6:
            return float(np.tan(np.pi / 4 + phi_r / 2) ** 2)
        cos_b = np.cos(beta_r)
        cos_p = np.cos(phi_r)
        num = cos_b + np.sqrt(cos_b ** 2 - cos_p ** 2)
        den = cos_b - np.sqrt(cos_b ** 2 - cos_p ** 2)
        return float(cos_b * num / den)
    if method == "coulomb":
        num = np.cos(phi_r + alpha_r) ** 2
        sin_term = (np.sqrt(np.sin(phi_r + delta_r)
                            * np.sin(phi_r + beta_r)
                            / (np.cos(alpha_r - beta_r)
                               * np.cos(alpha_r - delta_r))))
        den = (np.cos(alpha_r) ** 2 * np.cos(alpha_r - delta_r)
               * (1 - sin_term) ** 2)
        return float(num / den)
    raise ValueError("method must be 'rankine' or 'coulomb'.")

earth_pressure

Python
earth_pressure(gamma: float, H: float, phi: float, c: float = 0.0, kind: str = 'active', surcharge: float = 0.0, water_table: float = None, K_method: str = 'rankine') -> dict

Lateral pressure distribution on a vertical wall of height H.

Earth pressure at depth z: sigma_h = K * gamma * z (cohesionless, no surcharge) + K * surcharge - 2 c sqrt(K) (cohesion subtraction in active case)

Returns the resultant force P (kN/m) and its line of action.

PARAMETER DESCRIPTION
gamma

Bulk unit weight (kN/m^3).

TYPE: float

H

Wall height (m).

TYPE: float

phi

Friction angle (degrees).

TYPE: float

c

Cohesion (kPa).

TYPE: float DEFAULT: 0.0

kind

'active' | 'passive' | 'at_rest'.

TYPE: str DEFAULT: 'active'

surcharge

Uniform surface surcharge q (kPa).

TYPE: float DEFAULT: 0.0

water_table

Depth of water table below the top of the wall (m). If given, an additional hydrostatic pressure (gamma_w*(z - z_w)) is added and the effective unit weight below z_w is gamma - gamma_w.

TYPE: float DEFAULT: None

RETURNS DESCRIPTION
dict

{'K': ..., 'P_total': ..., 'P_water': ..., 'P_soil': ..., 'h_point': ...}.
Reference

Das (2014), Ch. 7-8.

Source code in geoeq/design/earth_pressure.py
Python
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
def earth_pressure(
    gamma: float, H: float, phi: float, c: float = 0.0,
    kind: str = "active", surcharge: float = 0.0,
    water_table: float = None, K_method: str = "rankine",
) -> dict:
    """Lateral pressure distribution on a vertical wall of height H.

    Earth pressure at depth z:
        sigma_h = K * gamma * z   (cohesionless, no surcharge)
                + K * surcharge
                - 2 c sqrt(K)     (cohesion subtraction in active case)

    Returns the resultant force P (kN/m) and its line of action.

    Parameters
    ----------
    gamma : float
        Bulk unit weight (kN/m^3).
    H : float
        Wall height (m).
    phi : float
        Friction angle (degrees).
    c : float
        Cohesion (kPa).
    kind : str
        'active' | 'passive' | 'at_rest'.
    surcharge : float
        Uniform surface surcharge q (kPa).
    water_table : float, optional
        Depth of water table below the top of the wall (m). If given,
        an additional hydrostatic pressure (gamma_w*(z - z_w)) is added
        and the effective unit weight below z_w is gamma - gamma_w.

    Returns
    -------
    dict
        ``{'K': ..., 'P_total': ..., 'P_water': ..., 'P_soil': ...,
        'h_point': ...}``.

    Reference
    ---------
    Das (2014), Ch. 7-8.
    """
    check_positive(H, "H")
    check_non_negative(c, "c")
    check_positive(gamma, "gamma")
    kind = kind.lower()

    if kind == "active":
        K = Ka(phi, method=K_method)
    elif kind == "passive":
        K = Kp(phi, method=K_method)
    elif kind in ("at_rest", "atrest", "rest"):
        K = K0(phi)
    else:
        raise ValueError("kind must be 'active', 'passive', or 'at_rest'.")

    # Cohesion subtraction (active) or addition (passive).
    sign = -1.0 if kind == "active" else (1.0 if kind == "passive" else 0.0)
    cohesion_term = sign * 2 * c * np.sqrt(K)  # subtracted from active

    # Soil pressure -- single layer for now.
    P_soil_top = K * surcharge + cohesion_term
    P_soil_bot = K * surcharge + K * gamma * H + cohesion_term
    P_soil = 0.5 * (P_soil_top + P_soil_bot) * H

    # Water (if any).
    if water_table is not None and water_table < H:
        from geoeq.core.constants import GAMMA_WATER
        h_w = H - water_table
        P_water = 0.5 * GAMMA_WATER * h_w ** 2
    else:
        P_water = 0.0
    P_total = P_soil + P_water

    # Line of action (from base): for a trapezoid, the centroid is at
    # H * (2*top + bot) / (3*(top+bot)).
    if (P_soil_top + P_soil_bot) > 0:
        h_point = H * (P_soil_top + 2 * P_soil_bot) / (
            3 * (P_soil_top + P_soil_bot))
    else:
        h_point = H / 3
    return {
        "K": float(K), "P_soil": float(P_soil),
        "P_water": float(P_water), "P_total": float(P_total),
        "h_point": float(h_point),
        "sigma_top": float(P_soil_top), "sigma_bot": float(P_soil_bot),
    }

tension_crack_depth

Python
tension_crack_depth(c: float, gamma: float, Ka_value: float = None, phi: float = 0.0) -> float

Depth of tension crack in cohesive soil behind a wall.

Text Only
z_c = 2 c / (gamma * sqrt(Ka))
PARAMETER DESCRIPTION
c

Cohesion (kPa).

TYPE: float

gamma

Unit weight (kN/m^3).

TYPE: float

Ka_value

Active earth pressure coefficient. If None, computed from phi.

TYPE: float DEFAULT: None

phi

Friction angle (degrees). Used only if Ka_value is None.

TYPE: float DEFAULT: 0.0

Reference

Das (2014) Eq. 7.18.

Source code in geoeq/design/earth_pressure.py
Python
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
def tension_crack_depth(c: float, gamma: float, Ka_value: float = None,
                        phi: float = 0.0) -> float:
    """Depth of tension crack in cohesive soil behind a wall.

        z_c = 2 c / (gamma * sqrt(Ka))

    Parameters
    ----------
    c : float
        Cohesion (kPa).
    gamma : float
        Unit weight (kN/m^3).
    Ka_value : float, optional
        Active earth pressure coefficient. If None, computed from ``phi``.
    phi : float
        Friction angle (degrees). Used only if Ka_value is None.

    Reference
    ---------
    Das (2014) Eq. 7.18.
    """
    check_positive(c, "c")
    check_positive(gamma, "gamma")
    if Ka_value is None:
        Ka_value = Ka(phi, method="rankine")
    return float(2 * c / (gamma * np.sqrt(Ka_value)))

earth_pressure_plot

Python
earth_pressure_plot(gamma: float, H: float, phi: float, c: float = 0.0, kind: str = 'active', save_as: str = None)

Visualize the lateral pressure distribution behind a wall.

Returns the Matplotlib figure.

Source code in geoeq/design/earth_pressure.py
Python
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
def earth_pressure_plot(gamma: float, H: float, phi: float, c: float = 0.0,
                        kind: str = "active", save_as: str = None):
    """Visualize the lateral pressure distribution behind a wall.

    Returns the Matplotlib figure.
    """
    import matplotlib.pyplot as plt
    z = np.linspace(0, H, 100)
    kind = kind.lower()
    K = (Ka(phi) if kind == "active" else
         Kp(phi) if kind == "passive" else K0(phi))
    sigma = K * gamma * z - 2 * c * np.sqrt(K) * (kind == "active") + \
        2 * c * np.sqrt(K) * (kind == "passive")
    fig, ax = plt.subplots(figsize=(5, 8))
    ax.fill_betweenx(z, 0, sigma, alpha=0.3, color="#c0392b")
    ax.plot(sigma, z, color="#c0392b", linewidth=1.8)
    ax.invert_yaxis()
    ax.set_xlabel("Lateral pressure $\\sigma_h$ (kPa)")
    ax.set_ylabel("Depth (m)")
    ax.set_title(f"{kind.title()} earth pressure (K = {K:.3f})")
    ax.grid(alpha=0.3)
    if save_as:
        fig.savefig(save_as, dpi=300, bbox_inches="tight")
    return fig

tension_crack_depth

Python
tension_crack_depth(c: float, gamma: float, Ka_value: float = None, phi: float = 0.0) -> float

Depth of tension crack in cohesive soil behind a wall.

Text Only
z_c = 2 c / (gamma * sqrt(Ka))
PARAMETER DESCRIPTION
c

Cohesion (kPa).

TYPE: float

gamma

Unit weight (kN/m^3).

TYPE: float

Ka_value

Active earth pressure coefficient. If None, computed from phi.

TYPE: float DEFAULT: None

phi

Friction angle (degrees). Used only if Ka_value is None.

TYPE: float DEFAULT: 0.0

Reference

Das (2014) Eq. 7.18.

Source code in geoeq/design/earth_pressure.py
Python
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
def tension_crack_depth(c: float, gamma: float, Ka_value: float = None,
                        phi: float = 0.0) -> float:
    """Depth of tension crack in cohesive soil behind a wall.

        z_c = 2 c / (gamma * sqrt(Ka))

    Parameters
    ----------
    c : float
        Cohesion (kPa).
    gamma : float
        Unit weight (kN/m^3).
    Ka_value : float, optional
        Active earth pressure coefficient. If None, computed from ``phi``.
    phi : float
        Friction angle (degrees). Used only if Ka_value is None.

    Reference
    ---------
    Das (2014) Eq. 7.18.
    """
    check_positive(c, "c")
    check_positive(gamma, "gamma")
    if Ka_value is None:
        Ka_value = Ka(phi, method="rankine")
    return float(2 * c / (gamma * np.sqrt(Ka_value)))

earth_pressure_plot

Python
earth_pressure_plot(gamma: float, H: float, phi: float, c: float = 0.0, kind: str = 'active', save_as: str = None)

Visualize the lateral pressure distribution behind a wall.

Returns the Matplotlib figure.

Source code in geoeq/design/earth_pressure.py
Python
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
def earth_pressure_plot(gamma: float, H: float, phi: float, c: float = 0.0,
                        kind: str = "active", save_as: str = None):
    """Visualize the lateral pressure distribution behind a wall.

    Returns the Matplotlib figure.
    """
    import matplotlib.pyplot as plt
    z = np.linspace(0, H, 100)
    kind = kind.lower()
    K = (Ka(phi) if kind == "active" else
         Kp(phi) if kind == "passive" else K0(phi))
    sigma = K * gamma * z - 2 * c * np.sqrt(K) * (kind == "active") + \
        2 * c * np.sqrt(K) * (kind == "passive")
    fig, ax = plt.subplots(figsize=(5, 8))
    ax.fill_betweenx(z, 0, sigma, alpha=0.3, color="#c0392b")
    ax.plot(sigma, z, color="#c0392b", linewidth=1.8)
    ax.invert_yaxis()
    ax.set_xlabel("Lateral pressure $\\sigma_h$ (kPa)")
    ax.set_ylabel("Depth (m)")
    ax.set_title(f"{kind.title()} earth pressure (K = {K:.3f})")
    ax.grid(alpha=0.3)
    if save_as:
        fig.savefig(save_as, dpi=300, bbox_inches="tight")
    return fig

Walls

wall_overturning

Python
wall_overturning(resisting_moments: Sequence[float], driving_moments: Sequence[float]) -> dict

Factor of safety against overturning.

Text Only
FS_overturning = sum(M_R) / sum(M_O)

Returns dict with totals and FS.

Reference

Das (2014) Eq. 8.4.

Source code in geoeq/design/walls.py
Python
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def wall_overturning(
    resisting_moments: Sequence[float],
    driving_moments: Sequence[float],
) -> dict:
    """Factor of safety against overturning.

        FS_overturning = sum(M_R) / sum(M_O)

    Returns dict with totals and FS.

    Reference
    ---------
    Das (2014) Eq. 8.4.
    """
    M_R = float(np.sum(np.atleast_1d(resisting_moments)))
    M_O = float(np.sum(np.atleast_1d(driving_moments)))
    if M_O <= 0:
        raise ValueError("Total driving moment must be > 0.")
    return {"M_resisting": M_R, "M_driving": M_O, "FS": M_R / M_O}

wall_sliding

Python
wall_sliding(horizontal_forces: Sequence[float], vertical_forces: Sequence[float], mu: float = None, delta: float = None, c_base: float = 0.0, B: float = None, Pp: float = 0.0) -> dict

Factor of safety against sliding at the wall base.

Text Only
FS_sliding = (sum_V * tan(delta) + c_a * B + Pp) / sum_H
PARAMETER DESCRIPTION
horizontal_forces

Forces pushing the wall horizontally (kN/m of wall).

TYPE: sequence

vertical_forces

Vertical forces resisting sliding (kN/m of wall).

TYPE: sequence

mu

Coefficient of friction between base and soil. If given, used as tan(delta). Otherwise computed from delta.

TYPE: float DEFAULT: None

delta

Base friction angle (degrees).

TYPE: float DEFAULT: None

c_base

Cohesion at base (kPa). Typically c_a ~ 0.5 to ⅔ of c'.

TYPE: float DEFAULT: 0.0

B

Base width (m). Required if c_base > 0.

TYPE: float DEFAULT: None

Pp

Passive resistance at the toe (kN/m). Default 0.

TYPE: float DEFAULT: 0.0

Reference

Das (2014) Eq. 8.6.

Source code in geoeq/design/walls.py
Python
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def wall_sliding(
    horizontal_forces: Sequence[float],
    vertical_forces: Sequence[float],
    mu: float = None, delta: float = None,
    c_base: float = 0.0, B: float = None, Pp: float = 0.0,
) -> dict:
    """Factor of safety against sliding at the wall base.

        FS_sliding = (sum_V * tan(delta) + c_a * B + Pp) / sum_H

    Parameters
    ----------
    horizontal_forces : sequence
        Forces pushing the wall horizontally (kN/m of wall).
    vertical_forces : sequence
        Vertical forces resisting sliding (kN/m of wall).
    mu : float, optional
        Coefficient of friction between base and soil. If given, used as
        tan(delta). Otherwise computed from ``delta``.
    delta : float, optional
        Base friction angle (degrees).
    c_base : float
        Cohesion at base (kPa). Typically c_a ~ 0.5 to 2/3 of c'.
    B : float
        Base width (m). Required if c_base > 0.
    Pp : float
        Passive resistance at the toe (kN/m). Default 0.

    Reference
    ---------
    Das (2014) Eq. 8.6.
    """
    H = float(np.sum(np.atleast_1d(horizontal_forces)))
    V = float(np.sum(np.atleast_1d(vertical_forces)))
    check_positive(H, "sum(H)")
    if mu is None:
        if delta is None:
            raise ValueError("Provide mu or delta.")
        mu = np.tan(np.radians(delta))
    R_friction = V * mu
    R_cohesion = c_base * B if (c_base > 0 and B) else 0.0
    R_total = R_friction + R_cohesion + Pp
    return {
        "sum_V": V, "sum_H": H,
        "R_friction": R_friction, "R_cohesion": R_cohesion,
        "Pp": Pp, "R_total": R_total, "FS": R_total / H,
    }

wall_bearing

Python
wall_bearing(V: float, M_net: float, B: float, q_ult: float = None) -> dict

Bearing-pressure distribution under a wall foundation.

Computes the eccentricity e, maximum and minimum toe pressures, and -- if q_ult is given -- FS_bearing = q_ult / q_max.

PARAMETER DESCRIPTION
V

Total vertical load (kN/m).

TYPE: float

M_net

Net moment about the centreline of the base (kN.m/m). Positive = overturning sense.

TYPE: float

B

Base width (m).

TYPE: float

q_ult

Ultimate bearing capacity (kPa). If given, returns FS.

TYPE: float DEFAULT: None

RETURNS DESCRIPTION
dict

{'e', 'q_max', 'q_min', 'q_ult', 'FS', 'within_kern'}.
Reference

Das (2014) Eq. 8.8-8.10.

Source code in geoeq/design/walls.py
Python
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
def wall_bearing(
    V: float, M_net: float, B: float, q_ult: float = None,
) -> dict:
    """Bearing-pressure distribution under a wall foundation.

    Computes the eccentricity ``e``, maximum and minimum toe pressures,
    and -- if ``q_ult`` is given -- FS_bearing = q_ult / q_max.

    Parameters
    ----------
    V : float
        Total vertical load (kN/m).
    M_net : float
        Net moment about the centreline of the base (kN.m/m).
        Positive = overturning sense.
    B : float
        Base width (m).
    q_ult : float, optional
        Ultimate bearing capacity (kPa). If given, returns FS.

    Returns
    -------
    dict
        ``{'e', 'q_max', 'q_min', 'q_ult', 'FS', 'within_kern'}``.

    Reference
    ---------
    Das (2014) Eq. 8.8-8.10.
    """
    check_positive(V, "V")
    check_positive(B, "B")
    e = abs(M_net) / V
    in_kern = e <= B / 6
    if in_kern:
        q_max = (V / B) * (1 + 6 * e / B)
        q_min = (V / B) * (1 - 6 * e / B)
    else:
        # Triangular distribution -- only 3*(B/2 - e) of the base is in contact.
        q_max = (4 * V) / (3 * (B - 2 * e))
        q_min = 0.0
    out = {"e": float(e), "q_max": float(q_max), "q_min": float(q_min),
           "within_kern": bool(in_kern)}
    if q_ult is not None:
        check_positive(q_ult, "q_ult")
        out["q_ult"] = float(q_ult)
        out["FS"] = float(q_ult / q_max)
    return out

sheet_pile

Python
sheet_pile(gamma: float, H: float, phi: float, c: float = 0.0, gamma_sub: float = None, water_table: float = None, kind: str = 'cantilever') -> dict

Embedment depth of a sheet pile by Blum's simplified method.

Cantilever wall in cohesionless soil: Embedment depth D ~ found by satisfying moment equilibrium about the point of rotation. Closed-form approximation:

Text Only
    D ~ 1.2 to 1.5 * Dmin    where Dmin from theoretical balance.

Practical (Das Eq. 14.6): for a cantilever sheet pile in dry cohesionless backfill of height H above the dredge line:

Text Only
sigma_a = Ka * gamma * H   at the dredge line
gamma_b = gamma below dredge line (gamma_sub if submerged)

Total embedment D ~ H * [ (Kp - Ka) / (gamma * Kp) ]^something

Returns the theoretical embedment D_theory and a design D = 1.3 * D_theory.

PARAMETER DESCRIPTION
gamma

Unit weight of backfill (kN/m^3).

TYPE: float

H

Free height above the dredge line (m).

TYPE: float

phi

Friction angle (degrees).

TYPE: float

c

Cohesion (kPa). For c-phi soils, applies Bell's solution; for c=0 uses the Blum equations.

TYPE: float DEFAULT: 0.0

gamma_sub

Submerged unit weight below the water table (kN/m^3).

TYPE: float DEFAULT: None

water_table

Depth of the water table from the top (m).

TYPE: float DEFAULT: None

kind

'cantilever' (only currently supported).

TYPE: str DEFAULT: 'cantilever'

RETURNS DESCRIPTION
dict

{'D_theory', 'D_design', 'Ka', 'Kp'}.
Reference

Blum (1931); Das (2014) Ch. 14.

Source code in geoeq/design/walls.py
Python
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
def sheet_pile(
    gamma: float, H: float, phi: float, c: float = 0.0,
    gamma_sub: float = None, water_table: float = None,
    kind: str = "cantilever",
) -> dict:
    """Embedment depth of a sheet pile by Blum's simplified method.

    Cantilever wall in cohesionless soil:
        Embedment depth D ~ found by satisfying moment equilibrium about
        the point of rotation. Closed-form approximation:

            D ~ 1.2 to 1.5 * Dmin    where Dmin from theoretical balance.

    Practical (Das Eq. 14.6): for a cantilever sheet pile in dry
    cohesionless backfill of height H above the dredge line:

        sigma_a = Ka * gamma * H   at the dredge line
        gamma_b = gamma below dredge line (gamma_sub if submerged)

        Total embedment D ~ H * [ (Kp - Ka) / (gamma * Kp) ]^something

    Returns the theoretical embedment D_theory and a design D = 1.3 * D_theory.

    Parameters
    ----------
    gamma : float
        Unit weight of backfill (kN/m^3).
    H : float
        Free height above the dredge line (m).
    phi : float
        Friction angle (degrees).
    c : float
        Cohesion (kPa). For c-phi soils, applies Bell's solution; for
        c=0 uses the Blum equations.
    gamma_sub : float, optional
        Submerged unit weight below the water table (kN/m^3).
    water_table : float, optional
        Depth of the water table from the top (m).
    kind : str
        'cantilever' (only currently supported).

    Returns
    -------
    dict
        ``{'D_theory', 'D_design', 'Ka', 'Kp'}``.

    Reference
    ---------
    Blum (1931); Das (2014) Ch. 14.
    """
    check_positive(gamma, "gamma")
    check_positive(H, "H")
    check_non_negative(c, "c")
    if kind.lower() != "cantilever":
        raise NotImplementedError("Only 'cantilever' is implemented.")
    Ka_ = Ka(phi)
    Kp_ = Kp(phi)
    # Driving force per metre at and below dredge line.
    # Cohesionless: sigma_a at dredge = Ka*gamma*H, sigma_p_resisting at base.
    # Using Bowles' simplified equation for D (no water for now).
    # D_theory derived from cubic moment balance; approximate via:
    #   D = H * sqrt(Ka / (Kp - Ka))   -- Bowles "rough first estimate"
    if Kp_ - Ka_ <= 0:
        raise ValueError("Invalid soil: Kp <= Ka.")
    D_theory = H * np.sqrt(Ka_ / (Kp_ - Ka_))
    D_design = 1.3 * D_theory
    return {
        "Ka": float(Ka_), "Kp": float(Kp_),
        "D_theory": float(D_theory),
        "D_design": float(D_design),
        "method": "Blum cantilever (approximate)",
    }

Pile design

pile_end_bearing

Python
pile_end_bearing(phi: float = None, sigma_v_eff: float = None, c: float = None, method: str = 'meyerhof', Su: float = None) -> dict

Unit tip resistance q_p (kPa) at the pile base.

Two regimes:

  • Sand (drained): q_p = sigma'_v * N_q* Meyerhof: N_q* from chart -- approximation N_q* = exp(pi tan phi) * tan^2(45+phi/2) but capped per Meyerhof (1976) Fig 9.11. Vesic: N_q* = (1 + 2 K_0) / 3 * (tan phi)^? -- simplified here.
  • Clay (undrained): q_p = 9 * Su (Skempton)
PARAMETER DESCRIPTION
phi

Friction angle (degrees). Required for sand.

TYPE: float DEFAULT: None

sigma_v_eff

Effective vertical stress at pile tip (kPa). Required for sand.

TYPE: float DEFAULT: None

Su

Undrained shear strength at tip (kPa). Required for clay.

TYPE: float DEFAULT: None

method

'meyerhof' | 'vesic' | 'skempton' (clay -- alias).

TYPE: str DEFAULT: 'meyerhof'

c

Optional drained cohesion (kPa). Adds c * Nc* contribution.

TYPE: float DEFAULT: None

RETURNS DESCRIPTION
dict

{'q_p': kPa, 'Nq*': ..., 'method': ...}.
Reference

Meyerhof (1976); Vesic (1977); Das (2014) Ch. 9.

Source code in geoeq/design/piles.py
Python
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def pile_end_bearing(
    phi: float = None, sigma_v_eff: float = None, c: float = None,
    method: str = "meyerhof", Su: float = None,
) -> dict:
    """Unit tip resistance q_p (kPa) at the pile base.

    Two regimes:

    * **Sand (drained):** q_p = sigma'_v * N_q*
        Meyerhof:  N_q* from chart -- approximation N_q* = exp(pi tan phi) * tan^2(45+phi/2) but
                   capped per Meyerhof (1976) Fig 9.11.
        Vesic:     N_q* = (1 + 2 K_0) / 3 * (tan phi)^? -- simplified here.
    * **Clay (undrained):** q_p = 9 * Su    (Skempton)

    Parameters
    ----------
    phi : float
        Friction angle (degrees). Required for sand.
    sigma_v_eff : float
        Effective vertical stress at pile tip (kPa). Required for sand.
    Su : float
        Undrained shear strength at tip (kPa). Required for clay.
    method : str
        'meyerhof' | 'vesic' | 'skempton' (clay -- alias).
    c : float
        Optional drained cohesion (kPa). Adds c * Nc* contribution.

    Returns
    -------
    dict
        ``{'q_p': kPa, 'Nq*': ..., 'method': ...}``.

    Reference
    ---------
    Meyerhof (1976); Vesic (1977); Das (2014) Ch. 9.
    """
    method = method.lower()
    if Su is not None:
        # Clay -- Skempton's classical.
        check_positive(Su, "Su")
        return {"q_p": float(9 * Su), "Nc_star": 9.0, "method": "skempton_clay"}

    check_range(phi, "phi", 0, 50)
    check_positive(sigma_v_eff, "sigma_v_eff")
    phi_r = np.radians(phi)
    if method == "meyerhof":
        Nq_star = np.exp(np.pi * np.tan(phi_r)) * \
            np.tan(np.pi / 4 + phi_r / 2) ** 2
        # Meyerhof's empirical cap.
        Nq_star = min(Nq_star, 500.0)
    elif method == "vesic":
        # Reduced/refined per Vesic (1977).
        Nq_star = np.exp((np.pi / 2 + phi_r) * np.tan(phi_r)) * \
            np.tan(np.pi / 4 + phi_r / 2)
    else:
        raise ValueError("method must be 'meyerhof', 'vesic', or pass Su for clay.")
    q_p = sigma_v_eff * Nq_star
    if c:
        # Add cohesion contribution (Das Eq. 9.21).
        Nc_star = (Nq_star - 1) / np.tan(phi_r) if phi > 0 else 9.0
        q_p += c * Nc_star
    return {"q_p": float(q_p), "Nq_star": float(Nq_star), "method": method}

pile_skin_friction

Python
pile_skin_friction(Su: float = None, sigma_v_eff: float = None, method: str = 'alpha', phi: float = None, alpha: float = None, beta: float = None, K: float = None, delta: float = None, lambda_: float = None, layer_thicknesses: Sequence[float] = None) -> dict

Unit shaft friction f_s along a pile.

Methods (drained vs total stress):

  • alpha (Tomlinson 1957) -- total stress, for clays: f_s = alpha * Su
  • beta (Burland 1973) -- effective stress, for clays or sands: f_s = beta * sigma'_v where beta = K * tan(delta). Default K = K0 = 1 - sin(phi); delta = phi.
  • lambda (Vijayvergiya & Focht 1972) -- offshore long piles: f_s,avg = lambda * (sigma'_v_avg + 2 * Su_avg)
PARAMETER DESCRIPTION
Su

Undrained shear strength (kPa). Required for alpha.

TYPE: float DEFAULT: None

sigma_v_eff

Effective vertical stress at depth (kPa). Required for beta/lambda.

TYPE: float DEFAULT: None

method

'alpha' | 'beta' | 'lambda'.

TYPE: str DEFAULT: 'alpha'

alpha

Method-specific parameters. Sensible defaults are computed.

TYPE: float DEFAULT: None

beta

Method-specific parameters. Sensible defaults are computed.

TYPE: float DEFAULT: None

K

Method-specific parameters. Sensible defaults are computed.

TYPE: float DEFAULT: None

delta

Method-specific parameters. Sensible defaults are computed.

TYPE: float DEFAULT: None

lambda_

Method-specific parameters. Sensible defaults are computed.

TYPE: float DEFAULT: None

layer_thicknesses

For lambda, weighting of multi-layer averages.

TYPE: sequence DEFAULT: None

Reference

Tomlinson (1957); Burland (1973); Vijayvergiya & Focht (1972).

Source code in geoeq/design/piles.py
Python
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def pile_skin_friction(
    Su: float = None, sigma_v_eff: float = None,
    method: str = "alpha", phi: float = None,
    alpha: float = None, beta: float = None,
    K: float = None, delta: float = None, lambda_: float = None,
    layer_thicknesses: Sequence[float] = None,
) -> dict:
    """Unit shaft friction f_s along a pile.

    Methods (drained vs total stress):

    * **alpha (Tomlinson 1957)** -- total stress, for clays:
        f_s = alpha * Su
    * **beta (Burland 1973)** -- effective stress, for clays or sands:
        f_s = beta * sigma'_v   where beta = K * tan(delta).
        Default K = K0 = 1 - sin(phi); delta = phi.
    * **lambda (Vijayvergiya & Focht 1972)** -- offshore long piles:
        f_s,avg = lambda * (sigma'_v_avg + 2 * Su_avg)

    Parameters
    ----------
    Su : float
        Undrained shear strength (kPa). Required for alpha.
    sigma_v_eff : float
        Effective vertical stress at depth (kPa). Required for beta/lambda.
    method : str
        'alpha' | 'beta' | 'lambda'.
    alpha, beta, K, delta, lambda_ : float
        Method-specific parameters. Sensible defaults are computed.
    layer_thicknesses : sequence, optional
        For ``lambda``, weighting of multi-layer averages.

    Reference
    ---------
    Tomlinson (1957); Burland (1973); Vijayvergiya & Focht (1972).
    """
    method = method.lower()
    if method == "alpha":
        check_positive(Su, "Su")
        if alpha is None:
            # API RP 2A (1993) tabulation; simplified.
            if Su <= 25:
                alpha = 1.0
            elif Su <= 75:
                alpha = 1.0 - 0.5 * (Su - 25) / 50  # linear 1.0 -> 0.5
            elif Su <= 175:
                alpha = 0.5 - 0.0014 * (Su - 75)    # linear 0.5 -> ~0.35
            else:
                alpha = 0.35
        return {"f_s": float(alpha * Su), "alpha": float(alpha),
                "method": "alpha"}
    if method == "beta":
        check_positive(sigma_v_eff, "sigma_v_eff")
        if beta is None:
            if K is None:
                if phi is None:
                    raise ValueError("beta method needs phi (or beta/K).")
                K = 1 - np.sin(np.radians(phi))
            if delta is None:
                if phi is None:
                    raise ValueError("beta method needs phi (or delta).")
                delta = phi
            beta = K * np.tan(np.radians(delta))
        return {"f_s": float(beta * sigma_v_eff), "beta": float(beta),
                "method": "beta"}
    if method == "lambda":
        check_positive(sigma_v_eff, "sigma_v_eff")
        check_positive(Su, "Su")
        if lambda_ is None:
            # Vijayvergiya & Focht (1972) chart, approximated.
            # lambda decreases with embedment; rough default 0.3.
            lambda_ = 0.3
        f_s = lambda_ * (sigma_v_eff + 2 * Su)
        return {"f_s": float(f_s), "lambda": float(lambda_),
                "method": "lambda"}
    raise ValueError("method must be 'alpha', 'beta', or 'lambda'.")

pile_capacity

Python
pile_capacity(D: float, L: float, q_p: float, f_s: float, area_base: float = None, perimeter: float = None, FS: float = 3.0) -> dict

Axial capacity Q_ult = Q_p + Q_s = q_p * A_p + f_s * As.

PARAMETER DESCRIPTION
D

Pile diameter (m).

TYPE: float

L

Pile length (embedded depth, m).

TYPE: float

q_p

Unit tip resistance (kPa) -- from pile_end_bearing.

TYPE: float

f_s

Unit shaft friction (kPa) -- from pile_skin_friction. Pass the depth-averaged value, or pre-multiply by length.

TYPE: float

area_base

Base area (m^2). Default pi*D^2/4 (closed-ended round pile).

TYPE: float DEFAULT: None

perimeter

Shaft perimeter (m). Default pi*D.

TYPE: float DEFAULT: None

FS

Factor of safety on Q_ult.

TYPE: float DEFAULT: 3.0

Reference

Das (2014) Eq. 9.20.

Source code in geoeq/design/piles.py
Python
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
def pile_capacity(
    D: float, L: float,
    q_p: float, f_s: float, area_base: float = None,
    perimeter: float = None, FS: float = 3.0,
) -> dict:
    """Axial capacity Q_ult = Q_p + Q_s = q_p * A_p + f_s * As.

    Parameters
    ----------
    D : float
        Pile diameter (m).
    L : float
        Pile length (embedded depth, m).
    q_p : float
        Unit tip resistance (kPa) -- from ``pile_end_bearing``.
    f_s : float
        Unit shaft friction (kPa) -- from ``pile_skin_friction``.
        Pass the depth-averaged value, or pre-multiply by length.
    area_base : float
        Base area (m^2). Default pi*D^2/4 (closed-ended round pile).
    perimeter : float
        Shaft perimeter (m). Default pi*D.
    FS : float
        Factor of safety on Q_ult.

    Reference
    ---------
    Das (2014) Eq. 9.20.
    """
    check_positive(D, "D")
    check_positive(L, "L")
    if area_base is None:
        area_base = np.pi * D ** 2 / 4.0
    if perimeter is None:
        perimeter = np.pi * D
    Q_p = q_p * area_base
    Q_s = f_s * perimeter * L
    Q_ult = Q_p + Q_s
    Q_all = Q_ult / FS
    return {
        "Q_p": float(Q_p), "Q_s": float(Q_s),
        "Q_ult": float(Q_ult), "Q_all": float(Q_all),
        "FS": float(FS),
    }

pile_group_efficiency

Python
pile_group_efficiency(n: int, m: int, D: float, s: float, method: str = 'converse_labarre') -> dict

Group efficiency factor eta for an n x m pile group.

Methods: * Converse-Labarre: eta = 1 - [theta * ((n-1)m + (m-1)n)] / (90 * n * m) with theta = atan(D/s) in degrees. * Feld: each pile loses 1/16 capacity for each adjacent pile.

Reference

Converse & Labarre (1947); Feld (1943); Das (2014) Eq. 9.74.

Source code in geoeq/design/piles.py
Python
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
def pile_group_efficiency(
    n: int, m: int, D: float, s: float, method: str = "converse_labarre",
) -> dict:
    """Group efficiency factor eta for an n x m pile group.

    Methods:
    * **Converse-Labarre**:
        eta = 1 - [theta * ((n-1)*m + (m-1)*n)] / (90 * n * m)
        with theta = atan(D/s) in degrees.
    * **Feld**: each pile loses 1/16 capacity for each adjacent pile.

    Reference
    ---------
    Converse & Labarre (1947); Feld (1943); Das (2014) Eq. 9.74.
    """
    if n < 1 or m < 1:
        raise ValueError("n, m must be >= 1.")
    check_positive(D, "D")
    check_positive(s, "s")
    method = method.lower()
    if method in ("converse_labarre", "cl"):
        theta = np.degrees(np.arctan(D / s))
        eta = 1 - theta * ((n - 1) * m + (m - 1) * n) / (90 * n * m)
    elif method == "feld":
        # Per pile: 1 - n_adj * 1/16. Average over the group.
        # Corner: 3 adj; edge: 5; interior: 8.
        total = 0.0
        for i in range(n):
            for j in range(m):
                adj = 0
                for di in (-1, 0, 1):
                    for dj in (-1, 0, 1):
                        if (di, dj) == (0, 0):
                            continue
                        if 0 <= i + di < n and 0 <= j + dj < m:
                            adj += 1
                total += max(0, 1 - adj / 16)
        eta = total / (n * m)
    else:
        raise ValueError("method must be 'converse_labarre' or 'feld'.")
    return {"eta": float(eta), "method": method,
            "n_piles": n * m, "spacing_ratio": s / D}

pile_settlement

Python
pile_settlement(Q_w: float, Q_p: float, Q_s: float, D: float, L: float, Es: float, Ep: float = 25000000.0, mu_s: float = 0.3, Cp: float = 0.03, Cs: float = None, qp_ult: float = None) -> dict

Vesic's three-component settlement of a single pile.

Text Only
s = s1 + s2 + s3
PARAMETER DESCRIPTION
Q_w

Working axial load (kN).

TYPE: float

Q_p

Tip load at working condition (kN).

TYPE: float

Q_s

Shaft load at working condition (kN).

TYPE: float

D

Pile diameter and length (m).

TYPE: float

L

Pile diameter and length (m).

TYPE: float

Es

Soil elastic modulus (kPa).

TYPE: float

Ep

Pile material modulus (kPa). Default 25 GPa (concrete).

TYPE: float DEFAULT: 25000000.0

mu_s

Soil Poisson's ratio.

TYPE: float DEFAULT: 0.3

Cp

Empirical tip-settlement coefficient (Das Table 9.5). Default 0.03.

TYPE: float DEFAULT: 0.03

Cs

Shaft-settlement coefficient. Default 0.93 + 0.16 sqrt(L/D) * Cp.

TYPE: float DEFAULT: None

qp_ult

Ultimate tip resistance (kPa). If given, used in s2 calc.

TYPE: float DEFAULT: None

Reference

Vesic (1977); Das (2014) Eq. 9.83-9.86.

Source code in geoeq/design/piles.py
Python
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
def pile_settlement(
    Q_w: float, Q_p: float, Q_s: float, D: float, L: float, Es: float,
    Ep: float = 25e6, mu_s: float = 0.3,
    Cp: float = 0.03, Cs: float = None, qp_ult: float = None,
) -> dict:
    """Vesic's three-component settlement of a single pile.

        s = s1 + s2 + s3

    Parameters
    ----------
    Q_w : float
        Working axial load (kN).
    Q_p : float
        Tip load at working condition (kN).
    Q_s : float
        Shaft load at working condition (kN).
    D, L : float
        Pile diameter and length (m).
    Es : float
        Soil elastic modulus (kPa).
    Ep : float
        Pile material modulus (kPa). Default 25 GPa (concrete).
    mu_s : float
        Soil Poisson's ratio.
    Cp : float
        Empirical tip-settlement coefficient (Das Table 9.5). Default 0.03.
    Cs : float, optional
        Shaft-settlement coefficient. Default 0.93 + 0.16 sqrt(L/D) * Cp.
    qp_ult : float
        Ultimate tip resistance (kPa). If given, used in s2 calc.

    Reference
    ---------
    Vesic (1977); Das (2014) Eq. 9.83-9.86.
    """
    check_positive(D, "D")
    check_positive(L, "L")
    check_positive(Es, "Es")
    check_positive(Ep, "Ep")
    # s1: elastic compression of pile shaft.
    # s1 = (Q_p + xi * Q_s) * L / (A_p * E_p)   with xi ~ 0.5 to 0.67 (avg).
    xi = 0.5
    A_p = np.pi * D ** 2 / 4
    s1 = (Q_p + xi * Q_s) * L / (A_p * Ep)
    # s2: tip settlement caused by Q_p.
    if qp_ult is None:
        qp_ult = (Q_p / A_p) * 3  # rough -- mobilization at FS=3
    s2 = (Q_p / A_p) * D * (1 - mu_s ** 2) / Es * Cp  # Vesic semi-empirical
    if Cs is None:
        Cs = (0.93 + 0.16 * np.sqrt(L / D)) * Cp
    # s3: shaft skin friction transferred load.
    s3 = (Q_s / (np.pi * D * L)) * D * (1 - mu_s ** 2) / Es * Cs
    s_total = s1 + s2 + s3
    return {
        "s1_elastic": float(s1), "s2_tip": float(s2), "s3_shaft": float(s3),
        "s_total": float(s_total),
    }

Slope stability

infinite_slope

Python
infinite_slope(phi: float, beta: float, c: float = 0.0, gamma: float = 18.0, H: float = 1.0, seepage: bool = False, gamma_sat: float = None) -> dict

Factor of safety of an infinite slope.

Cohesionless, dry/no-seepage: FS = tan(phi) / tan(beta)

Cohesive (c-phi) without seepage: FS = c / (gamma * H * cos(beta)^2 * tan(beta)) + tan(phi) / tan(beta)

Cohesionless with full seepage parallel to slope: FS = (gamma' / gamma_sat) * tan(phi) / tan(beta)

PARAMETER DESCRIPTION
phi

Friction angle (degrees).

TYPE: float

beta

Slope inclination (degrees).

TYPE: float

c

Cohesion (kPa). Default 0.

TYPE: float DEFAULT: 0.0

gamma

Bulk unit weight (kN/m^3).

TYPE: float DEFAULT: 18.0

H

Depth of failure plane below surface (m). Only matters if c > 0.

TYPE: float DEFAULT: 1.0

seepage

If True, assumes seepage parallel to slope to the surface.

TYPE: bool DEFAULT: False

gamma_sat

Saturated unit weight (kN/m^3). Required if seepage=True.

TYPE: float DEFAULT: None

RETURNS DESCRIPTION
dict

{'FS': ..., 'gamma_eff': ...}.
Reference

Das (2010) Eq. 13.4-13.10.

Source code in geoeq/design/slopes.py
Python
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
def infinite_slope(
    phi: float, beta: float, c: float = 0.0, gamma: float = 18.0,
    H: float = 1.0, seepage: bool = False, gamma_sat: float = None,
) -> dict:
    """Factor of safety of an infinite slope.

    Cohesionless, dry/no-seepage:
        FS = tan(phi) / tan(beta)

    Cohesive (c-phi) without seepage:
        FS = c / (gamma * H * cos(beta)^2 * tan(beta))
             + tan(phi) / tan(beta)

    Cohesionless with full seepage parallel to slope:
        FS = (gamma' / gamma_sat) * tan(phi) / tan(beta)

    Parameters
    ----------
    phi : float
        Friction angle (degrees).
    beta : float
        Slope inclination (degrees).
    c : float
        Cohesion (kPa). Default 0.
    gamma : float
        Bulk unit weight (kN/m^3).
    H : float
        Depth of failure plane below surface (m). Only matters if c > 0.
    seepage : bool
        If True, assumes seepage parallel to slope to the surface.
    gamma_sat : float
        Saturated unit weight (kN/m^3). Required if seepage=True.

    Returns
    -------
    dict
        ``{'FS': ..., 'gamma_eff': ...}``.

    Reference
    ---------
    Das (2010) Eq. 13.4-13.10.
    """
    check_range(phi, "phi", 0, 50)
    check_range(beta, "beta", 0.01, 89.99)
    check_non_negative(c, "c")
    check_positive(gamma, "gamma")
    beta_r = np.radians(beta)
    phi_r = np.radians(phi)

    if seepage:
        if gamma_sat is None:
            raise ValueError("gamma_sat required when seepage=True.")
        gamma_eff = gamma_sat - GAMMA_WATER
        FS = (gamma_eff / gamma_sat) * np.tan(phi_r) / np.tan(beta_r)
        if c > 0:
            FS += c / (gamma_sat * H * np.cos(beta_r) ** 2 * np.tan(beta_r))
    else:
        FS = np.tan(phi_r) / np.tan(beta_r)
        if c > 0:
            FS += c / (gamma * H * np.cos(beta_r) ** 2 * np.tan(beta_r))
    return {"FS": float(FS), "seepage": bool(seepage)}

culmann

Python
culmann(c: float, phi: float, gamma: float, beta: float) -> dict

Culmann's critical height H_cr for a planar-failure slope.

Text Only
H_cr = (4 c sin(beta) cos(phi)) / (gamma * (1 - cos(beta - phi)))
PARAMETER DESCRIPTION
c

Cohesion (kPa).

TYPE: float

phi

Friction angle (degrees).

TYPE: float

gamma

Bulk unit weight (kN/m^3).

TYPE: float

beta

Slope angle (degrees), > phi.

TYPE: float

Reference

Culmann (1875); Das (2010) Eq. 13.13.

Source code in geoeq/design/slopes.py
Python
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
def culmann(c: float, phi: float, gamma: float, beta: float) -> dict:
    """Culmann's critical height H_cr for a planar-failure slope.

        H_cr = (4 c sin(beta) cos(phi)) / (gamma * (1 - cos(beta - phi)))

    Parameters
    ----------
    c : float
        Cohesion (kPa).
    phi : float
        Friction angle (degrees).
    gamma : float
        Bulk unit weight (kN/m^3).
    beta : float
        Slope angle (degrees), > phi.

    Reference
    ---------
    Culmann (1875); Das (2010) Eq. 13.13.
    """
    check_positive(c, "c")
    check_range(phi, "phi", 0, 50)
    check_positive(gamma, "gamma")
    check_range(beta, "beta", 0.01, 89.99)
    if beta <= phi:
        raise ValueError("Culmann requires beta > phi.")
    beta_r = np.radians(beta)
    phi_r = np.radians(phi)
    H_cr = (4 * c * np.sin(beta_r) * np.cos(phi_r)) / \
        (gamma * (1 - np.cos(beta_r - phi_r)))
    return {"H_cr": float(H_cr), "method": "culmann"}

taylor_stability

Python
taylor_stability(phi: float, c: float, gamma: float, H: float, beta: float) -> dict

Factor of safety by Taylor's stability number m.

Text Only
m = c / (gamma * H_cr)         (Taylor 1937)
FS = c / (m * gamma * H)

Uses Taylor's chart (interpolated table for phi=0..25, beta=15..90).

Reference

Taylor (1937, 1948); Das (2010) Fig. 13.10.

Source code in geoeq/design/slopes.py
Python
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def taylor_stability(phi: float, c: float, gamma: float, H: float,
                     beta: float) -> dict:
    """Factor of safety by Taylor's stability number m.

        m = c / (gamma * H_cr)         (Taylor 1937)
        FS = c / (m * gamma * H)

    Uses Taylor's chart (interpolated table for phi=0..25, beta=15..90).

    Reference
    ---------
    Taylor (1937, 1948); Das (2010) Fig. 13.10.
    """
    check_non_negative(c, "c")
    check_range(phi, "phi", 0, 30)
    check_positive(gamma, "gamma")
    check_positive(H, "H")
    check_range(beta, "beta", 15, 90)

    # Taylor's stability number m for phi=0 (Das Fig 13.10b approximation):
    # For phi=0: m depends on beta and depth factor D=H'/H.
    # For phi > 0: m decreases. Simplified curve-fit (Das Table 13.1):
    # The standard "toe circle" chart for phi = 0..25, beta = 15..90 is hard
    # to reproduce closed-form; use a Janbu-style empirical fit:
    #   m approx = 0.181 - 0.0028 phi - 0.0006 beta + 0.000007 phi * beta
    # This matches Taylor's chart to ~0.01 in the practical range.
    m = max(0.01,
            0.181 - 0.0028 * phi - 0.0006 * beta + 7e-6 * phi * beta)
    FS = c / (m * gamma * H)
    return {"m": float(m), "FS": float(FS), "method": "taylor"}

bishop

Python
bishop(slices: Sequence[dict], R: float = None, max_iter: int = 50, tol: float = 0.0001) -> dict

Bishop's simplified method (1955) for circular slope failure.

Each slice is a dict with keys: b -- slice width (m) h -- slice height (m) alpha -- base inclination (degrees, +ve when uphill) c -- cohesion on slice base (kPa) phi -- friction angle (degrees) gamma -- unit weight (kN/m^3) u -- pore pressure at slice base (kPa), default 0

Iterates: FS_{k+1} = sum_i [ (c_i b_i + (W_i - u_i b_i) tan phi_i) / m_alpha_i ] / sum_i ( W_i sin alpha_i )

with m_alpha_i = cos alpha_i + sin alpha_i tan phi_i / FS_k.

Reference

Bishop (1955); Das (2010) Eq. 13.50.

Source code in geoeq/design/slopes.py
Python
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
def bishop(
    slices: Sequence[dict], R: float = None,
    max_iter: int = 50, tol: float = 1e-4,
) -> dict:
    """Bishop's simplified method (1955) for circular slope failure.

    Each ``slice`` is a dict with keys:
        b      -- slice width (m)
        h      -- slice height (m)
        alpha  -- base inclination (degrees, +ve when uphill)
        c      -- cohesion on slice base (kPa)
        phi    -- friction angle (degrees)
        gamma  -- unit weight (kN/m^3)
        u      -- pore pressure at slice base (kPa), default 0

    Iterates:
        FS_{k+1} = sum_i [ (c_i b_i + (W_i - u_i b_i) tan phi_i) / m_alpha_i ]
                   / sum_i ( W_i sin alpha_i )

    with m_alpha_i = cos alpha_i + sin alpha_i tan phi_i / FS_k.

    Reference
    ---------
    Bishop (1955); Das (2010) Eq. 13.50.
    """
    if not slices:
        raise ValueError("At least one slice required.")
    # Normalize the slice fields.
    def W(s):
        return s.get("gamma", 18.0) * s["b"] * s["h"]
    def alpha(s):
        return np.radians(s["alpha"])
    def phi_rad(s):
        return np.radians(s.get("phi", 0))
    def c_(s):
        return s.get("c", 0.0)
    def u_(s):
        return s.get("u", 0.0)
    def b_(s):
        return s["b"]

    driving = sum(W(s) * np.sin(alpha(s)) for s in slices)
    if abs(driving) < 1e-9:
        raise ValueError("Sum of W*sin(alpha) is zero -- check slice geometry.")

    FS = 1.0  # initial guess
    for it in range(max_iter):
        num = 0.0
        for s in slices:
            Wi = W(s)
            m_alpha = np.cos(alpha(s)) + \
                np.sin(alpha(s)) * np.tan(phi_rad(s)) / FS
            if m_alpha <= 0:
                # ill-conditioned; nudge FS upward.
                m_alpha = 0.01
            num += (c_(s) * b_(s)
                    + (Wi - u_(s) * b_(s)) * np.tan(phi_rad(s))) / m_alpha
        FS_new = num / driving
        if abs(FS_new - FS) < tol:
            FS = FS_new
            return {"FS": float(FS), "iterations": it + 1,
                    "converged": True, "method": "bishop_simplified"}
        FS = FS_new
    return {"FS": float(FS), "iterations": max_iter,
            "converged": False, "method": "bishop_simplified"}

taylor_chart_plot

Python
taylor_chart_plot(phi_values=(0, 5, 10, 15, 20, 25), beta_range=None, ax=None, save_as: str = None)

Taylor's stability number chart: m vs slope angle beta for a series of phi values.

Source code in geoeq/design/plots.py
Python
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
def taylor_chart_plot(phi_values=(0, 5, 10, 15, 20, 25),
                     beta_range=None, ax=None, save_as: str = None):
    """Taylor's stability number chart: m vs slope angle beta for
    a series of phi values."""
    import matplotlib.pyplot as plt
    if beta_range is None:
        beta_range = np.linspace(15, 90, 30)
    if ax is None:
        fig, ax = plt.subplots(figsize=(7, 5))
    else:
        fig = ax.figure
    for phi in phi_values:
        m_values = []
        for beta in beta_range:
            # We bypass FS by reading the m value directly.
            res = taylor_stability(phi=phi, c=1, gamma=1, H=1, beta=beta)
            m_values.append(res["m"])
        ax.plot(beta_range, m_values, label=f"$\\phi = {phi}°$",
                linewidth=1.7)
    ax.set_xlabel(r"Slope angle $\beta$ (°)")
    ax.set_ylabel(r"Stability number $m = c / (\gamma H_{cr})$")
    ax.set_title("Taylor's stability number (toe circles)")
    ax.grid(alpha=0.3)
    ax.legend(loc="upper right", fontsize=9)
    if save_as:
        fig.savefig(save_as, dpi=300, bbox_inches="tight")
    return fig