Uno script Python per la costruzione delle curve granulometriche

Come creare delle curve granulometriche con le librerie Matplotlib e NumPy

L’uso di Python nelle geoscienze è molto vantaggioso grazie al grande numero di librerie che permettono di eseguire operazioni molto elaborate.
Ad esempio in questo articolo mostro prima come come manipolare le giaciture geologiche e poi come fare analisi strutturali con gli stereonet attraverso il modulo mplstereonet. In questo articolo vedremo invece come costruire delle curve granulometriche con Python.

La mia tesi magistrale ha previsto una fase di lavoro in campagna nella quale ho anche prelevate campioni di terreno da analizzare presso il laboratorio di geotecnica del mio dipartimento. Qui ho stimato la granulometria dei campioni sia con il metodo della setacciatura che con quello della sedimentazione. Per lo sviluppo di questo script ho usato i dati granulometrici di uno di questi campioni.

setacciatrice

 

Un pò di teoria

Con il metodo per setacciatura, la curva granulometrica si costruisce conoscendo le percentuali di passante ad ogni setaccio. Il passante ad ogni setaccio si calcola con la seguente formula:

(1)   \begin{equation*} P_i=\frac{P_s-\sum P_t}{P_s}*100 \end{equation*}

setaccio

Con il metodo della della sedimentazione la curva granulometrica, nel quale si usa un densimetro per prendere delle misure di densità per gradini di tempo, si stima il diametro delle particelle con la legge di Stokes.

(2)   \begin{equation*} D=\sqrt{\frac{1800*\eta_L}{\gamma_s-\gamma_l}}*\sqrt{\frac{H_r}{\delta_t*60}} \end{equation*}

La percentuale di particelle (P%) passanti al setaccio ideale, di diametro uguale a quello del diametro calcolato D, è data dalla seguente espressione:

(3)   \begin{equation*} P(\%)= R''*\frac{100}{P_s}*\frac{\gamma_s}{\gamma_s-\gamma_l} \end{equation*}

sedimentazione

 

Per approfondire questi elementi, anche dal punto di vista operativo, vi rimando alla normativa che regola questo tipo di prove [AGI(1994,1996) e ASTM(2006) – D422-63] e a questo manuale.

Il codice Python

Innanzitutto vanno importare la libreria NumPy e matplotlib:

import numpy as np
import matplotlib.pyplot as plt

Dopo è stata implementata la parte di codice riguardante la porzione di curva ottenuta dalla setacciatura:

peso_camp =[777.04]
diametri=[12.5, 6.3, 4.75, 2, 1, 0.425, 0.250, 0.150, 0.106, 0.075, 0.045]
peso_grani=[16.01, 82.55, 55.08, 52.61, 101,98, 55.61, 29.05, 14.33, 11.05, 11.33, 28.71]
parziale=[]
trattenuto=[]
passante=[]


for i in peso_grani:
    print(parziale.append((i*100)/777.04))

val=0
val2=parziale[0]

for i in parziale:
    if i == parziale[0]:
        print(trattenuto.append(i))
    elif i == parziale[1]:
        print(trattenuto.append(val2+i))
    elif i == parziale[2]:
        print(trattenuto.append(val2+parziale[1]+i))
    elif i == parziale[3]:
        print(trattenuto.append(val2+parziale[1]+parziale[2]+i))
    elif i == parziale[4]:
        print(trattenuto.append(val2+parziale[1]+parziale[2]+parziale[3]+i))
    elif i == parziale[5]:
        print(trattenuto.append(val2+parziale[1]+parziale[2]+parziale[3]+parziale[4]+i))
    elif i == parziale[6]:
        print(trattenuto.append(val2+parziale[1]+parziale[2]+parziale[3]+parziale[4]+parziale[5]+i))
    elif i == parziale[7]:
        print(trattenuto.append(val2+parziale[1]+parziale[2]+parziale[3]+parziale[4]+parziale[5]+parziale[6]+i))
    elif i == parziale[8]:
        print(trattenuto.append(val2+parziale[1]+parziale[2]+parziale[3]+parziale[4]+parziale[5]+parziale[6]+parziale[7]+i))
    elif i == parziale[9]:
        print(trattenuto.append(val2+parziale[1]+parziale[2]+parziale[3]+parziale[4]+parziale[5]+parziale[6]+parziale[7]+parziale[8]+i))
    elif i == parziale[10]:
        print(trattenuto.append(val2+parziale[1]+parziale[2]+parziale[3]+parziale[4]+parziale[5]+parziale[6]+parziale[7]+parziale[8]+parziale[9]+i))
    

a = 100
for i in trattenuto:
    if i == trattenuto[0]:
        print(passante.append(a-i))
    elif i == trattenuto[1]:
        print(passante.append(a-i))
    elif i==trattenuto[2]:
        print(passante.append(a-i))
    elif i==trattenuto[3]:
        print(passante.append(a-i))
    elif i==trattenuto[4]:
        print(passante.append(a-i))
    elif i==trattenuto[5]:
        print(passante.append(a-i))
    elif i==trattenuto[6]:
        print(passante.append(a-i))
    elif i==trattenuto[7]:
        print(passante.append(a-i))
    elif i==trattenuto[8]:
        print(passante.append(a-i))
    elif i==trattenuto[9]:
        print(passante.append(a-i))
    elif i==trattenuto[10]:
        print(passante.append(a-i))
    else:
        print(False)

Invece la porzione di curva data dalla sedimentazione si ottiene con il seguente codice:

Tempi = [1,2,4,8,15,30,60,120,240,480,1440,2880]
Lettura = [25,24,23,22,20.5,18.5,17,15,13,12,10.5,8]
Temperatura = [17.8,17.7,17.8,17.9,17.9,18.1,18.2,18.3,18.4,18.5,17.9,18.3]
C_t = [-0.4, -0.4, -0.4, -0.35, -0.35, -0.3, -0.3, -0.3, -0.25, -0.25, -0.35, -0.3]

gamma_s= 2.7
gamma_l= 1.0
Ps=40
V=1000
P_set=44.03  #passante al setaccio 0.075

nL = []


for i in Temperatura:
    a=10**(-5)
    b=i**2
    print(nL.append((1.81*a)/(1+0.034*i+0.00022*b)))

Lettura_menisco=[j+0.5 for j in Lettura]


Hr=[16.43-(0.2747*i) for i in Lettura_menisco]

tempo_s=[k*60 for k in Tempi]

z=[(1800*i)/(gamma_s-1) for i in nL]

if len(Lettura_menisco) == len(C_t):
    Lettura_corr=[w-5+h for (w,h) in zip(Lettura_menisco,C_t)]
else:
    print("Errore: le lunghezze di Lettura_menisco e C_t sono diverse")
    print("len(Lettura_menisco) = ", len(lettura_menisco))
    print("len(C_t) = ", len(C_t))
    exit(1)


if len(Hr) == len(tempo_s):
    minnie=[a/b for (a,b) in zip(Hr,tempo_s)] #rapporto tra gli elementi di Hr e tempo_s
    print(minnie)
else:
    print("Errore: le lunghezze di Hr e tempo_s sono diverse")
    print("len(tempo_s) = ", len(tempo_s))
    print("len(Hr) = ", len(Hr))
    exit(1)
    
#DIAMETRO SEDIMENTAZIONE 

if len(z) == len(minnie):
    diametro=[(k*i)**0.5 for (k,i) in zip(z,minnie)]
    print(diametro)
else:
    print("Errore: le lunghezze di z e minnie sono diverse")
    print("len(z) = ", len(z))
    print("len(minnie) = ", len(minnie))
    exit(1)

La curva granulometrica totale viene costruita con il seguente blocco di codice:

a=[12.5, 6.3, 4.75, 2, 1, 0.425, 0.25, 0.15, 0.106, 0.075, 0.045, 0.04239596508749394, 0.030450234096941694, 0.02180705659519285, 0.015611315129263687, 0.011627809103480602, 0.008410321957201528, 0.006047781524118338, 0.004370905906703727, 0.003155765761466145, 0.0022526383784678143]
b=[97.93961700813342, 87.31596828992073, 80.22753011427983, 73.45696489241223, 60.45892103366622, 47.84695768557603, 40.69031195305261, 36.95176567486873, 35.107587768969424, 33.68552455472047, 32.22742715947699, 31.62176470588235, 29.962058823529404, 27.341470588235293, 23.93470588235294, 21.31411764705882, 17.819999999999997, 14.413235294117648, 12.666176470588235, 9.870882352941177, 5.590588235294117]

#creo una figura di 8x6 pollici con risoluzione di 80 punti/pollice
plt.figure(figsize=(15,10))

# creo un unico sistema di assi cartesiani nella figura
plt.subplot(1, 1, 1) #(1riga,1colonna,grafico1)

# Plottaggio grafico
grafico=plt.plot(a, b, color="blue", linewidth=1.5, linestyle="-", label='CURVA GRANULOMETRICA')

# Impostazioni asse x

plt.xscale('log')

# Impostazioni asse y
plt.ylim(0, 100)
plt.yticks([0,10,20,30,40,50,60,70,80,90,100])

#griglia
plt.grid(True)

#label assi
plt.xlabel('DIAMETRI')
plt.ylabel('PASSANTE (%)')

#legenda
plt.legend()

#salvo la figura
plt.savefig('curva.png')

Il risultato è il seguente:

curva granulometrica

Per avere chiarimenti, darmi qualche consiglio su implementazioni future o per avere il codice sorgente potete lasciare un commento nel form sottostante.

Ciao, mi chiamo Antonio Nirta e sono un geologo. Mi sono laureato all’Università di Pisa e dal 2017 svolgo la libera professione.

Attraverso i post presenti in questo blog cercherò di trasmettere la mia passione per le Scienze Geologiche e di fornire informazioni che riguardano il bellissimo settore della Geologia.

2 thoughts on “Uno script Python per la costruzione delle curve granulometriche

  1. Complimenti per lo script. Dovrebbe esserci però un errore: nella sezione riguardante la setacciatura, nella lista peso_grani, il valore dovrebbe essere 101.98 e non 101,98 (che sono due valori in quest’ultimo caso)

    1. Ciao Simone, grazie per i complimenti. Controllo quella porzione di codice per vedere se da qualche problema. Inoltre sto pensando ad implementare il codice, magari inserendo alla fine la classificazione granulometrica.
      A presto, Antonio.

Lascia un commento

Il tuo indirizzo email non sarà pubblicato. I campi obbligatori sono contrassegnati *