Examples

Example 01: Create a Soil Profile

"""This example demonstrates how to create a soil profile. The groundwater is at the ground surface."""

from pyvstress import SoilLayer, SoilProfile


# Create layers using the SoilLayer object, with the required parameters
# SoilLayer(layer_thickness, gamma_bulk, gamma_sat)
# layer_thickness and gamma_bulk are required soil parameters for each layer
# if a value for gamma_sat is not provided then the value for gamma_bulk is used
layer1 = SoilLayer(5, 110, 120)
layer2 = SoilLayer(7, 115, 125)
layer3 = SoilLayer(4, 115)
layer4 = SoilLayer(1, 112, 120)
layer5 = SoilLayer(10, 125)
layers = [layer1, layer2, layer3, layer4, layer5]

# create a soil profile by adding soil layers to the SoilProfile object
# soil layers can only be added as a list during initialization as shown below.
# Also required is the depth to groundwater, positive depth is below ground surface
# and negative depth above the ground surface indicate standing water.
# The default unit weight of water is 62.4 and thus gammaw can be omitted.
# No units is assumed in the calculations, so the unitweight of water determines the units
# of calculation
soilprofile = SoilProfile(layers, zw=0.0, gammaw=62.4)

print("The soil layers:")
print(
    f"{'ztop [ft]':>10s}, {'zbot [ft]':>10s}, {'layer_thickness':>20s}, {'gamma_bulk':>15s}, {'gamma_sat':>15s}"
)
for layer in soilprofile:
    print(
        f"{layer.ztop:>10.2f}, {layer.zbot:>10.2f}, {layer.thickness:>20.2f}, {layer.gamma_bulk:>15.2f}, {layer.gamma_sat:>15.2f}"
    )

The above code should return:

The soil layers:
ztop [ft],  zbot [ft],      layer_thickness,      gamma_bulk,       gamma_sat
     0.00,       5.00,                 5.00,          110.00,          120.00
     5.00,      12.00,                 7.00,          115.00,          125.00
    12.00,      16.00,                 4.00,          115.00,          115.00
    16.00,      17.00,                 1.00,          112.00,          120.00
    17.00,      27.00,                10.00,          125.00,          125.00

Example 02: Soil profile with watertable inbetween layer boundary

"""This example demonstrates the splitting of the layer at the depth of the groundwater.""""

from pyvstress import SoilLayer, SoilProfile


# create soil layers
layer1 = SoilLayer(5, 110, 120)
layer2 = SoilLayer(7, 115, 125)
layer3 = SoilLayer(4, 115)
layer4 = SoilLayer(1, 112, 120)
layer5 = SoilLayer(10, 125)
layers = [layer1, layer2, layer3, layer4, layer5]

# create soil profile with groundwater within a soil layer
soilprofile = SoilProfile(layers, zw=10)

print("The soil layers:")
print(
    f"{'ztop [ft]':>10s}, {'zbot [ft]':>10s}, {'layer_thickness':>20s}, {'gamma_bulk':>15s}, {'gamma_sat':>15s}"
)
for layer in soilprofile:
    print(
        f"{layer.ztop:>10.2f}, {layer.zbot:>10.2f}, {layer.thickness:>20.2f}, {layer.gamma_bulk:>15.2f}, {layer.gamma_sat:>15.2f}"
    )

The above code should return:

The soil layers:
ztop [ft],  zbot [ft],      layer_thickness,      gamma_bulk,       gamma_sat
    0.00,       5.00,                 5.00,          110.00,          120.00
    5.00,      10.00,                 5.00,          115.00,          125.00
   10.00,      12.00,                 2.00,          115.00,          125.00
   12.00,      16.00,                 4.00,          115.00,          115.00
   16.00,      17.00,                 1.00,          112.00,          120.00
   17.00,      27.00,                10.00,          125.00,          125.00

Example 03: Compute vertical stress profile

"""This example demonstrates vertical stresss profile. The vertical stresses are computed at the top of the soil profile, at each soil layer boundary, at the groundwater table if it is not one of the layer boundaries, and at the bottom of the soil profile. The stress is linear between these points.
file: examples/example_03.py""""

from pyvstress import SoilLayer, SoilProfile
import pyvstress.utility_functions as utf


# Create layers
layer1 = SoilLayer(5, 100, 110, phi=26, c=100)
layer2 = SoilLayer(7, 115, 125, phi=28, c=100)
layer3 = SoilLayer(4, 115, 122, phi=30)
layer4 = SoilLayer(1, 110, 115, phi=25, c=25)
layer5 = SoilLayer(10, 125, 135, phi=35)
layers = [layer1, layer2, layer3, layer4, layer5]

# create soil profile, with watertable at 7 feet from the ground surface
soilprofile = SoilProfile(layers, zw=7.0, gammaw=62.4)

# get the points for creating vertical stresse profile
# Point objects are returned
pts = utf.stress_profile_points(soilprofile)

# compute vertical stresses at the above points, the stresses computed are stored in the pts
# and can be retrieved with pt.total_stress, pt.pore_pressure, pt.effective_stress
soilprofile.vertical_stresses(pts)
# results = []
# for pt in pts:
#     results.append([pt.z, pt.total_stress, pt.pore_pressure, pt.effective_stress])

print("Vertical stresses:")
print(
    f"{'z':>5s}, {'Total stress':>15s}, {'Pore press.':>15s}, {'Eff. stress':>15s}"
)
for pt in pts:
    print(
        f"{pt.z:>5.2f}, {pt.total_stress:>15.2f}, {pt.pore_pressure:>15.2f}, {pt.effective_stress:>15.2f}"
    )

The above code should return.

Vertical stresses:
z    ,    Total stress,     Pore press.,     Eff. stress
 0.00,            0.00,            0.00,            0.00
 5.00,          500.00,            0.00,          500.00
 7.00,          730.00,            0.00,          730.00
12.00,         1355.00,          312.00,         1043.00
16.00,         1843.00,          561.60,         1281.40
17.00,         1958.00,          624.00,         1334.00
27.00,         3308.00,         1248.00,         2060.00

Example 04: Compute lateral earth pressure profile

"""This example demonstrates lateral earth pressure profile. The lateral earth pressures are computed at the top of the profile, two calculations at each of the boundary layers, one for the above layer and one for the bottom layer, at the groundwater depth if the groundwater is not one of the soil layer boundaries, at at the bottom of the soil profile. Two calculations are required at layer boundaries because soil parameters are different for the above and bottom layer at the boundary. The stress is linear between these points.
file: examples/example_04.py"""

# Create layers with optional arguments for lateral earth pressure calculations
layer1 = SoilLayer(5, 100, 110, phi=26, c=100)
layer2 = SoilLayer(7, 115, 125, phi=28, c=100)
layer3 = SoilLayer(4, 115, 122, phi=30, c=0)
layer4 = SoilLayer(1, 110, 115, phi=25, c=25)
layer5 = SoilLayer(10, 125, 135, phi=35, c=0)
layers = [layer1, layer2, layer3, layer4, layer5]

# create soil profile, with watertable at 7 feet from the ground surface
soilprofile = SoilProfile(layers, zw=7.0, gammaw=62.4)

# use the utility function to compute lateral earth pressures for utility
vertstresses, latstresses = utf.profile_lateral_earth_pressures(soilprofile)

print("Vertical stresses:")
print(
    f"{'z':>5s}, {'Total stress':>15s}, {'Pore press.':>15s}, {'Eff. stress':>15s}"
)
for vstress in vertstresses:
    print(
        f"{vstress[0]:>5.2f}, {vstress[1]:>15.2f}, {vstress[2]:>15.2f}, {vstress[3]:>15.2f}"
    )

print("Lateral earth pressures:")
print(
    f"{'z':>5s}, {'Eff. Stress':>15s}, {'Active Pressure':>15s}, {'Passive Pressure':>15s}"
)
for vertstress, latstress in zip(vertstresses, latstresses):
    print(
        f"{latstress[0]:>5.2f}, {vertstress[3]:>15.2f}, {latstress[1]:>15.2f}, {latstress[2]:>15.2f}"
    )

The above code should return. Since the primary interest is lateral earth pressure profile there are two points at each internal layer boundary.

Vertical stresses:
    z,    Total stress,     Pore press.,     Eff. stress
 0.00,            0.00,            0.00,            0.00
 5.00,          500.00,            0.00,          500.00
 5.00,          500.00,            0.00,          500.00
 7.00,          730.00,            0.00,          730.00
12.00,         1355.00,          312.00,         1043.00
12.00,         1355.00,          312.00,         1043.00
16.00,         1843.00,          561.60,         1281.40
16.00,         1843.00,          561.60,         1281.40
17.00,         1958.00,          624.00,         1334.00
17.00,         1958.00,          624.00,         1334.00
27.00,         3308.00,         1248.00,         2060.00

Lateral earth pressures:
    z,     Eff. stress, Active Pressure, Passive Pressure
 0.00,            0.00,         -124.97,          320.07
 5.00,          500.00,           70.26,         1600.60
 5.00,          500.00,           60.34,         1717.77
 7.00,          730.00,          143.38,         2354.83
12.00,         1043.00,          256.39,         3221.78
12.00,         1043.00,          347.67,         3129.00
16.00,         1281.40,          427.13,         3844.20
16.00,         1281.40,          488.21,         3235.74
17.00,         1334.00,          509.56,         3365.34
17.00,         1334.00,          361.50,         4922.69
27.00,         2060.00,          558.24,         7601.76

Example 05: Compute vertical stresses and lateral earth pressures without the use of utility functions

In this example we are going to demonstrates that this simple package is mainly intented to compute vertical stresses at discrete points and to return soil parameters at those points. If the soil parameter falls on a layer boundary it is assumed to belong to the above layer. We can then use the vertical stresses and the soil parameters to perform subsequent further analyses. This is example_09.py in the examples folder. Check example_10.py for an alternative method.

Steps that we will follow in this example are:

  1. Create a soil profile. Supply all soil parameters required for the analysis during layer initialization. e.g. if we will use c and phi for computing active and passive earth pressures, then the soil layers have to be initialized with c and phi in addition to the three required parameters, layer_thickness, bulk_unit_weight, and saturated_unit_weight. If we want to use the c parameter for some layers and not for other layers do not omit c in the layers we do not want to use c, assign c = 0. When extracting soil layer parameters, all soil layers must have those parameters.

# Create layers with optional arguments for lateral earth pressure calculations
layer1 = SoilLayer(5, 100, 110, phi=26, c=50)
layer2 = SoilLayer(10, 110, 120, phi=30, c=0)
layer3 = SoilLayer(7, 120, 130, phi=34, c=0)
layers = [layer1, layer2, layer3]

# create soil profile, with watertable at 5 feet from the ground surface and surcharge
soilprofile = SoilProfile(layers, zw=6.0, q=250, gammaw=62.4)
  1. Create class::Point instances at the required depths by initializing the class::Point. In the above examples we only required points at the layer boundaries because we were interested in drawing stress profiles and the slope of the vertical and lateral stresses change (also the value) at the layer boundaries. However, for this example we will only consider points inbetween the layer boundaries to demonstrate the use of this package. If the analysis being done requires lateral stresses then use stress_profile_points(soilprofile, analysis_type=”lateral”) for two-points at layer boundaries and then add intermediate points.

  2. Determine vertical stresses at those points.

# depth at which to calculate vertical and lateral stresses
# again do not choose layer boundaries for lateral earth pressures
# but if you are only interested in vertical stresses then it is alright to choose layer
# boundaries depths
# the second layer will be subdivided into two layers at the groundwater level but that ok
# because both the layers obove and below that boundary have the same parameters so one pont
# at that boundary is alright
zs = [0, 3.0, 6.0, 10.0, 16.0, 20.0, 22.0]

# convert the depth to points
pts = [Point(z) for z in zs]

# compute vertical total stress, porewater pressure, and effective stresses
soilprofile.vertical_stresses(pts)

# if you are only interested in vertical stresses then you can view them
# for pt in pts:
#     print(pt.z, pt.total_stress, pt.pore_pressure, pt.effective_stress)
# and if you want it as an numpy aarray
vstress = np.asarray(
    [[pt.z, pt.total_stress, pt.pore_pressure, pt.effective_stress] for pt in pts]
)
  1. Extract layer properties at those points. Since, in this example, c and phi are required for lateral earth pressure calculations we will extract those for each layer. First of all, we need to make sure we have supplied these parameters during layer initialization for all the layers. If a parameter is zero, e.g. c = 0, then initialize the layer as c = 0, rather than omitting the parameter. The parameter names must match in every layer. However, the points do not know in which layer they are as yet, so first we need to assign the layer index of the layer the points are in.

soilprofile.assign_layer_index(pts)
keylist = ["c", "phi"]
params = soilprofile.get_params(pts, keylist)
# params is a list, convert into numpy array, the shape of params will be (7,2), seven rowws
# for the seven points and two columns for the two parameters
params = np.asarray(params)
# Let us check if the correct parameters are extracted for the layers
for p in params:
    print(p)

The above should return:

[50 26]
[50 26]
[ 0 30]
[ 0 30]
[ 0 34]
[ 0 34]
[ 0 34]

The first two points are from layer_1, the second two from layer_3 and the last three from layer_3. So it is all good.

  1. That is it, that was the purpose of this package. Now use custom functions to perfom further analysis. Here we will compute lateral earth pressures. The functions are in the example_09.py.

latstress = np.zeros((vstress.shape[0], 3))
# the first column are the depths of the points
latstress[:, 0] = vstress[:, 0]
# the second column is the active earth pressure values at those points
latstress[:, 1] = calc_active_earth_pressure(
    vstress[:, 3], params[:, 0], params[:, 1]
)
# the third column is the passive earth pressure values at those points
latstress[:, 2] = calc_passive_earth_pressure(
    vstress[:, 3], params[:, 0], params[:, 1]
)

print("Vertical stresses:")
print(
    f"{'z':>5s}, {'Total stress':>15s}, {'Pore press.':>15s}, {'Eff. stress':>15s}"
)
for vstress in vstresses:
    print(
        f"{vstress[0]:>5.2f}, {vstress[1]:>15.2f}, {vstress[2]:>15.2f}, {vstress[3]:>15.2f}"
    )

print("Lateral stresses:")
print(
    f"{'z':>5s}, {'Active Pressure':>15s}, {'Passive Pressure':>15s}"
)
for latstress in latstresses:
    print(
        f"{latstress[0]:>5.2f}, {latstress[1]:>15.2f}, {latstress[2]:>15.2f}"
    )

The above code should return:

Vertical stresses:
    z,    Total stress,     Pore press.,     Eff. stress
 0.00,          250.00,            0.00,          250.00
 3.00,          550.00,            0.00,          550.00
 6.00,          860.00,            0.00,          860.00
10.00,         1340.00,          249.60,         1090.40
16.00,         2070.00,          624.00,         1446.00
20.00,         2590.00,          873.60,         1716.40
22.00,         2850.00,          998.40,         1851.60

Lateral stresses:
    z, Active Pressure, Passive Pressure
 0.00,           14.19,         2030.03
 3.00,           53.93,         4294.62
 6.00,          800.00,          920.00
10.00,         1030.40,         1150.40
16.00,         1378.00,         1514.00
20.00,         1648.40,         1784.40
22.00,         1783.60,         1919.60