Before diving into the problem, I want to acknowledge Marco Trezzi, you can connect with him on on LinkedIn and Medium. Marco is a brilliant geotechnical engineer and is keen on Python coding. Much of the article’s merit is due to him since he first discovered the method, proposed to write something about and developed a big part of the code. Sadly it is not possible to co-author on Medium, so this content is being published on my account.

Usually, in common engineering practice, there is the habit of designing the length of piles in group taking into consideration the most loaded one. In most cases, there is a back and forth between the geotechnical and structural engineers in which the first give to the second piles springs and bearing capacity values, and the second extract the highest load as the dimensioning one. Recently Marco stepped into this article from Di Laora, De Sanctis and Aversa from “Acta Geotecnica”. We found the contribution valuable, and we thought it could be interesting for professionals to develop a tool to implement this formulation. Reference shall be made to the article for a detailed explanation of the theoretical background involved.

The underlying principle assumes the axial capacity increased by the contribution of the piles in traction on the opposite side of load directions, as in the following figure.

The final objective is to create an interaction diagram that provides a distance from failure for each load combination.

Each node of the domain represents a different load configuration where, according to the direction of the moment and the axial action acting of the raft, the piles are loaded differently. Just to give an example: if one look at the right side of the diagram, when no moment is acting on the foundation, of course the ultimate bearing capacity of the group is n times the bearing capacity of a single pile (n is the number of piles). The same in traction (left side of the diagram). To compute the bearing capacity for each “node” of the diagram, the following formula is used:

The first formula refers to the plot’s upper side (moment major than 0 or clockwise rotation) and the second to the left side of the diagram (moment minor than 0 or anti-clockwise rotation). ξ value is defined as follows:

α is the angle between the moment in x an y directions acting on the raft foundation

After some introduction, let’s dive into the code and start importing the packages we are going to use for our purpose.

import pandas as pd
import numpy as np
from matplotlib.lines import Line2D
import matplotlib.pyplot as plt

First of all, we create 2 CSV containing respectively the piles’ positions with the bearing capacity (compression and tension) and one with the global loads acting on the foundation. Then, the CSVs are imported as data frames.

df_piles = pd.read_csv('piles_group.csv', delimiter=',')
display(df_piles)
df_loads = pd.read_csv('piles_loads.csv', delimiter=',')
display(df_loads)

Let’s calculate the resulting moment and the angle α for each load combination.

def calc_alfa_M(df): #REQUIRES DF LOADS
    mtot = []
    alpha = []
    for mx,my in zip(df['Medx [kNm]'],df['Medy [kNm]']):
        m = np.sqrt(mx**2+my**2)
        
        if mx>0:
            a = np.degrees(np.arccos(my/m))
            mtot.append(m)
        else:
            a = 360 - np.degrees(np.arccos(my/m))
            mtot.append(-m)
        alpha.append(a)
        
    df['alpha [°]'] = alpha
    df['Medtot [kNm]'] = mtot
    
    return df

df_loads = calc_alfa_M(df_loads)
display(df_loads)

The strategy is to build as much domains as the values of α available in our load combinations. The domain will then be built by fixing α and plotting the loads values. First of all, we define how to compute ξ value.

def calc_xi (a,df): #REQUIRES DF PILES
    xi = []
    for x,y in zip(df['X [m]'],df['Y [m]']):
        xi_j = x * np.cos(np.radians(a)) - y * np.sin(np.radians(a))
        xi.append(xi_j)
    df['xi [-]'] = xi
    return df

Now we are going to apply the previous formula to the domain calculation. Note the p and k indexes refer respectively to the upper and lower side of the domain.

def calc_domain(a,df): #REQUIRES DF PILES

    Nu = []
    Mu = []
    
    df['xi_c'] = [i*j for i,j in zip(df['Rcd [kN]'],df['xi [-]'])]
    df['xi_t'] = [i*j for i,j in zip(df['Rtd [kN]'],df['xi [-]'])]
    
    df_sort = df.sort_values(['xi [-]'], ascending = True)
    
    rcd_somm = df_sort['Rcd [kN]']
    rtd_somm = df_sort['Rtd [kN]']
    rcd_psi_somm = df_sort['xi_c']
    rtd_psi_somm = df_sort['xi_t']
    
    #UPPER SIDE DOMAIN
    for i in range(len(rcd_somm)):
    
        Nu.append(sum(rcd_somm[: i]) - sum(rtd_somm[i :]))
        Mu.append(sum(- rcd_psi_somm[: i]) + sum(rtd_psi_somm[i:]))
    
    #LOWER SIDE DOMAIN
    for i in range(len(rcd_somm)):    
 
        Nu.append(sum(rcd_somm[i:]) - sum(rtd_somm[:i]))
        Mu.append(sum(rcd_psi_somm[: i]) - sum(rtd_psi_somm[i:]))
    
    #LAST POINT TO CLOSE THE DOMAIN
    Nu.append(Nu[0])
    Mu.append(Mu[0])
    Nu = np.array(Nu)
    Mu = np.array(Mu)
    
    return Mu,Nu

Now we iterate over the different α values and plot the domain:

# ALL AVAILABLE MARKERS OF MATPLOTLIB
markers = [k for k in Line2D.markers.keys()]

# COLOR FUNCTION 
def get_cmap(n, name='hsv'):
    return plt.cm.get_cmap(name, n)

alpha_calc = df_loads['alpha [°]'].unique()
cmap = get_cmap(len(alpha_calc) + 1)

fig, ax = plt.subplots(figsize=(16,8))

# ITERATION OVER DIFFERENT VALUES OF ALPHA
for n,alpha in enumerate(alpha_calc):
    
    #CALULATE XI
    df_piles = calc_xi(alpha,df_piles)
    
    #CALCULATE DOMAIN
    Mu,Nu = calc_domain(alpha,df_piles)
    
    #PLOT
    ax.plot(Nu, Mu, label='$α$ = ' + str(int(alpha)) + ' °', marker=markers[n], color=cmap(n), markersize=10)
    loads_plot = df_loads[df_loads['alpha [°]'] == alpha]
    ax.scatter(loads_plot['Ned [kN]'], loads_plot['Medtot [kNm]'], color=cmap(n), marker=markers[n], s=100)

ax.legend()
plt.grid()
fig.tight_layout()
ax.set_xlabel('N [kN]')
ax.set_ylabel('M [kNm]')
plt.show()

For each load combination, we can then compute the utilization ratio. The ratio is meant as the distance of the point from the domain considering constant eccentricity, the slope of the segment connecting the origin to the load combination. For this scope, we are going to use the following package.

import shapely
from shapely.geometry import LineString, Point

A function is defined to iterate over all the segments componing the domain to find the minimum distance of the load combination from the
boundary.

def calc_FS(M,N,Mu,Nu):
    
    # INITIALIZE THE FIRST POINT OF THE SEGMENT
    x0,y0 = Nu[0],Mu[0]
   
    # FORCE LINE (COSTANT ECCENTRICITY)
    A = (min(Nu), min(Nu)*M/N)
    B = (max(Nu), max(Nu)*M/N)
    
    # DUMMY DISTANCE
    dist = (max(Nu)-min(Nu))*10
    
    #DISTANCE OF THE LOAD COMBINATION FROM THE ORIGIN
    dist_force = np.sqrt((M)**2+(N)**2)
    
    #ITERATE OVER ALL THE SEGMENTS COMPONING THE DOMAIN
    for i in range(1,len(Mu)):
             
        # SEGMENT LINE 
        C = (x0, y0)
        D = (Nu[i], Mu[i])
        
        line1 = LineString([A, B])
        line2 = LineString([C, D])
        
        # POINT OF INTERSECTION
        int_pt = line1.intersection(line2)
        
        #CHECK IF THE TWO LINES INTERSECT
        if str(int_pt) != 'LINESTRING EMPTY':

            dist_t = np.sqrt((int_pt.x-N)**2+(int_pt.y-M)**2)
            
            # UPDATE THE MINIMUM DISTANCE
            if dist_t<dist:
    
                dist = dist_t
                x = int_pt.x
                y = int_pt.y
        
                dist_tot = np.sqrt((x)**2+(y)**2)
                UR = dist_force/dist_tot
                
        x0,y0 = Nu[i],Mu[i]
                
    return UR,x,y,dist

Iterating, we apply the function to all the load combination in our data frame.

ur = []

for i in df_loads.index:
    Mu,Nu = calc_domain(df_loads.loc[i,'alpha [°]'],df_piles)
    UR,x,y,dist = calc_FS(df_loads.loc[i,'Medtot [kNm]'],df_loads.loc[i,'Ned [kN]'],Mu,Nu)
    ur.append(float(UR))

df_loads['UR [-]'] = ur
display(round(df_loads,1))

Concluding, we developed a promising method using Python to compute the bearing capacity of a group of piles to generate substantial savings in foundation design.

All the code and data frames are available in the following GitHub repository.