Fattore di sicurezza con il “Metodo Monte Carlo”

Quando si ha a che fare con un sistema fisico si deve tener presente che esso può “instabilizzarsi”. Si può affermare che uno dei parametri che quantizzi al meglio e “a colpo d’occhio” quanto tale sistema sia distante dell’instabilità è il Fattore di Sicurezza. Esso può essere definito, ad esempio, come rapporto tra “forze resistenti” e “forze instabilizzanti”:

$F_S = \frac{F_R}{F_I}$

Se $F_S < 1$ vuol dire che il denominatore è maggiore del numeratore, ovvero che $F_R < F_I$ . In altri termini, le forze instabilizzanti sono superiori a quelle che si oppongono a tale instabilizzazione. Pertanto, a seconda dell’incertezza dei parametri, è bene mantenere tale rapporto superiore all’unità.

Notare che sarebbe possibile utilizzare anche il reciproco di tale rapporto ma, a quel punto, si deve puntare ad avere un valore inferiore all’unità (complementarità).

Tale concetto è applicabile a più ambiti. Si consideri, ad esempio, un piano inclinato.

La forza instabilizzante Fi è la componente della forza peso parallela alla superficie di scivolamento. Tale forza sarà contrastata dalla forza di attrito che assumerà le vesti di forza resistente Fr.

Tralasciando la notazione vettoriale (essendo la direzione esplicitata dal disegno), note la massa m del corpo e la geometria del problema si ha:

$F_P =mg$

$F_i =F_P sin(\theta)=mg sin(\theta)$

$P =F_P cos(\theta)=mg cos(\theta)$

La forza resistente Fr sarà data dal prodotto tra la componente P (normale alla superficie) e il coefficiente di attrito statico.

$F_r =P \mu =mg \mu cos(\theta)$

La componente P viene bilanciata dalla reazione vincolare offerta dal piano (se la struttura è sufficientemente robusta, altrimente cede). Il moto lungo la direzione tangenziale del piano avviene se sussiste una forza dinamica, ovvero se la differenza tra la forza ressitente e la forza instabilizzante non è nulla. Il corpo “resterà fermo” se:

$F_r =F_i \rightarrow mg\mu cos(\theta) = mg sin(\theta)$

Ovvero l’angolo critico (oltre il quale avviene il moto) sarà quello che soddisfa la seguente condizione:

$\mu cos(\theta) = sin(\theta) \rightarrow \mu = \frac{sin(\theta)}{cos(\theta)} \rightarrow \mu = tg(\theta)$

Rammentando la definizione di Fattore di sicurezza Fs, applicandola al caso in esame, si ha che:

$F_S = \frac{F_R}{F_I}= \frac{mg \mu cos(\theta)}{mg sin(\theta)} = \frac{\mu }{tg(\theta)}$

Si osservi che con Fs = 1 (come accennato sopra, equilibrio perfetto tra le forze) si ottiene nuovamente la relazione che definisce l’angolo critico:

$F_S = \frac{\mu }{tg(\theta)} = 1 \rightarrow \mu = tg(\theta)$

A questo punto è d’uopo un chiarimento. Qualora il fattore di sicurezza fosse maggiore di uno, matematicamente ne risulterebbe che la forza di attrito è maggiore della componente tangenziale. Questo NON corrisponde alla realtà fisica. Se così fosse, il blocco inizierebbe a muoversi verso l’alto. Prima di contattare le autorità ecclesiastiche, si rammenti che la forza di attrito, definita come forza normale per un coefficiente di attrito statico, è quella “massima”. Supponento infatti che l’angolo $\theta$ sia nullo.

$F_r =P \mu =mg \mu cos(\theta) = mg \mu cos(0) = mg \mu$

Per muovere il blocco è necessario comunque applicare una forza dall’esterno che superi F_r, ma se ciò non avviene non si ha moto.

Tale risultato lo si ritrova plottando le componenti della forza peso (supposta, per questo calcolo, unitaria) al variare di theta, scrivendo il seguente codice:

import numpy as np
import matplotlib.pyplot as plt
import math
n= 100 #numero di campioni
us = 0.61 #coefficiente attrito statico
theta = np.linspace(0, math.pi ,n) #angolo theta_ da 0 a 90°
Fp = 1 #Forza peso
#vettori delle componenti al variare di theta
Fi = np.zeros(n)
Fr = np.zeros(n)
#calcolo componenti
for i in range(0,n):
    Fi[i] = Fp * math.sin(theta[i])
    Fr[i] = us*Fp * math.cos(theta[i])
plt.plot(theta,Fi, label="Fi")
plt.plot(theta,Fr, label="Fr")
plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left', ncol=2, mode="expand", borderaxespad=0.)
plt.xlabel('Angolo')
plt.ylabel('Forze')
plt.grid()
plt.show()

Le curve si incrociano proprio nel valore di theta di 0,54rad: oltre il quale, la forza resistente Fr non riuscirà più a contrastare la forza motrice Fi.

Quindi, idealmente, fino a che l’inclinazione non supera tale valore, non si corre il rischio che il blocco scivoli via. Il problama “reale” è che i suddetti valori non sono noti senza incertezza.

Metodo Monte Carlo per la stima del fattore di sicurezza

Fino a che si è in un contesto deterministico, ovvero il caso ideale, quanto detto è sufficiente. Nella realtà si hanno delle “incertezze” dei parametri misurati. Essi infatti non saranno definiti da uno scalare ma da una distribuzione probabilistica.

A seguire un esempio in cui i valori numerici sono scelti di autorità.

Si supponga di disporre un blocco di Alluminio su un piano di acciaio, quindi, idealmente :

$\mu = 0.61$

Da un semplice calcolo si ha che l’angolo limite vale:

$\theta = atan(\mu) = 0,54 rad$

Quindi, per sicurezza,si predispone un blocco che non consenta di andare oltre i 0,5 rad. Bisogna ora verificere se tale “azione per la sicurezza” è sufficiente.

Il problema è che il dispositivo di regolazione dell’inclinazione presenta una tolleranza del 5%. Il piano, inoltre, potrebbe essere sporco, cosa che potrebbe variare il coefficiente di attrito statico del 20%.

Partendo dalle distribuzioni probabilistiche cumulative (C.D.F.) , sia del coefficiente di attrito che dell’angolo:

Dalle quali, partendo dall’estrazione casuale di un valore da una distribuzione probabilistica uniforme (0-1), si estrarranno i due valori necessari. Per ogni estrazione si calcolerà il fattore di sicurezza (relativo all’estrazione i-esmima):

$F_Si = \frac{\mu_i }{tg(\theta_i)}$

Procedendo quindi con un numero di estrazioni pari a n=1000, si “contano” le volte in cui $F_Si$ è minore ad un certo valore limite $F_{min}$ che per questioni cautelative poniamo pari a 1.05, 5% in più del necessario, per poi calcolare la frequenza con cui ciò accade.

Il codice in Python:

import numpy as np
import matplotlib.pyplot as plt
import math 
import random 

#Valori centrali
us = 0.61 #Coefficiente di attrito statico 
theta = 0.5 #angolo di inclinazione [rad]

tol_theta= 0.05 #tolleranza sull'angolo
tol_u = 0.2 #tolleranza sul coefficiente di attrito

theta_min = theta(1-tol_theta) 
theta_max = theta(1+tol_theta)
us_min = us(1-tol_u) 
us_max = us(1+tol_u)

K = 300 #numero di passi
#Parametri distribuzione normale 
sigma_us = 0.025 #Deviazione standard (sqrt)
sigma_theta = 0.007 #Deviazione standard 
(sqrt)


#Funzione per la costruzione della PDF
def PDF(x,sig,xc,k):#(x, sigma, valore centrale, passi)
     a = 1/(pow(2*np.pi*pow(sig,2),0.5))
     pdf = np.zeros(k)
     for i in range(0,k):
         pdf[i] = a*math.exp(-0.5*pow(x[i]-xc,2)/(pow(sig,2)))
     return pdf

#Funzione per la costruzione della CDF
def CDF(pdf,k):
     cdf = np.zeros(k)
     cdf[0] = pdf[0]
     for i in range(1,k):
         cdf[i] = cdf[i-1]+pdf[i]# Normalizzazione C.D.F.
     cdf = cdf /max(cdf)
     return cdf

#Costruzione PDF angolo
vec_theta = np.linspace(theta_min,theta_max,K)#vettore angoli
PDF_theta = PDF(vec_theta,sigma_theta,theta,K)

#Costruzione PDF "u"
vec_us = np.linspace(us_min,us_max,K)#vettore coefficienti
PDF_us = PDF(vec_us,sigma_us,us,K)

#Costruzione CDF 
CDF_theta = CDF(PDF_theta,K)
CDF_us = CDF(PDF_us,K)
#Estrazione valore da DISTRIBUZIONE UNIFORME

#Funzione per l'estrazione (distribuzione uniforme 0-1)
def dado(): return np.random.uniform(0,1)#Random

#Funzione di estrazione
def estrazione(cdf, k):
    r = dado()
    if(r<= cdf[0]):
        r_v = 0 
    for i in range(1,k):
        if(r<= cdf[i] and r > cdf[i-1]):
            r_v = i
     return r_v

 #ESTRAZIONI
 n = 1000 #numero di estrazioni
 Fs_min = 1.05
 f = np.zeros(K)
 for i in range(0,K):
     f[i] = Fs_min*math.tan(vec_theta[i])
 cont = 0
 for i in range(0,n):       
     us_es = estrazione(CDF_us,K)#Estrai u
     theta_es = estrazione(CDF_theta,K)#Estrai theta
     #Plotta il punto (coppia di valori) estratto
     plt.scatter(vec_theta[theta_es],vec_us[us_es],color ="red",s = 1) 
     #Verifica se il Fsi è inferiore al valore limite
     if(vec_us[us_es]/math.tan(vec_theta[theta_es])< Fs_min):
         cont += 1 #Incrementa contatore (Fs <Fs_min)


 plt.plot(vec_theta,f,color ="black") 
 plt.ylabel(" coefficiente di attrito")
 plt.xlabel("Angolo")
 plt.grid()
 plt.show()

 Prob = 100*cont/n #Calcola probabilità
 risultato = "La probabilità che Fs sia inferiore a " + str(Fs_min) + " è del: " + str(Prob) + "%"
 print(risultato)

Il risultato ottenuto:

“La probabilità che Fs sia inferiore a 1.05 è del: 9.0%”

Quindi, seppur l’angolo è stato impostato al di sotto della soglia limite, ovvero quella che “matematicamente” garantiva la staticità del blocco, dal Metodo Monte Carlo si evince che si ha comunque un rischio elevato.

Ciò che si può fare, ora, è ridurre il valore di angolo massimo impostabile. Provando con 0,45 rad (quindi cambiando solo il valore centrale della P.D.F. di partenza):

“La probabilità che Fs sia inferiore a 1.05 è del: 0.1%”

Risultato che consente di lavorare in sicurezza. Si rammenta che, qualora fosse possibile, sia per questioni etiche che per Legge, qualora fosse possibile prendere ulteriori provvedimenti che consentano di ridurre tale probabilità, questi dovranno essere messi in campo.

Si rammenta che il rischio è definito come prodotto tra la magnitudo del “danno” che l’evento nefasto causerebbe e la probabilità che ciò accada. Seppur 0,1% può risultare un valore irrisorio, gli effetti di un blocco di qualche quintale che scivola via sono piuttosto ingenti, pertanto il rischio rimane elevato.

Matteo Gentileschi

Lune Elettriche

"Amo colui che tutto dona e nulla restituisce"

Fattore di sicurezza con il “Metodo Monte Carlo”

Lascia un commento Cancella risposta

Condividi:

Correlati

Lascia un commento Cancella risposta