Come costruire un curva granulometrica con Python

Usare 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.

Che cosa è una curva granulometrica

La curva granulometrica è la rappresentazione grafica dei risultati ottenuti con un’analisi granulometrica. Una analisi granulometrica raggruppa, in diverse classi di grandezza, le particelle costituenti il terreno, e di determinare successivamente le percentuali in peso di ciascuna classe, riferendolo al peso secco del campione iniziale.

Per approfondire come si fa un’analisi granulometrica ti rimando a questo post.

Il codice Python

Innanzitutto vanno importare la libreria NumPy e matplotlib:

[python]
import numpy as np
import matplotlib.pyplot as plt
[/python]

La porzione di codice per costruire la parte di curva che si ricava dalla setacciatura è la seguente:

[python]
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:
    trattenuto = list(np.cumsum(parziale))[:-1]

print(trattenuto)

passante = [100-i for i in trattenuto]

print(passante)   
[/python]

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

[python]
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)
[/python]

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

[python]
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')
[/python]

Il risultato è il seguente:

Approfondimenti

Per approfondire l’argomento trattato in questo post vi consiglio i seguenti testi: