Robitex's Blog

Ideas in the web

Archivi Categorie: Antisismica

L’Arte della pericolosità sismica – Parte 2


Sommario

Con un nuovo tocco di colore, torniamo sulla mappa della pericolosità sismica italiana per produrre un’altra pagina artistica. Questa volta useremo una scala di colori con quelli della nostra bandiera…

Colore italiano

Riprendendo il codice illustrato in questo post, è sufficiente creare una nuova funzione Lua, riportata in seguito, per ottenere la mappa sottostante:

La mappa di pericolosità sismica in versione tricolore

La mappa di pericolosità sismica in versione tricolore

Il colore, funzione dell’accelerazione, varierà da verde per il valore minimo a bianco per l’accelerazione media, fino ad arrivare al rosso per l’accelerazione massima.
In LaTeX, un colore può essere miscelato facilmente con il bianco specificandone dopo il nome ed un punto esclamativo, un valore variabile da 0 a 100. Il valore nullo corrisponde al colore bianco, mentre il valore 100 lascia intatto il colore.
Così, l’espressione ‘green!50’ corrisponde al colore bianco con il 50% di verde.

Il codice seguente mostra la funzione che genera l’espressione del colore in funzione dell’accelerazione. Valori bassi di pericolosità saranno rappresentati da un colore verde mescolato con sempre meno bianco a mano a mano che l’accelerazione decresce, mentre valori elevati saranno rappresentati da un colore rosso mescolato ad una componente decrescente di bianco:

-- funzione di interpolazione lineare
local function linear(x,x1,x2,y1,y2)
   return y1 + (y2-y1)*(x-x1)/(x2-x1)
end

-- funzione scala di colori dal verde al rosso
-- passando per il bianco
local function green2red(ag)
   local agmean = (agmax+agmin)/2
   local basecol, val
   if ag < agmean then
      basecol = 'green'
      val = linear(ag, agmin, agmean, 100, 0)
   else
      basecol = 'red'
      val = linear(ag, agmean, agmax, 0, 100)
   end
   
   return string.format('%s!%d', basecol, val)
end

Ho stampato anche questa nuova mappa, naturalmente a colori, e l’ho appesa in ufficio a beneficio della collezione primavera estate 2012.
Ciao.

L’Arte della pericolosità sismica


Sommario

Attraverso i software del mondo TeX, creeremo una visualizzazione artistica della mappa di pericolosità sismica italiana contenuta nella norma D.M. 14/01/2008. L’elaborazione grafica sarà opera del potente pacchetto TikZ progettato per LaTeX, mentre un piccolo programma in Lua si occuperà della costruzione logica.

Pericolosità sismica S1

Il progetto S1 di INGV, l’Istituto Nazionale di Geofisica e Vulcanologia, recepito poi dalla normativa tecnica, ha definito la pericolosità sismica nel territorio italiano attraverso una griglia di punti distanti uno dall’altro poco più di 5 km.
L’idea è quella di replicare le rappresentazioni grafiche diffuse dall’INGV con gli strumenti del mondo TeX, realizzando mappe a colori con i valori dell’accelerazione di picco al suolo PGA, Peak Ground Accelaration (un risultato è riportato in fondo all’articolo).

Immagine della mappa S1 tratta dal sito INGV

Immagine della mappa S1 tratta dal sito INGV

Struttura del codice

Ecco una descrizione di come verranno elaborati i dati:

  • uno script principale in Lua, efficiente linguaggio di scripting, leggerà i dati della pericolosità producendone un file di testo contenente comandi grafici per TikZ;
  • un sorgente LaTeX elaborerà i comandi producendo la mappa grafica in formato pdf.

I dati

Ho tradotto la tabella scaricabile dal sito del Consiglio dei Lavori Pubblici, contentente i dati di pericolosità sismica, in una forma direttamente comprensibile da Lua: ecco come appare il contenuto del file di dati che ho chiamato ‘GridS1.lua’:

GridS1 = {
[13111] = {lon=6.5448, lat=45.134, Tr30={0.0263, 2.50, 0.18}, Tr50={0.0340, 2.51, 0.21}, ...ecc
[13333] = {lon=6.5506, lat=45.085, Tr30={0.0264, 2.49, 0.18}, Tr50={0.0341, 2.51, 0.21}, ...ecc
[13555] = {lon=6.5564, lat=45.035, Tr30={0.0264, 2.50, 0.18}, Tr50={0.0340, 2.51, 0.20}, ...ecc
...ecc

Creare i comandi grafici

Disegneremo per ciascun punto della griglia S1, un quadratino colorato secondo una data scala di colori associata al valore della PGA, ma per farlo dobbiamo risolvere due problemi:

  • tradurre la coppia longitudine, latitudine dei punti in coordinate nel piano x, y;
  • associare un colore al valore di accelerazione.

Il sorgente Lua mostrato nel seguito dapprima carica in memoria i dati strutturati della griglia, determina i valori minimi e massini dell’accelerazione tra tutti i 10751 punti (per il tempo di ritorno di 475 anni) utili in seguito, poi definisce le due funzioni fondamentali di trasformazione delle coordinate (trasf()) e di associazione valore->colore (wavecolor).

Preparato il terreno, lo script si conclude con la funzione go(), che scrive su un file esterno chiamato defpoint.tex, le istruzioni TikZ di disegno.

Notate che la funzione go() accetta come argomenti anche delle funzioni. Il vantaggio immediato di questo tipo di programmazione disponibile in Lua è che la funzione go() diviene indipendente sia dal tipo di trasformazione di coordinate che dalla scala di colore (io ho scelto di lavorare con la definizione di colore basato sulla lunghezza d’onda della luce con rapporto lineare rispetto alla PGA).

Per i più curiosi, nei commenti del codice sono riportati i due punti della maglia S1 corrispondenti al punto sismicamente più pericoloso e quello meno pericoloso sul territorio della penisola italiana…

-- carica la tabella dati dal file GridS1.lua
dofile "GridS1.lua"

-- funzione di massimo o minimo valore
local function S1(funval, key)
   local rif = funval(GridS1[13111])
   local id

   for idx, t in pairs(GridS1) do
      local val = funval(t)
      if key == "max" then
         if  val > rif then
            rif = val
            id = idx
         end
      elseif key == "min" then
         if val < rif then
             rif = val
             id = idx
          end
       else
          error("The key '"..key"' is not supported.")
       end
    end
    return rif, id
end
-- funzioni di estrazione dati singolo nodo
local function Lat(t)
    return t.lat
end
local function Lon(t)
    return t.lon
end
local function Ag475(t)
    return t.Tr475[1]
end
local function countNode()
    local i = 0
    for _, t in pairs(GridS1) do
       i = i + 1
    end
    return i
end
-- valori minimi e massimi dei dati di mappa
local agmax, idmax = S1(Ag475, 'max')  --> 0.278 , 49640
local agmin, idmin = S1(Ag475, 'min')  --> 0.0364, 12693

-- latitudine
local latmax = S1(Lat, 'max')   --> 47.178
local latmin = S1(Lat, 'min')   --> 36.573

-- longitudine
local lonmax = S1(Lon, 'max')   --> 18.594
local lonmin = S1(Lon, 'min')   --> 6.5448

-- numero nodi S1 (solo per saperlo...)
local numnode = countNode()     --> 10751

-- funzione di servizio
local function makecmdfile(fn, t)
   local outputfile = assert(io.open(fn..".tex","w"))
   outputfile:write(table.concat(t,"\n"))
   outputfile:close()
end

-- restituisce le coordinate di mappa
-- da quelle geografiche
local function trasf(lon, lat)
   -- fattore di conversione gradi->cm
   local f = 4.1

   return 0.80*lon*f, lat*f
end

-- funzione scala di colori wavelets
local function wavecolor(ag)
   local w1 = 363
   local w2 = 814
   local m = 80
   w1 = w1 + m
   w2 = w2 - m
   return w1 + (w2-w1)*(ag-agmin)/(agmax-agmin)
end

-- maschera di formato disegno
local cmask = '\\convertcolorspec{wave}{%d}{rgb}\\col\\color[rgb]{\\col}'
local dmask = '\\draw (%0.3f,%0.3f) node {};'

-- filename = nome senza estensione del file risultato
-- fcoord = funzione di trasformazione coordinate
-- fcolor = funzione di trasformazione colore
local function go(filename, fcoord, fcolor)
   -- memorizzazione comandi
   local t = {}
   for id, dati in pairs(GridS1) do
      local ag = Ag475(dati)
      local lt = Lat(dati)
      local ln = Lon(dati)

      t[#t+1] = string.format(cmask, fcolor(ag))
      t[#t+1] = string.format(dmask, fcoord(ln,lt))
   end

   makecmdfile(filename, t)
end

-- esecuzione
go('defpoint', trasf, wavecolor)

Creare la mappa

Per creare la mappa è sufficiente compilare con pdfLaTeX questo minuscolo file:

\documentclass{standalone}
\usepackage{tikz}
\usepackage{xcolor}
\tikzstyle{every node}=[
    fill,
    minimum size=0.1pt
]

\begin{document}
\begin{tikzpicture}
\input{defpoint}
\end{tikzpicture}
\end{document}

Risultato

La mappa risultato è questa, insieme di 10751 quadratini colorati in modo diverso uno dall’altro (fate un click sull’immagine per scaricarla, un’esperienza da brivido…):

Mappa della pericolosità sismica italiana

La mappa della pericolosità sismica italiana risultato dell’elaborazione descritta nell’articolo

Note tecniche

Per riprodurre la mappa eseguendo il codice mostrato nel post, occorre disporre di una installazione TeX, come TeX Live o MiKTeX, e l’interprete Lua, tutti strumenti gratuiti, o meglio, liberi.
Avremo potuto lavorare con un unico strumento (Lua incorporata in TeX, ovvero LuaTeX, ma lascio il facile esercizio al lettore che non perderà l’occasione di inventarsi nuove scale di colore… 🙂
Alla prossima!

Il terremoto come fenomeno poissoniano: il tempo di ritorno


Scarica l’articolo in formato PDF per la stampa

Probabilità poissoniana di accadimento

La distribuzione poissoniana delle probabilità ipotizza che se un evento è del tutto casuale, allora la probabilità che accada n volte nel periodo T è data dalla relazione seguente, noto il valore di y, il numero medio di accadimenti nel periodo considerato:

\displaystyle p_{n}=\frac{y^n}{n!}e^{-y}.

Se per esempio un evento casuale capita mediamente una volta l’anno (y=1), la probabilità che capiti due volte l’anno è di soli 18,4%.

La probabilità dipende quindi dall’intervallo di tempo di riferimento V_R, periodo in cui è noto il numero medio di accadimenti dell’evento casuale.
Se dividiamo il numero medio di accadimenti con il periodo ad esso correlato, otteniamo la frequenza media di accadimento dell’evento \lambda.

\displaystyle y=\lambda V_R.

Se invece al contrario, dividiamo il periodo di riferimento con il numero medio di accadimenti, otteniamo l’intervallo di tempo medio tra gli accadimenti dell’evento. Questo intervallo è chiamato periodo di ritorno T_R.
Il periodo di ritorno e la frequenza sono uno l’inverso dell’altro.

Chiediamoci ora la probabilità che l’evento casuale non si verifichi nel periodo di riferimento. La relazione di Poisson per nessun evento diventa (zero fattoriale vale 1):

\displaystyle p_{0}=\frac{(\lambda V_R)^0}{0!}e^{-\lambda V_R}=e^{-\lambda V_R}.

Il complemento ad uno di p_0 è la probabilità che almeno un evento si verifichi nel periodo considerato, in formule:

\displaystyle p_{n \geq 1}=p_{VR}=1-e^{-\lambda V_R}.

Questa relazione ci permette di ottenere la relazione tra il periodo di ritorno e la probabilità di superamento, fissato il periodo di riferimento. Con qualche semplice passaggio:

\displaystyle T_R=-\frac{V_R}{\log(1-p_{VR})}.

Le curve di frequenza/probabilità

Mettiamo su un grafico semi-logaritmico, quello che abbiamo appena trattato con il pacchetto PGFPlots, un componente software utilizzabile all’interno di un documento LaTeX, con l’ausilio di Gnuplot per sopperire alle capacità di calcolo di TeX.
Operativamente, occorre scrivere il testo sintatticamente corretto (il codice) che successivamente viene elaborato da pdflatex per produrre il grafico nel formato pdf (fate click sull’immagine sottostante per scaricare il risultato).

Probabilità poissoniana di un evento puramente casuale

Dal punto di vista della dotazione software, occorre una distribuzione TeX in configurazione abbastanza completa e recente (TeX Live 2010), gnuplot ed (opzionale) la suite Image Magick per tradurre il pdf in un file immagine. Tutti i programmi citati sono liberi e disponibili per vari sistemi operativi, compreso Windows.
Una volta copiato il codice in un file di testo, chiamato per fissare le idee pvr.tex, la sequenza completa dei comandi da terminale sarà:

$ pdflatex -shell-escape pvr
$ pdfcrop --margin=3 pvr.pdf pvr.pdf
$ convert -density 120 pvr.pdf pvr.png

Un breve commento al codice: si può arrivare abbastanza velocemente a scriverlo con un po’ di lavoro di rifinitura come vedete dalla sfilza di opzioni passate all’ambiente semilogxaxis. Una volta che si sono decisi i dettagli del grafico, risulta molto semplice plottare le funzioni. Basta costruire il comando \addplot gnuplot dando le opzioni di identificazione (id) e l’intervallo di dominio e scrivere l’espressione della funzione voluta. Anzi, a ben vedere si tratta di un’ottima interfaccia a gnuplot, per costruire grafici in formato vettoriale (pdf).

\documentclass{minimal}
\usepackage{pgfplots}
\thispagestyle{empty}

\begin{document}
\begin{tikzpicture}
\begin{semilogxaxis}[
    xlabel={Frequenza media nel periodo, $\lambda=1/T_R$},
    ylabel={Probabilit\'a per $n\geq1$, $p_{V_R}$},
    ytick={0,0.2,0.4,0.6,0.8,1},
    yticklabels={0\%,20\%,40\%,60\%,80\%,100\%},
    grid=both,
    xmin=0.001,
    xmax=0.1,
    ymin=0,
    ymax=1,
    no markers,
    line width=1pt,
    width=15cm,
    height=9cm,
    smooth,
    legend pos=north west]
    
\addplot gnuplot[id=pvriv,domain=0.001:0.1] {1-exp(-x*200)};
\addlegendentry{$V_R = 200$ anni};

\addplot gnuplot[id=pvriii,domain=0.001:0.1]{1-exp(-x*100)};
\addlegendentry{$V_R = 100$ anni};

\addplot gnuplot[id=pvrii,domain=0.001:0.1] {1-exp( -x*50)};
\addlegendentry{$V_R = 50$ anni};

\addplot gnuplot[id=pvri,domain=0.001:0.1]  {1-exp( -x*10)};
\addlegendentry{$V_R = 10$ anni};

\end{semilogxaxis}
\end{tikzpicture}
\end{document}

Probabilità dell’evento sismico

Il terremoto non è un evento puramente casuale per il semplice fatto che il susseguirsi degli eventi sismici di una stessa regione, narra la storia evolutiva delle strutture di faglia e quindi non sono eventi indipendenti uno dall’altro ne eventi puramente casuali.

Tuttavia la nuova normativa tecnica D.M. 14 gennaio 2008, fa l’assunzione che la probabilità di accadimento del terremoto sia poissoniana, usufruendo della semplicità matematica del modello ma commettendo un errore considerato comunque accettabile.
Infatti in una stessa regione, per intervalli di tempo dell’ordine di qualche migliaia d’anni, è ragionevole supporre che le faglie attive continuino ad evolversi costantemente nel tempo mantenendo immutate le velocità relative ed i meccanismi d’interazione.
Nel lungo periodo invece i mutamenti di placche e microplacche della crosta comportano che le faglie aumentino o diminuiscano la loro pericolosità. Nel futuro ci attenderà un assetto completamente mutato, ma per raggiungerlo i passi sono talmente piccoli che possiamo approssimarli a nessun movimento.

Allora possiamo rispondere alla domanda: qual’è la probabilità che almeno un evento sismico con periodo di ritorno di 475 anni si verifichi nell’intervallo di 50 anni (intervallo in cui una costruzione può dirsi efficiente)?

\displaystyle p_{VR}=1-e^{-50/475}=0.10

La risposta è quindi il 10% (corrisponde allo Stato limite di salvaguardia della vita per il periodo di riferimento di 50 anni), ed è chiamata nel linguaggio della norma probabilità di superamento (o di eccedenza).
Terremoti più intesi avranno un tempo di ritorno più lungo e di conseguenza probabilità più piccole di verificarsi nel periodo di riferimento, ed è appunto questo il fenomeno che misuriamo con la legge delle probabilità di Poisson.

Conclusioni: le norme tecniche e la ricerca

Le nuove normative tecniche (Eurocodici ed NTC 2008) mettono in risalto ancor di più la tendenza non tanto ad essere prestazionali nell’approccio alla sicurezza, ma soprattutto a discendere direttamente ed esclusivamente dalle trattazioni scientifiche dei problemi. In conseguenza è possibile ritrovare l’origine e le motivazioni nascoste delle prescrizioni e ciò rende il processo normativo maggiormente trasparente, chiarendo quali ipotesi sono state assunte ed il loro contesto di validità.

Con l’avanzamento delle conoscenze, le norme tecniche divengono sempre più aderenti alla complessità dei fenomeni e questo esige una sempre maggiore preparazione di base da parte dei progettisti e, poiché la ricerca scientifica è un fatto internazionale, tendono ad essere sempre più somiglianti una con l’altra.

Dunque l’ampliarsi di quelle che sono le capacità di fare, in altri termini l’ampliarsi dell’ingegneria, induce di conseguenza alla specializzazione dei progettisti, a cui è affidato il difficile compito di mantenere l’equilibrio tra profondità ed ampiezza della conoscenza, ed un uniformarsi delle normative che svolgono una vera e propria selezione naturale dei risultati delle ricerche prodotte nel mondo.

La pericolosità sismica nella rete dell’INGV


Scarica l’articolo in formato PDF per la stampa

Database sismologici aperti

Non tutti sanno che l’Istituto Nazionale di Geofisica e Vulcanologia mette a disposizione sul web non solo le informazioni generali sull’origine di terremoti e vulcani in Italia e nel mondo, ma anche i database sismologici storici ed attuali con i dati provenienti direttamente dalla rete dei sismografi.
Questo post vi aiuterà in qualche modo a consultare queste risorse ed a rendervi meglio conto di cosa siano le conseguenze dei fenomeni tettonici in evoluzione da sempre sul nostro pianeta.

Iside

Connessi alla home page dell’INGV trovate nel menù a sinistra la voce Banche Dati che vi apre lo sconfinato mondo dei dati geofisici più diversi.

Quello che ci interessa è per esempio ricercare i terremoti che si sono verificati in Italia nei primi due giorni di febbraio accedendo ad Iside, l’accesso ai dati strumentali della rete di rilevamento (circa 250 stazioni).
Cliccate in alto alla voce Terremoti->Localizzazioni. Si aprirà il form per eseguire l’interrogazione. Inseriamo l’intervallo di date dal 1/02/2011 al 2/02/2011, inseriamo il codice di controllo antibot (consiglio caldamente che vi registriate sul sito anche per sostenere le attività di INGV oltre che per evitare di reinserire ogni volta il codice) e premiamo il pulsante ‘Cerca’.

Il sistema mi restituisce ben 43 eventi sismici mostrati nella mappa sottostante, di cui solo 4 superano la soglia della magnitudo 2. Clicchiamo il pulsante “Mostra mappa” per constatare che tutti questi piccoli sismi, difficilmente rilevati dalla popolazione, sono accaduti sull’appenino centrale.
Tra l’altro la mappa è disponibile anche in formato vettoriale Postscript, in formato testuale CSV (Comma Separated Value), ed in kml il formato gestito da Google Earth.

Earthquakes from 2011/02/01 to 2011/02/03 in Italy

Earthquakes from 2011/02/01 to 2011/02/03 in Italy

Emidius

Nel sito emidius, sempre di INGV, potete invece consultare il database macrosismico che riporta le registrazioni storiche dei terremoti italiani tratte da documenti d’archivio, e che quindi si basano su testimonianze degli effetti prodotti dai terremoti.

Al solito, facciamo click sul link di menù sulla sinistra chiamato “Consultazione per località”, per elencare i terremoti storicamente registrati nella nostra capitale: Roma.
Selezionando la lettera R in alto e la città di Roma nell’elenco alfabetico che di conseguenza compare a sinistra, si vedrà la tabella con i dati delle registrazioni ed il grafico che vi riporto qui sotto:

Rome earthquake history

Rome earthquake history

Il progetto S1

Il progetto S1 è essenzialmente lo studio sulla pericolosità sismica su cui si fonda la nuova normativa tecnica italiana D.M. 14 gennaio 2008. Se volete visualizzare i dati della vostra aerea collegatevi al sito esse1. Cercate il comune digitandone il nome in basso a destra nel campo “Search Municipality”, poi cambiate la scala nel campo soprastante, diciamo 200.000 e fate click su “Change scale/center”. Poi attivate il basso il flag “Show the grid points with the value of:” e date un “Redraw the map” cliccando sul pulsante con la freccia circolare.
Quello che vi appare qui sotto è l’applicazione del procedimento descritto per il comune di Roma.

Seismic hazard for the Rome city

Seismic hazard for the city of Rome

Adesso prendete in alto a destra della pagina lo strumento d’informazione chiamato “Graph on the grid point” (click sul radio bottom) e selezionate un punto di griglia che compare nella mappa (quadratino colorato). Vi appariranno le curve dell’accelerazione sismica PGA con le probabilità di accadimento, per tre diversi valori del percentile (derivati dall’analisi statistica dei modelli adottati (sedici scenari)).

Infine, selezionate, sempre nello stesso menù, lo strumento “Disaggregation graph” e cliccate sul punto di griglia d’interesse. Appare un grafico che mostra per il sito in esame come si distribuisce la pericolosità con la distanza. Se il grafico mostra zone colorate concentrate entro una breve distanza significa che il sito si trova nei pressi di una faglia attiva, se invece vi sono aeree colorate più ampie, vorrà dire che la pericolosità deriva da regioni sismogenetiche lontane.

INGVterremoti

L’Istituto INGV ha perfino un canale YuoTube dove brevi filmati illustrano la genesi di terremoti e Tsunami ed le ricerche più recenti, come quella del video linkato qui sotto, relativo al sisma dell’Aquila, dove in tempi molto brevi utilizzando dati dai satelliti dei sistemi GPS, non solo si è riusciti a stimare lo spostamento del terreno causato dal sisma dell’Aquila (che si è abbassato di 20 centimetri in accordo al tipo di faglia di tipo normale) ma, soprattutto, si è riusciti a mettere a punto un modello confermato sperimentalmente della faglia di Paganica stessa.

Cosa significa tutto ciò?

La comunicazione scientifica via web è uno strumento importante almeno da tre punti di vista:
1- si tratta di informazioni utili immediatamente consultabili dagli utenti (studenti, ricercatori, amministratori, progettisti, incaricati di protezione civile e semplici cittadini);
2- è un incentivazione al lavoro scientifico stesso affinché non manchino i fondi necessari e non si chiudano le possibilità di attivare nuovi progetti di ricerca;
3- è la partecipazione e collaborazione tra sorgenti del sapere e comunità che sortisce effetti senza dubbio importanti sulle scelte.

Se si parla di rischio sismico, abbiamo gli strumenti per valutare la pericolosità ed intervenire opportunamente in prevenzione sul patrimonio edilizio. Solo in questo modo il paese si svilupperà nel verso giusto avvantaggiandosi opportunamente delle conoscenze scientifiche e tecnologiche invece di ignorarne l’esistenza.

Pensare al pianeta in termini della sua struttura tettonica complessiva è una dimensione nuova che convive con la quotidianità fatta dei problemi dello studio, del lavoro e della famiglia. Forse ci stiamo avvicinando ad essere un umanità mutata che vede il pianeta non più come una sconfinata distesa oltre i dintorni in cui ognuno vive, ma come un luogo dove si è tutto e parte, dove l’orizzonte si sposta fuori da noi stessi per abbracciare un intero pianeta.

Un quadro strutturale


Come volevasi dimostrare

Un’idea è sempre un’idea.
Sto osservando lo schermo che propone viste tridimensionali di una struttura a telaio i cui elementi hanno colori casuali, a volte scelti in automatico senza un criterio dal software di modellazione durante le operazioni di modifica.

Osservo lo schermo e, ad un tratto, mi sembra di vedere una composizione. Eccola (intanto che ammirate vi saluto, ciao):

Quadro strutturale n. 1

Quadro strutturale n. 1

Diagramma d’acciaio


Caratterizzare le prove di trazione

Spesso un diagramma presenta i dati molto più efficacemente dei corrispondenti valori numerici, dunque l’idea di graficare le prove a trazione delle barre in acciaio per armatura.
Si tratta delle prove secondo lo standard UNI EN 10002-1:2004, che si fanno eseguire in un laboratorio autorizzato, per la certificazione dei materiali impiegati nella costruzione della struttura in calcestruzzo armato. Tali prove sono caratterizzate dal valore della tensione di snervamento, dalla percentuale di allungamento a rottura e della tensione di rottura (riferimento norma D.M. 14/01/2008).

Di solito il certificato di prova che giunge dal laboratorio, non include un diagramma della prova ma solo i parametri principali. Se tuttavia assumiamo che l’allungamento allo snervamento abbia un valore dello 0.2%, un diagramma di caratterizzazione della prova si può disegnare congiungendo con linee rette i punti di snervamento e rottura nel classico piano (tensione,deformazione).

Ulteriori informazioni possono essere scaricate dal sito del consorzio Sismic.

A test for reinforced steel bar drawing by the pgfplots package

Un diagramma d’acciaio

Potevamo farci sfuggire l’occasione? Naturalmente no!
Affidiamoci al pacchetto LaTeX di turno, in questo caso il meraviglioso pacchetto pgfplots, per ottenere un diagramma d’acciaio di alta qualità.

Il risultato visibile nell’immagine riportata è scaricabile in formato pdf (fai click su questo link), dimostra che in breve tempo, consultando il manuale di pgfplots, sono stato in grado di ottenere tutte le personalizzazioni del grafico che desideravo.
Ecco il codice LaTeX completo pronto per la compilazione:

\documentclass[a5paper,landscape]{article}
\usepackage[italian]{babel}
\usepackage[T1]{fontenc}

\usepackage{amssymb}
\usepackage{pgfplots}
\usepackage[decimalsymbol=comma]{siunitx}
\usepackage[margin=10mm]{geometry}

\thispagestyle{empty}

\begin{document}

\begin{tikzpicture}
% defining personal drawing style
\pgfplotsset{blue style/.style={line width=1pt,mark=x,color=blue}}
\pgfplotsset{red style/.style={line width=1pt,mark=x,color=red}}
\pgfplotsset{green style/.style={line width=1pt,mark=x,color=green}}

\begin{axis}[legend pos=south east,
             font=\footnotesize,
             ymin=0,ymax=730,xmin=0,xmax=14.3,
             xlabel=Allungamento (\%),ylabel=\si{N/mm^2},
             grid=major,
             width=15cm,height=9cm]
             
             
% diametro 8
\addplot[blue style] coordinates {(0,0)(0.2,584.0)(8.3,658.4)};
\addplot[blue style] coordinates {(0,0)(0.2,569.9)(11.4,656.6)};
\addplot[blue style] coordinates {(0,0)(0.2,573.9)(8.0,657.0)};
%
% diametro 10
\addplot[red style] coordinates {(0,0)(0.2,497.6)(13.3,592.0)};
\addplot[red style] coordinates {(0,0)(0.2,501.0)(11.5,595.2)};
\addplot[red style] coordinates {(0,0)(0.2,501.4)(12.0,597.9)};
%
% diametro 12
\addplot[green style] coordinates {(0,0)(0.2,545.4)(12.5,633.4)};
\addplot[green style] coordinates {(0,0)(0.2,558.0)(9.8,647.7)};
\addplot[green style] coordinates {(0,0)(0.2,558.5)(10.5,651.5)};
%
\legend{\(\varnothing{} 8\) arm. platea,,,
        \(\varnothing{} 10\) arm. platea,,,
        \(\varnothing{} 12\) arm. platea}
\end{axis}
\end{tikzpicture}

\end{document}

Bel risultato vero? Abbiamo fatto uso anche di stili di tracciamento personalizzati, grazie alla modernità del linguaggio di PGF.

Alla prossima prova d’acciaio allora! Ciao.

Promenade structural


Progetto di Caina

In questo breve video potete constatare come la rappresentazione 3D dei modelli strutturali sia un mezzo formidabile per la comprensione della geometria delle opere. Numerose indicazioni possono essere ottenute per risolvere le complesse relazioni tra aste, l’ottimizzazione dei telai, ed i rapporti con gli elementi architettonici.

Alla prossima!

Spettri sismici – Parte 2


Un codice più flessibile

Continuiamo il perfezionamento del codice che avevamo scritto in precedenza per il disegno degli spettri elastici di progetto secondo il capitolo 3 delle nuove norme sulle costruzioni (D.M. 14/01/2008), dipendente dai tre parametri base a_g/g, F_0, T_C^* e relativamente agli stati limite di operatività SLO, di danno SLD, di salvaguardiadella vita SLV e di collasso SLC.

In Asymptote, il software che useremo per svolgere i calcoli e plottare i grafici, esiste la possibilità di creare nuovi tipi di dato che incorporano dati e funzioni, ovvero l’invenzione della programmazione ad oggetti.

Questa possibilità è offerta da struct, il costrutto che assomiglia a quello del C++, che ci permetterà di rappresentare le Categorie di sottosuolo e le categorie topografiche del sito.

E cominciando da quest’ultima. Serve un oggetto che chiameremo TopoCat che memorizzi il valore del coefficiente topografico e che lo restituisca per mezzo di una funzione:

// oggetto per le categorie topografiche
struct TopoCat {
   private real st;
   void operator init(real st){
      this.st = st;
   }
   real St(){
      return this.st;
   }
}

Abbiamo dichiarato come privato il campo st così che solo le funzioni pubbliche ne costituiscano l’interfaccia, ed inserito lo speciale operatore init che definisce il cosidetto costruttore, con cui l’utente creerà nuove istanze dell’oggetto. In questo caso il costruttore inizializza il campo interno accessibile con il riferimento this in dot notation.

Per instanziare un oggetto per ciascuna categoria topografica scriveremo così, creando i quattro tipi classici previsti dalla tabella 3.2.VI della norma:

// categorie topografiche classiche
TopoCat T1 = TopoCat(1.0);
TopoCat T2 = TopoCat(1.2);
TopoCat T3 = TopoCat(1.2);
TopoCat T4 = TopoCat(1.4);

Anche per le categorie di sottosuolo creiamo un oggetto apposito che viene instanziato con i coefficienti della tabella 3.2.V per fornire, sottoforma di funzioni naturalmente, i valori di Ss e Cc.

Da notare come il codice gestisce l’espressione di Ss: invece di riferirsi a campi interni dell’oggetto, i valori di ag/g e F0 devono semplicemente essere passati alle funzioni pubbliche (dette anche ‘metodi’ nel gergo della Object Oriented Programming).
Credo anche che fare la stessa cosa con un foglio di calcolo non sia così immediato.

// Categoria di sottosuolo
struct SoilCat {
   private real minSs;
   private real maxSs;
   private real k1;
   private real k2;
   private real h1;
   private real h2;

   // constructor
   void operator init(real mn, real mx, real k1, real k2, real h1, real h2){
      this.minSs = mn;
      this.maxSs = mx;
      this.k1 = k1;
      this.k2 = k2;
      this.h1 = h1;
      this.h2 = h2;
   }

   real Ss(real agg, real f0){// return the Ss value
      private real s;
      s = this.k1 - this.k2 * f0 * agg;
      if ( s < this.minSs ) {
         s = this.minSs;
      }
      if ( s > this.maxSs ) {
         s = this.maxSs;
      }
      return s;
   }
   
   real Cc(real tc){// return the Cc value
      return this.h1 * tc^this.h2;
   }
}

//
// istanziamento oggetti categorie sottosuolo (Tab. 3.2.V)
//
SoilCat A = SoilCat( 1.00 , 1.00 , 1.00 ,  0.00 , 1.00 ,  0.00);
SoilCat B = SoilCat( 1.00 , 1.20 , 1.40 , -0.40 , 1.10 , -0.20);
SoilCat C = SoilCat( 1.00 , 1.50 , 1.70 , -0.60 , 1.05 , -0.33);
SoilCat D = SoilCat( 0.90 , 1.80 , 2.40 , -1.50 , 1.25 , -0.50);
SoilCat E = SoilCat( 1.00 , 1.60 , 2.00 , -1.10 , 1.15 , -0.40);

Bene, notate come ovunque nel codice ci siano le dichiarazioni di tipo (il tipo corrispondente viene premesso al nome di ogni variabile), come richiede il linguaggio fortemente tipizzato di Asymptote.

Il codice finora scritto consente di rappresentare due diverse caratteristiche del sito. E’ naturale inglobare le due categorie in un ulteriore nuovo oggetto per rappresentare unitariamente il sito ed esporre i metodi per il disegno degli spettri.

Chiamiamo questo oggetto Spettro:

// Struct contenenti i parametri sismici del
// sito per lo Stati Limite considerato
// e funzioni grafiche e numeriche
//
// agrel = ag/g  ,  f0= F0  , tcstar = TC*
//
// q = fattore di struttura
//
struct Spettro {
   private real agrel;
   private real f0;
   private real tcstar;
  
   private real q;
   
   private real eta;
   
   private real smax;
   
   // periodi d'intervallo
   private real tb;
   private real tc;
   private real td;
   
   // ecco quì le caratteristiche del sito
   // rappresentate come campi privati
   private SoilCat soilC;
   private TopoCat topoC;
   
  // array dei punti notevoli
   pair[] keydot;
   
   // funzioni di definizione per intervalli dello spettro
   real s1(real T) {
      return  this.smax * ((T/this.tb) + (1 - T/this.tb)/(this.eta * this.f0));
    }

   real s3(real T){
      return this.smax * this.tc / T;
   }

   real s4(real T){
      return this.smax * (this.tc * this.td ) / T^2;
   }
   
  // costruttore
   void operator init(SoilCat s, TopoCat t, real q, real ag, real f0, real tc){
      this.soilC = s;
      this.topoC = t;
      this.q = q;
      this.agrel = ag;
      this.f0 = f0;
      this.tcstar = tc;

      this.eta = 1/q;

      this.tc = s.Cc(tc)*tc;
      this.tb = this.tc/3;
      this.td = 4.0 * ag + 1.60;

      this.smax = ag * s.Ss(ag,f0) * t.St() * (1/q) * f0;

      this.keydot = new pair[]{
        (      0 , this.s1(0)),
        (this.tb , this.smax),
        (this.tc , this.smax),
        (this.td , this.s4(this.td))
      };
   }
   guide plot(real xmax){
     return graph(this.s1 , 0 , this.tb)                   &
            (this.tb , this.smax) -- (this.tc , this.smax) &
            graph(this.s3 , this.tc , this.td)             &
            graph(this.s4 , this.td , xmax);
   }
}

Un commento anche per questo ultimo codice: la struct inizia con la dichiarazione dei campi privati e prosegue con quella delle funzioni matematiche dei rami dello spettro sismico ricavate dalla norma. La loro definizione deve precedere l’operatore init, poiché il costruttore ne fa uso per inizializzare le coordinate dei punti estremi dei rami dello spettro ai periodi Tb, Tc, e Td, salvandoli in un array chiamato keydot. Infine, il metodo plot che restituisce l’intero “percorso” del grafico.

L’oggetto non memorizza ulteriori metodi per facilitare il disegno dei grafici perché così facendo si lascia totale libertà di configurazione dei dettagli. In fondo la curva dello spettro è un oggetto indipendente dalla “tela” sulla quale verrà plottata.

Il codice “utente”

Ecco come può essere utilizzato il codice appena scritto. Consideriamo un sito su un terreno tipo B in pendio (T2). I quattro spettri elastici (q=1) con i vari parametri ag/g, F0 e Tc* saranno per esempio:

//
// un "oggetto struct" per singolo spettro
// istanziati con i dati sismici del sito
//
Spettro slo = Spettro(B, T2, 1.00, 0.04873, 2.462, 0.239);
Spettro sld = Spettro(B, T2, 1.00, 0.06107, 2.484, 0.251);
Spettro slv = Spettro(B, T2, 1.00, 0.15182, 2.402, 0.290);
Spettro slc = Spettro(B, T2, 1.00, 0.19475, 2.379, 0.298);

Il risultato

Non rimane che eseguire i plottaggi con il seguente codice, per ottenere in formato PDF il grafico completo pronto per essere inserito in relazione:

//
// plot degli spettri
//
size(16cm, 10cm, IgnoreAspect);
xlimits( 0 , 4.32);
ylimits( 0 , 0.72);

xaxis("$T$ (s)" ,
      BottomTop(),
      LeftTicks( Label(fontsize(8)) ,
               ptick=lightgray ,
               pTick=gray ,
               extend = true
               )
      );

yaxis("$S_{e}/g$" ,
      LeftRight(),
      RightTicks( Label(fontsize(8)),
                  pTick=lightgray,
                  extend = true
               )
      );

// plottaggio spettri (l'ordine è dall'alto verso il basso
// per la corrispondenza con la legenda
draw(slc.plot(4.2), red    + 1.80, "SLC");
draw(slv.plot(4.2), blue   + 1.80, "SLV");
draw(sld.plot(4.2), orange + 1.80, "SLD");
draw(slo.plot(4.2), green  + 1.80, "SLO");

// disegno dei punti notevoli
pen punto = linewidth(1.28bp);
dot(slc.keydot, punto);
dot(slv.keydot, punto);
dot(sld.keydot, punto);
dot(slo.keydot, punto);

// legenda
add(legend(), (3.1 , 0.60) , UnFill );
Spettri Sismici, fare click per scaricare il PDF

Spettri Sismici, fare click per scaricare il PDF

Cosa ci aspetta il futuro

Potremo aggiungere al grafico alcune informazioni di input per dare un tocco di professionalità, e pensare ad un modulo che calcoli i parametri ag/g, F0, Tc* del sito inserendone le coordinate geografiche. Chissà staremo a vedere.
Ciao

Spettri sismici Parte 1


Spettri sismici

Seismic Spectrum

Seismic Spectrum

La nuova normativa D.M. 14/01/2008 “Norme Tecniche per le Costruzioni”, in breve NTC 2008, entrata in vigore a metà del 2009 poco dopo il sisma in Abruzzo, ha introdotto una nuova mappatura del rischio sismico basata sul progetto S1 dell’INGV – Istituto Nazionale di Geofisica e Vulcanologia.

La pericolosità sismica dipende ora non più dai confini amministrativi dei comuni, ma dal punto d’interesse stesso, individuato in coordinate geografiche nel datum ED50.

Se anche a voi piace determinare gli spettri per proprio conto vi propongo un percorso tra gli strumenti del mondo (La)TeX, in particolare con Asymptote, con cui ottenere il grafico degli spettri sismici in formato vettoriale, pronto da inserire nella relazione tecnica del progetto strutturale.

Asymptote e gli spettri sismici

Oltre che ad imparare la programmazione con Asymptote (è consigliabile consultare il manuale di Asymptote presente nella vostra installazione o reperibile in rete), costruiremo gli spettri passo passo leggendo direttamente il paragrafo 3.2.3 della NTC 2008.

Questo post introduce l’argomento concentrandosi solo sullo spettro elastico di progetto allo stato limite di salvaguardia della vita per la sola componente orizzontale, e per un terreno di pianura di tipo B. Nei prossimi vedremo come il linguaggio di Asymptote possa implementare una pratica e completa interfaccia ai dati ed alle elaborazioni.

Installazione di Asymptote

Se avete una distribuzione LaTeX recente già installata sul vostro PC, qualsiasi sia il sistema operativo, tutto il software necessario è già pronto per l’uso, altrimenti potete rivolgervi indifferentemente a TeX Live 2009 (per Linux, Mac OS X, Windows, ecc) o MiKTeX 2.8 (per Windows), oppure potete installare esclusivamente Asymptote scaricandolo da questo link ed inserire nella relazione il grafico in uno dei nel formati grafici per esempio gif o png.

Calcolare lo spettro

Cominciamo!!! Avviate un editor di testo (Blocco Note, SciTE, Emacs, Gedit, ecc ecc), e scriviamo il seguente codice per instanziare le prime variabili:

// loading modulo graph
import graph;

// accelerazione di gravità (valore standard) in m/s2
real g = 9.80665;

// coefficiente di amplificazione topografica
real st = 1.0;

// fattore di struttura (spettri elastici)
real q = 1.0;

// eta per spettri di progetto
real eta = 1/q;

La sintassi è simile a Java, con le dichiarazione di tipo (real per le variabili numeriche in virgola mobile), il punto e virgola finale e le righe di commento che iniziano con il doppio slash.

Vediamo ora la prima funzione che restituisce il valore del coefficiente stratigrafico per un suolo di tipo B dalla tabella 3.2.V della NTC:

\displaystyle 1{,}00 \leq 1{,}40 - 0{,}40 F_0 \frac{a_g}{g}\leq 1{,}20

e la sua implementazione in Asymptote:

// funzione del coefficiente di amplificazione
// stratigrafica per il suolo tipo B
real ss( real agsug, real f0 ) {
	real tmp = 1.40 - 0.40 * f0 * agsug;
	if ( tmp < 1.00 ) {
		tmp = 1.0;
	}

	if ( tmp > 1.20 ) {
		tmp = 1.20;
	}
	return tmp;
}

Scriviamo adesso il codice per determinare i periodi caratteristici dello spettro a partire dai parametri fondamentali del sito a_g/g, F_0, T_C^*, per un periodo di ritorno corrispondente allo SLV, e per un sito di esempio:

// parametri fondamentali del sito
// agrel=ag/g, f0=F0, tcstar=Tc*
//
real agrel  = 0.1946;
real f0     = 2.423;
real tcstar = 0.28;

real tc = 1.10 * tcstar^(-0.20) * tcstar;
real tb = tc / 3 ;
real td = 4.0 * agrel + 1.60;

L’espressione dei quattro tratti dello spettro sono le seguenti (3.2.3.2.1 NTC 2008):

\displaystyle S_1 = a_g S \eta F_0 \left(\frac{T}{T_B}+\frac{1}{\eta F_0}\left(1-\frac{1}{T_B}\right)\right)

\displaystyle S_2 = a_g S \eta F_0

\displaystyle S_3 = a_g S \eta F_0\frac{T_c}{T}

\displaystyle S_4 = a_g S \eta F_0 \frac{T_cT_D}{T^2}

assegnamo quindi altrettante funzioni che restituiscono il valore dello spettro rapportato all’accelerazione di gravità:

// costante moltiplicativa
real s0 = agrel * ss(agrel, f0) * st * f0;

// primo tratto dello spettro
real s1(real T) {
  return  s0 * ((T/tb) + 1/(eta*f0) * (1 - (T/tb)));
}

// secondo tratto
real s2(real T) {
   return s0;
}

// terzo tratto
real s3(real T) {
   return s0 * tc / T;
}

// quarto tratto
real s4(real T) {
   return s0 * tc * td / T^2;
}

Graficare lo spettro

Tutto è pronto per disegnare il grafico. Basta scrivere una funzione chiamata appunto plot, che restituisce un oggetto di tipo guide, ovvero una curva a sua volta restituita dalla funzione di libreria graph. Il simbolo & è l’operatore di concatenazione delle singole curve, mentre il tratto orizzontale dello spettro è stato ottenuto tracciando direttamente un segmento tra il punto (tb, s0) ed il punto (tc, s0), utilizzando il simbolo -- secondo una sintassi molto comune tra i pacchetti grafici in LaTeX. La funzione graph accetta tre argomenti: il primo è la funzione che determina la curva mentre il secondo ed il terzo, delimitano l’intervallo di tracciamento sull’asse x.
La nostra funzione plot, accetta un numero inteso come limite complessivo del grafico sull’asse x, in modo da poter controllare la “coda”:

guide plot(real xmax){
	return 	graph(s1, 0, tb) &
		(tb, s0) -- (tc, s0) &
		graph(s3, tc, td)    &
		graph(s4, td, xmax);
}

Prepariamo l’area del grafico con l’ausilio di una libreria apposita integrata in Asymptote chiamata graph (modulo spiegato in dettaglio nel manuale di Asymptote):

// dimensione del grafico
size(16cm, 10cm, IgnoreAspect);

xlimits( 0 , 4.32);
ylimits( 0 , 0.70);

// etichettatura assi e caratterizzazione griglia
xaxis("$T$ (s)",
	BottomTop(),
	LeftTicks( Label(fontsize(8)),
		ptick=lightgray,
		pTick=gray,
		extend = true
	)
		);
yaxis("$S_{e}/g$" ,
	LeftRight(),
	RightTicks( Label(fontsize(8)),
		pTick=lightgray,
		extend = true
		)
	);

Plottiamo lo spettro passando alla funzione draw la nostra curva ed aggiungiamo una legenda:

// plottaggio spettro
draw(plot(4.2) , blue + 1.80 , "SLV");

// aggiungi legenda
add(legend(), (2.8 , 0.6) , UnFill );

Da notare come il linguaggio di Asymptote gestisca le proprietà di spessore e colore del tratto, semplicemente sommando il nome del colore con il valore dello spessore. Forse non sarà così elegante e qualche programmatore Java storcerà il naso, ma è senza dubbio comodo.
Tutte le funzioni sono spiegate nel manuale di Asymptote, in particolare la funzione add() e la funzione legend() utilizzate nell’ultima parte del codice.

Confezioniamo il tutto

Se avete seguito fin qui il codice digitato nel vostro editor dovrebbe essere questo:

// loading modulo graph
import graph;

// accelerazione di gravità (valore standard) in m/s2
real g = 9.80665;

// coefficiente di amplificazione topografica
real st = 1.0;

// fattore di struttura (spettri elastici)
real q = 1.0;

// eta per spettri di progetto
real eta = 1/q;

// funzione del coefficiente di amplificazione
// stratigrafica per il suolo tipo B
real ss( real agsug, real f0 ) {
	real tmp = 1.40 - 0.40 * f0 * agsug;
	if ( tmp < 1.00 ) {
		tmp = 1.0;
	}

	if ( tmp > 1.20 ) {
		tmp = 1.20;
	}
	return tmp;
}

//
// parametri fondamentali del sito
// agrel=ag/g, f0=F0, tcstar=Tc*
//
real agrel  = 0.1946;
real f0     = 2.423;
real tcstar = 0.28;

real tc = 1.10 * tcstar^(-0.20) * tcstar;
real tb = tc / 3 ;
real td = 4.0 * agrel + 1.60;

// costante moltiplicativa
real s0 = agrel * ss(agrel, f0) * st * f0;

// primo tratto dello spettro
real s1(real T) {
  return  s0 * ((T/tb) + 1/(eta*f0) * (1 - (T/tb)));
}

// secondo tratto
real s2(real T) {
   return s0;
}

// terzo tratto
real s3(real T) {
   return s0 * tc / T;
}

// quarto tratto
real s4(real T) {
   return s0 * tc * td / T^2;
}

guide plot(real xmax){
	return 	graph(s1, 0, tb)     &
		(tb, s0) -- (tc, s0) &
		graph(s3, tc, td)    &
		graph(s4, td, xmax);
}

// dimensione del grafico
size(16cm, 10cm, IgnoreAspect);

xlimits( 0 , 4.32);
ylimits( 0 , 0.70);

// etichettatura assi e caratterizzazione griglia
xaxis("$T$ (s)",
	BottomTop(),
	LeftTicks( Label(fontsize(8)),
		ptick=lightgray,
		pTick=gray,
		extend = true
	)
		);
yaxis("$S_{e}/g$" ,
	LeftRight(),
	RightTicks( Label(fontsize(8)),
		pTick=lightgray,
		extend = true
		)
	);

// plottaggio spettro
draw(plot(4.2) , blue + 1.80 , "SLV");

// aggiungi legenda
add(legend(), (2.8 , 0.6) , UnFill );

Bene. Salvato il file con l’estensione *.asy per esempio con il nome spettro.asy, basta aprire una finestra di comando (in gergo *nix il Terminale) e scrivere il comando di compilazione (spostatevi con il comando cd nella directory dove risiede il file):

asy -f pdf spettro.asy

Se non avviene automaticamente, aprite il file pdf di output e controllate il risultato. Confrontatelo anche con questo compilato da me.

Per oggi chiudiamo qui. La prossima volta ci occuperemo di generalizzare per le altre situazioni il codice con il paradigma di programmazione object oriented, concettualizzando lo spettro in un oggetto.
Alla prossima.

%d blogger hanno fatto clic su Mi Piace per questo: