"""
This module provides classes and functions for calculating the capacity of soils using various
models for a given pipe and soil configuration.
**Features:**
- The `PSI` class implements soil capacity calculations for different pipe and soil properties,
supporting vectorized and model-based approaches.
- Designed for use in subsea pipeline and riser engineering, but general enough for any
pipe-soil interaction analysis.
- All calculations are vectorized using NumPy for efficiency and flexibility.
.. raw:: html
<hr style="height:6px; background-color:#888; border:none; margin:1.5em 0;" />
"""
import numpy as np
[docs]
class PSI: # pylint: disable=too-many-arguments
"""
Class for calculating the capacity of soils using various models
for a given pipe and soil configuration.
Parameters
----------
total_outer_diameter : float, optional
Total outer diameter of the pipe, including coating.
surface_roughness : str, optional
Surface roughness of the pipe ('Smooth' or 'Rough').
undrained_shear_strength_depth_array : array_like, optional
Per-pipe depth values, in m, defining the undrained shear strength profile.
undrained_shear_strength_value_array : array_like, optional
Per-pipe undrained shear strength values, in Pa, matching
``undrained_shear_strength_depth_array``.
submerged_unit_weight : float, optional
Submerged unit weight of the soil, gamma', in N/m3.
"""
def __init__(
self,
*,
total_outer_diameter=0.0,
surface_roughness=None,
undrained_shear_strength_depth_array=None,
undrained_shear_strength_value_array=None,
submerged_unit_weight=0.0
):
"""
Initialize a PipeSoilInteraction object with pipe and soil properties.
"""
self.total_outer_diameter = np.asarray(total_outer_diameter, dtype = float)
self.surface_roughness = np.asarray(surface_roughness, dtype = object)
# Profiles are kept as per-pipe lists so each pipe may use its own number of points.
self.undrained_shear_strength_depth_array = [
np.asarray(profile, dtype = float)
for profile in (undrained_shear_strength_depth_array or [])
]
self.undrained_shear_strength_value_array = [
np.asarray(profile, dtype = float)
for profile in (undrained_shear_strength_value_array or [])
]
self.submerged_unit_weight = np.asarray(submerged_unit_weight, dtype = float)
@staticmethod
def _depth_geometry(outer_diameter, maximum_depth):
'''
Calculate the depth, width, and penetrated area for a given pipe outer diameter.
'''
# Create an array of depth values from 0.001 m to the maximum depth of the
# undrained shear strength profile at 0.001 m spacing.
depth = np.arange(1, int(round(maximum_depth / 0.001)) + 1) * 0.001
# Calculate the width for a given pipe outer diameter (B).
width = np.where(
depth < outer_diameter / 2,
2 * np.sqrt(np.maximum(outer_diameter * depth - depth ** 2, 0)),
outer_diameter
)
# Calculate the penetrated area for a given pipe outer diameter (Abm)
penetrated_area = np.where(
depth < outer_diameter / 2,
np.arcsin(width / outer_diameter) * outer_diameter ** 2 / 4
- width * outer_diameter / 4 * np.cos(np.arcsin(width / outer_diameter)),
np.pi * outer_diameter ** 2 / 8
+ outer_diameter * (depth - outer_diameter / 2)
)
return depth, width, penetrated_area
def _model_inputs(self):
'''
Yield shared per-pipe inputs and depth geometry for the undrained models.
'''
for i, outer_diameter in enumerate(self.total_outer_diameter):
depth_profile = self.undrained_shear_strength_depth_array[i]
value_profile = self.undrained_shear_strength_value_array[i]
depth, width, penetrated_area = self._depth_geometry(
outer_diameter, np.max(depth_profile)
)
yield (
outer_diameter,
self.surface_roughness[i],
depth_profile,
value_profile,
self.submerged_unit_weight[i],
depth,
width,
penetrated_area
)
@staticmethod
def _shear_strength_gradient(depth, depth_profile, value_profile):
'''
Calculate the local undrained shear strength gradient along the depth array.
'''
return np.gradient(np.interp(depth, depth_profile, value_profile), depth)
@staticmethod
def _reference_shear_strength(
outer_diameter,
depth,
width,
depth_profile,
value_profile
):
'''
Calculate the reference depth and shear strength for Model 1.
'''
# Calculate the reference depth for Model 1 (zsu,0).
reference_depth = np.where(
depth < outer_diameter / 2 * (1 - np.sqrt(2) / 2),
0.0,
depth + outer_diameter / 2 * (np.sqrt(2) - 1) - width / 2
)
# Calculate the reference shear strength for Model 1 (su,0).
reference_shear_strength = np.interp(
reference_depth, depth_profile, value_profile
)
return reference_depth, reference_shear_strength
@staticmethod
def _interpolated_friction_factor(
width,
shear_strength_gradient,
seabed_shear_strength,
surface_roughness
):
'''
Calculate the interpolated pipe-soil friction factor.
'''
friction_factor = np.divide(
shear_strength_gradient * width,
seabed_shear_strength,
out=np.zeros_like(width),
where=seabed_shear_strength != 0
)
friction_factor_reference = np.linspace(0, 16, 17)
friction_factor_values = np.array(
[1.00, 1.12, 1.19, 1.24, 1.28, 1.31, 1.34, 1.36, 1.38,
1.39, 1.40, 1.41, 1.42, 1.43, 1.44, 1.445, 1.45]
if surface_roughness == 'Smooth' else
[1.00, 1.23, 1.36, 1.44, 1.50, 1.55, 1.59, 1.61, 1.64,
1.66, 1.67, 1.69, 1.70, 1.71, 1.72, 1.730, 1.74]
)
return np.interp(
friction_factor,
friction_factor_reference,
friction_factor_values
)
@staticmethod
def _depth_correction_factor(
seabed_shear_strength,
reference_shear_strength,
bearing_capacity,
reference_depth,
width,
bearing_capacity_factor=5.14
):
'''
Calculate the depth correction factor for Model 1.
'''
average_shear_strength_above = (
seabed_shear_strength + reference_shear_strength
) / 2
average_shear_strength_below = (
bearing_capacity / width / bearing_capacity_factor
)
return 0.3 * np.divide(
average_shear_strength_above,
average_shear_strength_below,
out=np.zeros_like(width),
where=average_shear_strength_below != 0
) * np.arctan2(reference_depth, width)
[docs]
def downward_undrained_model1(self): # pylint: disable=too-many-locals
"""
Calculate vertical penetration resistance using undrained Model 1.
The calculation follows the bearing-capacity formulation in the attached
PDF. For each pipe, the method evaluates depths at 0.001 m spacing from
0.001 m to the maximum depth of ``undrained_shear_strength_depth_array``
and prepends the zero-depth point to the returned arrays. The
reference shear strength is linearly interpolated at ``zsu,0`` from the
profile defined by ``undrained_shear_strength_depth_array`` and
``undrained_shear_strength_value_array``, and the buoyancy contribution
uses the supplied submerged unit weight.
Returns
-------
depth_arrays : np.ndarray
Depth arrays with shape ``(n_pipes, n_depths)``.
vertical_bearing_capacity_arrays : np.ndarray
Vertical penetration resistance arrays with shape
``(n_pipes, n_depths)``, in N/m.
Examples
--------
>>> psi = PSI(
... total_outer_diameter=[0.2731],
... surface_roughness=['Smooth'],
... undrained_shear_strength_depth_array=[[0.0, 0.3, 0.5, 1.0, 2.0]],
... undrained_shear_strength_value_array=[[1200, 5300, 3000, 4200, 6300]],
... submerged_unit_weight=[5500.0]
... )
>>> depths, capacities = psi.downward_undrained_model1()
>>> depths.shape, capacities.shape, float(depths[0, 0]), round(float(depths[0, -1]), 5), float(capacities[0, 0])
((1, 2001), (1, 2001), 0.0, 2.0, 0.0)
"""
depth_arrays = []
vertical_bearing_capacity_arrays = []
for (
outer_diameter,
surface_roughness,
shear_strength_depth_profile,
shear_strength_value_profile,
submerged_unit_weight,
depth,
width,
penetrated_area
) in self._model_inputs():
# Calculate reference depth and shear strength for Model 1.
reference_depth, reference_shear_strength = self._reference_shear_strength(
outer_diameter,
depth,
width,
shear_strength_depth_profile,
shear_strength_value_profile
)
# Local gradient and mudline strength derived from the profile.
shear_strength_gradient = self._shear_strength_gradient(
depth,
shear_strength_depth_profile,
shear_strength_value_profile
)
seabed_shear_strength = np.interp(
0.0, shear_strength_depth_profile, shear_strength_value_profile
)
# Calculate the interpolated pipe-soil friction factor based on surface roughness and shear strength.
friction_factor = self._interpolated_friction_factor(
width,
shear_strength_gradient,
seabed_shear_strength,
surface_roughness
)
# Define the bearing capacity factor for Model 1.
bearing_capacity_factor = 5.14
# Calculate the bearing capacity for Model 1.
bearing_capacity = friction_factor * (
bearing_capacity_factor *
reference_shear_strength + shear_strength_gradient * width / 4
) * width
# Calculate the depth correction factor for Model 1.
depth_correction_factor = self._depth_correction_factor(
seabed_shear_strength,
reference_shear_strength,
bearing_capacity,
reference_depth,
width,
bearing_capacity_factor
)
# Calculate the vertical bearing capacity for Model 1.
vertical_bearing_capacity = (
bearing_capacity * (1 + depth_correction_factor)
+ submerged_unit_weight * penetrated_area
)
# Append the calculated depth and vertical bearing capacity arrays to the results lists.
depth_arrays.append(np.append(0, depth))
vertical_bearing_capacity_arrays.append(
np.append(0, vertical_bearing_capacity)
)
return np.array(depth_arrays), np.array(vertical_bearing_capacity_arrays)
[docs]
def downward_undrained_model2(self): # pylint: disable=too-many-locals
"""
Calculate vertical penetration resistance using undrained Model 2.
Model 2 follows the alternative formulation in the attached PDF. It
combines the minimum of the two resistance-factor terms with the soil
buoyancy contribution. For each pipe, the method evaluates depths at
0.001 m spacing from 0.001 m to the maximum depth of
``undrained_shear_strength_depth_array`` and prepends the zero-depth
point to the returned arrays. The shear strength at the pipe invert is
linearly interpolated from the profile defined by
``undrained_shear_strength_depth_array`` and
``undrained_shear_strength_value_array``.
Returns
-------
depth_arrays : np.ndarray
Depth arrays with shape ``(n_pipes, n_depths)``.
vertical_bearing_capacity_arrays : np.ndarray
Vertical penetration resistance arrays with shape
``(n_pipes, n_depths)``, in N/m.
Examples
--------
>>> psi = PSI(
... total_outer_diameter=[0.2731],
... surface_roughness=['Smooth'],
... undrained_shear_strength_depth_array=[[0.0, 0.3, 0.5, 1.0, 2.0]],
... undrained_shear_strength_value_array=[[1200, 5300, 3000, 4200, 6300]],
... submerged_unit_weight=[5500.0]
... )
>>> depths, capacities = psi.downward_undrained_model2()
>>> depths.shape, capacities.shape, float(depths[0, 0]), round(float(depths[0, -1]), 5), float(capacities[0, 0])
((1, 2001), (1, 2001), 0.0, 2.0, 0.0)
"""
depth_arrays = []
vertical_bearing_capacity_arrays = []
for (
outer_diameter,
_,
shear_strength_depth_profile,
shear_strength_value_profile,
submerged_unit_weight,
depth,
_,
penetrated_area
) in self._model_inputs():
# Calculate the shear strength at the pipe invert for Model 2.
shear_strength = np.interp(
depth,
shear_strength_depth_profile,
shear_strength_value_profile
)
# Calculate the resistance factor for Model 2 based on depth and outer diameter.
resistance_factor = np.minimum(
6 * (depth / outer_diameter) ** 0.25,
3.4 * (10 * depth / outer_diameter) ** 0.5
)
# Calculate the vertical bearing capacity for Model 2.
vertical_bearing_capacity = (
resistance_factor
+ np.divide(
1.5 * submerged_unit_weight * penetrated_area,
outer_diameter * shear_strength,
out=np.zeros_like(depth),
where=shear_strength != 0
)
) * outer_diameter * shear_strength
# Append the calculated depth and vertical bearing capacity arrays to the results lists.
depth_arrays.append(np.append(0, depth))
vertical_bearing_capacity_arrays.append(
np.append(0, vertical_bearing_capacity)
)
return np.array(depth_arrays), np.array(vertical_bearing_capacity_arrays)
