Update suction capacity model and add profile module for anchors in homogeneous soil profiles

Author: Moreno

famodel/anchors/anchors_famodel/capacity_suction.py

@@ -40,7 +40,7 @@ def getCapacitySuction(D, L, zlug, H, V, soil_type, gamma, Su0=None, k=None, phi
         Maximum vertical capacity [kN]
     '''  
 
-    lambdap = L/D; m = -2/3;         # Suction pile slenderness ratio
+    lambdap = L/D; m = 2/3;         # Suction pile slenderness ratio
     t = (6.35 + D*20)/1e3            # Suction pile wall thickness (m), API RP2A-WSD
     rlug = D/2                       # Radial position of the lug
     thetalug = 5                     # Angle of tilt misaligment, default is 5. (deg)
@@ -83,7 +83,7 @@ def calc_delta(Dr_val):
     if soil_type == 'clay':
         # Definitions for cohesive soils
         Nc = min (6.2*(1 + 0.34*np.arctan(lambdap)),9)   # End-bearing capacity factor
-        ez = (Su0*L**2/2 + k*L**3/3)/(Su0*L + k*L**2/2)
+        ez = (Su0*L**2/2 + k*L**3/3)/(Su0*L + k*L**2/2); print(ez)
         Np_fixed = 10.25; Np_free = 4                    # From Np vs L/D chart from CAISSON_VHM
         Su_av_L = Su0 + k*zlug                           # Undrained shear strength values (average) 
         Su_tip = Su0 + k*L                               # Undrained shear strength values (tip)
@@ -104,7 +104,7 @@ def calc_delta(Dr_val):
         M = - V*rlugTilt(rlug,zlug,thetalug) - H*(zlugTilt(rlug,zlug,thetalug) - ez)
         def f(Hmax):
              return m*(Hmax/(L*D*(Su0 + k*zlug)) - Np_fixed) + M*(Hmax/(L*D*(Su0 + k*zlug))/(Hmax*L))
-        Hmax = fsolve(f,5);
+        Hmax = fsolve(f,5)
 
         # Torsion capacity
         Fo = PileSurface(L, D)*alpha*Su_av_L
@@ -126,6 +126,15 @@ def f(Hmax):
         # "Leaking"        
         Vmax3 = (PileWeight(L, D, t, rhows) + PileSurface(L, D)*alphastar*Su_av_L + SoilWeight(L, D, t, gamma))                  
         Vmax = min(Vmax1, Vmax2, Vmax3)
+        
+        print("\n--- Parameter-Based Version ---")
+        print(f"Su_av_L        = {Su_av_L:.3f} kPa")
+        print(f"sigma'_v(zlug) = {sigma_v_eff:.3f} kPa")
+        print(f"psi_val        = {psi_val:.3f}")
+        print(f"alpha (API)    = {alpha:.3f}")
+        print(f"Hmax           = {Hmax[0]:.2f} kN")
+        print(f"Vmax           = {Vmax:.2f} kN")
+
 
     elif soil_type == 'sand':
         # Definition for non-cohesive soils
@@ -159,15 +168,22 @@ def f(Hmax):
             # return np.e**(-depth) - 1 + depth
         # Ze = D/(4*7); Zi = D/(4*5)    
         # Vmax = 7*gamma*Ze**2*y(L/Ze)*PileSurface(L, D)/L + 5*gamma*Zi**2*y(L/Zi)*PileSurface(L,(D - 2*t))/L
-     
+        
+        print("\n--- Parameter-Based (Sand) ---")
+        print(f"phi        = {phi:.2f} deg")
+        print(f"gamma      = {gamma:.2f} kN/m3")
+        print(f"deltastar  = {deltastar:.2f} -")
+        print(f"sigma_av_L = {sigma_av_L:.2f} kN")
+        print(f"sigma_tip  = {sigma_tip:.2f} kN")
+   
     # Pile weight (inc. stiffening plus vent) assessed as a factor
     Wp = 1.10*PileWeight(L, D, t, (rhows + rhow)) 
     # Submerged weight of the soil plug
     Ws = SoilWeight(L, D, t, gamma)
 
     # Capacity envelope
     aVH = 0.5 + lambdap; bVH = 4.5 + lambdap/3 
-    # print('Env. exp = ' +str(aVH)+'   '+str(bVH))
+    print('Env. exp = ' +str(aVH)+'   '+str(bVH))
     UC = (H/Hmax)**aVH + (V/Vmax)**bVH      
     x = np.cos(np.linspace (0, np.pi/2, 100))
     y = (1 - x**bVH)**(1/aVH)
@@ -194,7 +210,7 @@ def f(Hmax):
         resultsSuction['UC'] = UC[0]                # Unity check in clay
     elif soil_type == 'sand':
         resultsSuction['UC'] = UC                   # Unity check in sand
-    resultsSuction['Weight'] = Wp              # Dry weight of the suction pile (kN)
+    resultsSuction['Weight'] = Wp                   # Dry weight of the suction pile (kN)
     resultsSuction['Weight Soil'] = Ws              # Submerged weight of the soil plug (kN)
     resultsSuction['t'] = t                         # Pile thikness in [m]
 
@@ -291,8 +307,27 @@ def SoilWeight(Len, Dia, tw, gamma_soil):
     resultsSuctionSimp['Horizontal max.'] = Hmax    # Capacity at specified loading angle
     resultsSuctionSimp['Vertical max.'] = Vmax      # Capacity at specified loading angle
     resultsSuctionSimp['UC'] = UC                   # Unity check
-    resultsSuctionSimp['Weight'] = Wp          # Pile weight in kN
+    resultsSuctionSimp['Weight'] = Wp               # Pile weight in kN
     resultsSuctionSimp['Weight Soil'] = Ws          # in kN
     resultsSuctionSimp['t'] = t
 
-    return resultsSuctionSimp    
+    return resultsSuctionSimp  
+
+if __name__ == '__main__':
+    
+    D = 2.0          # Diameter in meters
+    L = 10.0         # Length in meters
+    zlug = (2/3)*L   # Padeye depth
+    H = 1500.0       # Horizontal load in kN
+    V = 1000.0       # Vertical load in kN
+    
+    gamma = 8
+    Su0 = 0
+    k = 10
+    
+    phi = 50
+    Dr = 75
+    
+    results_clay = getCapacitySuction(D, L, zlug, H, V, 'clay', gamma, Su0=Su0, k=k, phi=phi, Dr=Dr, plot=True) 
+    
+    results_sand = getCapacitySuction(D, L, zlug, H, V, 'sand', gamma, Su0=Su0, k=k, phi=phi, Dr=Dr, plot=True)   

famodel/anchors/anchors_famodel_profile/capacity_dandg.py

@@ -0,0 +1,245 @@
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.interpolate import interp1d
+from scipy import linalg
+from capacity_soils import rock_profile
+from capacity_pycurves import py_Lovera
+from capacity_plots import plot_pile
+
+def getCapacityDandG(profile, soil_type, L, D, zlug, H, V, plot=True):
+    '''Models a laterally loaded pile using the p-y method. The solution for
+    lateral displacements is obtained by solving the 4th order ODE, 
+    EI*d4y/dz4 - V*d2y/dz2 + ky = 0 using the finite difference method.
+    EI*d4y/dz4 - V*d2y/dz2 + K*z*dy/dz + ky = 0 using the finite difference method.
+
+    Assumes that EI remains constant with respect to curvature i.e. pile
+    material remains in the elastic region.
+
+    Parameters
+    ----------
+    profile : array
+        Rock profile as a 2D array: (z (m), UCS (MPa), Em (MPa))
+    soil_type : string
+        Select soil condition, 'rock'
+    L : float 
+        Pile length (m)
+    D : float 
+        Pile diameter (m)
+    zlug : float
+        Load eccentricity above the mudline or depth to mudline relative to the pile head (m)
+    H : float       
+        Horizontal load at pile lug elevation (N)
+    V : float          
+        Vertical load at pile lug elevation (N)
+    plot : bool
+        Plot the p-y curve and the deflection pile condition if True
+
+    Returns
+    -------
+    y : array
+        Lateral displacement at each node (n+1 real + 4 imaginary)
+    z : array
+        Node location along pile (m)
+    resultsDandG : dict
+        Dictionary with lateral, rotational, vertical and pile weight results
+    '''
+
+    n = 50; iterations = 10; loc = 2
+    nhuc = 1; nhu = 0.3          # Resistance factor (-)
+    delta_r = 0.08               # Mean roughness height (m)
+    
+    t = (6.35 + D*20)/1e3        # Pile wall thickness (m), API RP2A-WSD
+    E = 200e9                    # Elastic modulus of pile material (Pa)
+    rhows = 66.90e3              # Submerged steel specific weight (N/m3)
+    rhow = 10e3                  # Water specific weight (N/m3) 
+    
+    # Pile geometry
+    I = (np.pi/64.0)*(D**4 - (D - 2*t)**4)
+    EI = E*I
+    h = L/n                      # Element size
+    N = (n + 1) + 4              # (n+1) Real + 4 Imaginary nodes
+    
+    # Dry and wet mass of the pile    
+    def PileWeight(Len, Dia, tw, rho):
+        Wp = ((np.pi/4)*(Dia**2 - (Dia - 2*tw)**2)*Len)*rho
+        return Wp 
+
+    # Array for displacements at nodes, including imaginary nodes.
+    y = np.ones(N)*(0.01*D)      # An initial value of 0.01D was arbitrarily chosen
+
+    # Initialize and assemble array/list of p-y curves at each real node
+    z = np.zeros(N)
+    k_secant = np.zeros(N)
+    py_funs = []
+    DQ = []
+    z0, f_UCS, f_Em = rock_profile(profile)
+
+    for i in [0, 1]:              # Top two imaginary nodes
+        z[i] = (i - 2)*h
+        py_funs.append(0)
+        k_secant[i] = 0.0
+   
+    for i in range(2, n+3):       # Real nodes
+        z[i] = (i - 2)*h       
+        if z[i] < z0:
+             # No p-y curve above mudline
+             py_funs.append(lambda y_val: np.zeros_like(y_val))
+             k_secant[i] = 0.0
+             DQ.append(0.0)
+        else:
+             py_funs.append(py_Lovera(z[i], D, f_UCS, f_Em, zlug, z0, plot=True))            
+             UCS = f_UCS(z[i])
+             Em = f_Em(z[i])
+             SCR = nhuc*Em/(UCS*(1 + nhu))*delta_r/D
+             alpha = 0.36*SCR - 0.0005
+             fs = alpha*UCS
+             Dq = np.pi*D*fs*z[i]
+             DQ.append(Dq)
+             Vmax = PileWeight(L, D, t, rhows) + DQ[-1]
+             
+             k_secant[i] = py_funs[i](y[i])/y[i]
+
+    for i in [n+3, n+4]:         # Bottom two imaginary nodes
+        z[i] = (i - 2)*h
+        py_funs.append(0)
+        k_secant[i] = 0.0
+    
+    for j in range(iterations):
+        # if j == 0: print 'FD Solver started!'
+        y = fd_solver(n, N, h, EI, V, H, zlug, z0, k_secant)
+        if plot:
+            plt.plot(y[loc], k_secant[loc]*y[loc])
+
+        for i in range(2, n+3):
+            k_secant[i] = py_funs[i](y[i])/y[i]
+
+    if plot:
+        y1 = np.linspace(-2.*D, 2.*D, 500)
+        plt.plot(y1, py_funs[loc](y1))
+        plt.xlabel('y (m)'), plt.ylabel('p (N/m)')
+        plt.grid(True)
+        
+        fig, ax = plt.subplots(figsize=(3, 5))
+        y0 = np.zeros_like(z[2:-2])
+        # Plot original vertical pile
+        ax.plot(y0, z[2:-2], 'k', label='Original pile axis')
+        # Plot deflected shape
+        ax.plot(y[2:-2], z[2:-2], 'r', label='Deflected shape')
+        # Padeye marker
+        ax.plot(0, zlug, 'ko', label=f'Padeye (zlug = {zlug:.2f} m)')
+        # Mudline marker 
+        ax.axhline(z0, color='blue', linestyle='--', label=f'Mudline (z0 = {z0:.2f} m)')
+        
+        ax.set_xlabel('Lateral displacement (m)')
+        ax.set_ylabel('Depth (m)')
+        ax.set_xlim([-0.1*D, 0.1*D])
+        ax.set_ylim([L + 5, -2])  # Downward is positive z 
+        ax.grid(ls='--')
+        ax.legend()  
+            
+    resultsDandG = {
+        'Lateral displacement' : y[2],
+        'Rotational displacement' : np.rad2deg((y[2] - y[3])/h),
+        'Vertical max.' : Vmax,
+        'Weight pile' : PileWeight(L, D, t, (rhows + rhow))
+    }
+    
+    return y[2:-2], z[2:-2], resultsDandG
+
+def fd_solver(n, N, h, EI, H, V, zlug, z0, k_secant):
+    '''Solves the finite difference equations from 'py_analysis_1'. This function should be run iteratively for
+    non-linear p-y curves by updating 'k_secant' using 'y'. A single iteration is sufficient if the p-y curves
+    are linear.
+
+    Parameters
+    ----------
+    n : int
+        Number of elements (-)
+    N : int
+        Total number of nodes (-)
+    h : float
+        Element size (m)
+    EI : float
+        Flexural rigidity of the pile (Nm²)
+    H : float
+        Horizontal load at padeye (N)
+    V : float
+        Vertical load at padeye (N)
+    zlug : float
+        Padeye depth from pile head (m)
+    z0 : float
+        Mudline elevation from pile head (m)
+    k_secant : array
+        Secant stiffness at each node
+
+    Returns
+    -------
+    y : array
+        Lateral displacement at each node (n+1 real + 4 imaginary)
+    '''
+    
+    # Initialize and assemble matrix
+    X = np.zeros((N, N))
+
+    # (n+1) finite difference equations for (n+1) real nodes
+    for i in range(0,n+1):
+        X[i, i]   =  1.0
+        X[i, i+1] = -4.0 + V*h**2/EI
+        X[i, i+2] =  6.0 - 2*V*h**2/EI + k_secant[i+2]*h**4/EI
+        X[i, i+3] = -4.0 + V*h**2/EI
+        X[i, i+4] =  1.0
+
+    # Curvature at pile head
+    X[n+1, 1] =  1.0
+    X[n+1, 2] = -2.0
+    X[n+1, 3] =  1.0
+
+    # Shear at pile head
+    X[n+2, 0] = -1.0
+    X[n+2, 1] =  2.0 - V*h**2/EI
+    X[n+2, 2] =  0.0
+    X[n+2, 3] = -2.0 + V*h**2/EI
+    X[n+2, 4] =  1.0
+
+    # Curvature at pile tip
+    X[n+3, -2] =  1.0
+    X[n+3, -3] = -2.0
+    X[n+3, -4] =  1.0
+
+    # Shear at pile tip
+    X[n+4, -1] =  1.0
+    X[n+4, -2] = -2.0 + V*h**2/EI
+    X[n+4, -3] =  0.0
+    X[n+4, -4] =  2.0 - V*h**2/EI
+    X[n+4, -5] = -1.0
+
+    # Initialize vector q
+    q = np.zeros(N)
+    
+    # Always apply shear
+    # Index of the node where the horizontal load is applied (padeye)
+    zlug_index = int(zlug/h)
+    q[zlug_index] = 2*H*h**3
+
+    # Solve for displacement
+    y = linalg.solve(EI*X, q)
+
+    return y
+   
+if __name__ == '__main__':
+    
+    profile_rock = np.array([
+        [ 2.0, 1, 1e3],
+        [ 5.0, 2, 2e4],
+        [ 9.0, 4, 2e4],
+        [30.0, 6, 5e4]
+    ])
+    
+    D = 3.0           # Diameter (m)
+    L = 10.0          # Length (m)
+    zlug = 1          # Padeye elevation (m)
+    
+    y, z, results = getCapacityDandG(profile_rock, 'rock', L, D, zlug, H=9.5e12, V=3.0e6, plot=True)
+    
+    plot_pile(profile_rock, 'rock', y, z, D, L, profile_rock[0, 0], zlug)