La simulazione (parte I)

Eccoci di nuovo qui dopo una breve pausa (causa esame di Sistemi di Controllo (mentale della popolazione) a Più Variabili). Nell'ultimo articolo relativo alla modellistica del comportamento di una particella singola avevamo alla fine di tutto elaborato un modello completo del comportamento di una particella immersa in un fluido incomprimibile con un flusso laminare.
Ora non resta che mettere insieme i risultati ricavati ed esposti nei vari articoli, e infilarli tutti nel modellino originario per potersene avvalere in fase di simulazione.







Ricordiamo per praticità che il modello in questione è definito dalla seguente:

m_{p}\frac{dv_{z}}{dt}=\frac{6\pi \mu(T(z)) R_{eq}\chi}{C_{c}(R_{eq})}\cdot (u_{z}-v_{z})+
V_{p}\rho(z,T(z))\frac{d u_{z}}{dt}+\frac{V_{p}}{2}\rho(z,T(z))\left(\frac{d u_{z}}{dt}-\frac{d v_{z}}{dt}\right)+

+6R_{eq}^{2}\sqrt{\pi\mu(T(z)) \rho(z,T(z))}\int\limits_{0}^{t}{\frac{\frac{du_{z}}{d\tau}-\frac{dv_{z}}{d\tau}}{\sqrt{t-\tau}}d\tau}+\sum\limits_{i}Fe_{i}

con μ viscosità del mezzo, ρ densità del mezzo, Cc Slip Correction Factor, Req il raggio equivalente della particella, mp massa della particella (considerando il galleggiamento), Vp volume della particella, uz velocità del vento lungo l'asse z, vz velocità della particella lungo l'asse z, T temperatura del mezzo e χ il fattore dinamico di forma.
Arrivati a questo punto, abbiamo precedentemente ricavato ed enunziato ogni relazione possibile per poter calcolare i parametri della formula a partire da pressione, temperatura ed umidità relativa. Ci si basa su queste grandezze perchè sono quelle di maggior interesse fornite dalle misurazioni delle sonde metereologiche, dal mommento che l'intenzione è quella di operare una simulazione numerica in un ambiente che abbia dei profili di pressione, umidità relativa e temperatura reali e non approssimati ai minimi termini (anche se le cose variano veramente di poco).
Iniziando dai parametri più semplici, vediamo che per ottenere la viscosità μ, basta usare la relazione di Sutherland:

\mu = 18.27\cdot 10^{-8}\frac{411.15}{T(z)+120}\sqrt{\left(\frac{T(z)}{291.15}\right)^{3}}

Qui non c'è molto da aggiungere, dal momento che basta sostituire la temperatura all'altezza z, e si ricava banalmente la viscosità dell'aria all'altrzza z.
Per la densità c'è da macchinare un poco con il vapore e le formule di Van der Waals, infatti occorre risolvere le seguenti due equazioni di terzo grado (indipendenti) rispettivamente per ρa e ρv:

\left(p(z)-e(z,T)\cdot u(z)+a\frac{\rho_{a}^{2}}{M^{2}}\right)\cdot(1-\frac{\rho_{a}}{M}\cdot b)=\rho_{a}\frac{R}{M}T

con b = 3.64·10^(-5)[m³/mol] ed a = 0.001358 [hPa (m^6)/mol²] e:

\left(e(z,T)\cdot u(z)+a\frac{\rho_{v}^{2}}{M^{2}}\right)\cdot(1-\frac{\rho_{v}}{M}\cdot b)=\rho_{v}\frac{R}{M}T

con b = 30.52·10^(-6)[m³/mol] ed a = 0.005536 [hPa (m^6)/mol²] ed e(z,T) definito dal modello di Murhpy-Koop:

\log(100\cdot e) = 54.842763 - \frac{6763.22}{T(z)} - 4.21\cdot\log(T(z))+ 0.000367T(z) + \tanh{(0.0415(T(z) - 218.8))}\cdot(53.878 - \frac{1331.22}{T(z)} - 9.44523\cdot\log(T(z)) + 0.014025 T(z))

che formulato in questo modo fornisce la pressione di saturazione del vapore sull'acqua in [hPa] a partire da una temperatura in [K].
La densità cercata, quella da mettere nell'equazione differenziale per il moto della particella sarà data dalla somma dei sue risultati trovati (ρa + ρv = ρ(z)).
Fin qui nessun problema, nè numerico, nè analitico (le formule esplicite per le radici di un'equazione di terzo grado esistono anche se sono un po' complicate). Unica accortezza per chi volesse usare le formule esplicite: innanzitutto assicuratevi di usare la soluzione giusta (quella positiva), secondo tutto attenzione all'ordine in cui mettete le numerose operazioni (potreste aver bisogno di rigirarvi un po' la formula) cercando di minimizzare la propagazione dell'errore. In questo caso basta un metodo di Newton-Raphson non ammortizzato per avere una buona e veloce convergenza.

Ora ultimo ma non per importanza il Cc. Dall'analisi del moto Browniano di una particella, abbiamo ricavato in quest'articolo una relazione tra lo Slip Correction Factor ed il cammino libero medio:

\lambda_{p} = \frac{C_{c}}{6\mu}\sqrt{\frac{2\rho k T R_{p}}{3}}

Ancora prima però in quest'altro articolo avevamo proposto una formulazione dello Slip Correction Factor in funzione del cammino libero medio:

C_{c}=1+\frac{\lambda}{R_{p}}\left[1.257+0.4\cdot e^{-\frac{1.1R_{p}}{\lambda}}\right]

Sostituento la prima relazione nella seconda otteniamo:

C_{c}=1+\frac{\frac{C_{c}}{6\mu}\sqrt{\frac{2\rho k T R_{p}}{3}}}{R_{p}}\left[1.257+0.4\cdot e^{-\frac{1.1R_{p}}{\frac{C_{c}}{6\mu}\sqrt{\frac{2\rho k T R_{p}}{3}}}}\right]

Ora questa brutta bestia di equazione in Cc sfido chiunque a risolversela a mano ed a tirare fuori un'espressione esplicita per Cc. Fortuna nostra che tutti o quasi abbiamo un pc che fa il lavoro per noi in modo numerico. Ora questa stronza di equazione se non si prendono le dovute precauzioni manda in pappa (non lo fa convergere dove vogliamo noi) il nostro adorato metodo di Newton-Rahpson. Innanzitutto se lo si fa girare con un risolutore commerciale (MATLAB, Mpale, Mathematica) occorre stringere assai le tolleranze rispetto ai valori di default: per esempio per MATLAB qualcosa come 10^(-16) rispetto al 10^(-6) di default. Sostanzialmente questi numeretti dicono una serie di cose all'algoritmo:
1) quanto deve essere piccola la funzione per poter dire che vale effettivamente zero e terminare la ricerca dello zero della funzione;
2) quanto deve essere piccola la variazione per passo dell'algoritmo, dello zero che si sta trovando, per poter dire che ormai ci si sta stabilizzando ad un valore, e che quindi non vale la pena di continuare ad iterare, perchè tanto il valore di cui già si dispone non cambierà più di tanto.
Ricordo per semplicità che il metodo di Newton-Raphson classico per la risoluzione del problema f(x) = 0, è impostato come segue:

x(0) = x_{0}

x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})}

dove f'(xn) è la derivata rispetto a x di f(x), calcolata in xn, ed x0 è ovviamente è la stima iniziale dello zero della funzione, per il contro caso del Cc è buona cosa metterlo pari ad 1.
Altra accortezza (che però non necessariamente è da mettere in pratica) è la valutazione dello Jacobiano (derivata della funzione se stiamo nel caso monodimensionale): spesso i metodi iterativi per la risoluzione di equazioni, per andare un po' più veloci e fare meno conti risparmiano sul calcolo della derivata. Insomma per più di qualche passo molti algoritmi sostanzialmente non ricalcolano la derivata della funzione e continuano ad usare il valore trovato qualche iterazione fà. Questa però è un'arma a doppio taglio:
1) Dal momento che non dobbiamo far girare il risolutore da solo fine a sè stesso ma dobbiamo integrare il tutto in un risolutore numerico di equazioni differenziali (un sacco di altre iterazioni e conti), per ogni passo di integratore numerico, occorrerà fare x passi del metodo di Newton, quindi se facciamo meno calcoli per ogni passo, alla fine la simulazione ci metterà molto di meno, quindi è buona cosa (se si può) evitare di calcolarsi la derivata ad ogni passo.
2) Se non calcoliamo la derivata ad ogni passo potremmo perdere in precisione, cosa questa, non certo auspicabile.
Pertanto occorrerebbe trovare un buon compromesso tra le due specifiche contrastanti, specialmente se il risolutore numerico utilizzato è già un po' pesantello di suo.
Vediamo un po' ora che cosa succede setralasciamo gradienti di pressione e tutte le altre schifezze che ci sono nel modello e ci limitiamo a simulare:

\frac{dv_{z}}{dt}=-\frac{(6\pi\mu \cdot R_{eq}\Chi)}{(C_{c} V_{p} (\rho_{p}-\rho_{a}))}\cdot (v_{z}-u_{z})-g;

In questo caso si è scelto l'asse z positivo verso l'alto (opposto rispetto alla formulazione originaria, ma per rimediare basta cambiare il segno a v ed a dv/dt)
con condizioni iniziali v0=-60 [m/s] (verso il basso), z0=5000 [m] ed uz=6 [m/s].
Uniche due accortezze prima di lanciarci nella simulazione:
1) Usare un metodo di integrazione implicito (se usate MATLAB ode23tb, o un'altro risolutore per problemi stiff) come per esempio i vari Runge-Kutta impliciti (quadrature di Radau e Lobatto) o eulero implicito, o metodo dei trapezi implicito, perchè a causa della grande variazione della velocità iniziale in tempi brevissimi, i metodi espliciti (per chi usa MATLAB ode45 et similia) con passo di integrazione adattativo, tendono ad inchiodare velocemente ed all'inizio, e risultano quindi inutilmente lenti.
2) Riportare l'equazione differenziale del secondo ordine ad un sistema del primo ordine, facendo la sostituzione dz/dt = v, per ottenere invece della prima equazione, le due equazioni seguenti:

\left\{ \begin{array}{l}
\frac{\displaystyle dv_{z}}{\displaystyle dt}=-\frac{\displaystyle (6\pi\mu(z) \cdot R_{eq}\Chi)}{\displaystyle (C_{c}(z) V_{p} (\rho_{p}-\rho_{a}(z)))}\cdot (v_{z}-u_{z})-g\\
\\
\frac{\displaystyle dz}{\displaystyle dt}=v_{z}
\end{array} \right.

Ora mettendo dentro vari numeretti, ci possiamo fare le nostre simulazioni (tralasciando per ora la forza di Basset, il gradiente di pressione e la massa apparente della particella) ottenendo i seguenti risultati:
Per le velocità :

e per le posizioni:

Con gli zoom delle regioni di interesse per le due grandezze:
Per le velocità:

e per le posizioni:

Come si può vedere dalle legende, sono stati simulati i comportamenti di varie schifezze, le cui particelle ovviamente tutte con dimensione intorno al nanometro. La sigla HDPE (Hihg Density PolyEthilene) sicuramente non è nuova a chi abbia letto due cose sull'argomento Morgellons.
Le dimensioni sono state scelte come estremamente piccole (sostanzialmente con un fibra dello spessore del nm, possiamo infilzare una particella virale, come una cigliegia con un ago) appunto per minimizzare il fenomeno della separazione della scia in varie fasi. Insomma se spargiamo delle schifezze, ci dobbiamo assicurare che queste con il vento non si allontanino troppo dalla scia di vapore per non dare l'effetto a scia multi-stratificata che si vedrebbe da terra, e si mascherino bene tra le particelle di ghiaccio e di vapore. Come si vede appunto dal grafico, la dispersione è minima, anche se non nulla. Pertanto occorrono delle particelle non troppo grandi (altrimenti non rimarrebbero in aria per troppo poco tempo contrastando con l'ipotesi della durata eccessiva delle scie) nè troppo piccole (altrimenti se ne andrebbero in giro con il vento senza mai cadere a terra e provocando lo spiacevole e vistoso effetto della scia multi-stratificata). Come avremo modo di dimostrare però in seguito, risolvendo le equazioni di condensazione, queste dimensioni però stridono con l'ipotesi di rilevamento a terra di ingenti quantità di particelle nanometriche dopo tempi di permanenza in aria compatibili con l'ipotesi della grande durata delle scie chimiche (questo è dovuto al fatto che la condensazione e la flocculazione tendono a spostare a destra e comprimere la distribuzione di raggio delle particelle in funzione del loro numero con l'avanzare del tempo, facendo quindi sì che dopo un po' non si trovino quasi più particelle al di sotto di un certo diametro).
Nella prossima parte invece vedremo di condurre delle simulazioni più dettagliate considerando tutti i termini che compaiono nell'equazione e vedremo come fare ad impostare la risoluzione di un'equazione integro-differenziale che per giunta ha pure un kernel di integrazione assai brutto (guardate un po' che succede se tau è pari proprio a t, all'estremo superiore di integrazione).

Nota: nel mettere i dati prelevati dalla sonda nel modello, occorre fare attenzione alle unità di misura: la temperatura va convertita da [C°] in [K], l'umidità relativa va divisa per 100, la pressione e le altezze, non vanno modificate e possono essere usate tranquillamente nei modelli così come sono.

Appleman v2.1

Ordunque, nell'articolo Appleman 2.0 abbiamo sostanzialmente delineato i tratti per poter fare una previsione delle scie di condensazione permanenti (temperatura minore della temperatura critica). In questa sede invece evidenzieremo un'aspetto che riguarda la previsione di formazione di scie permanenti che si espandono fino a diventare strutture cirriformi e faremo alcune considerazioni circa la sensibilità della temperatura critica alle perturbazioni.

Quello che si è evidenziato nel precedente articolo è il fatto che uno degli ingredienti fondamentali per la costruzione del grafico delle temperature critiche è la pressione di saturazione del vapore sull'acqua. Ora affinchè la scia invece di "evaporare" o sublimare (passaggio dallo stato solido a quello gassoso) formi una bella nuvoletta di cristalli di ghiaccio occorre che l'umidità relativa presente in atmosfera sia superiore all'umidità relativa rispetto al ghiaccio.
Cosa sono queste grandezze ?
La prima è una legge che descrive l'equilibrio dinamico di un'interfaccia ghiaccio vapore (pressione di saturazione del vapore rispetto al ghiaccio). La seconda è una formula che definisce l'umidità relativa rispetto al ghiaccio.
Usiamo ancora una volta il modello di Murphy-Koop per la pressione di saturazione del vapore sul ghiaccio:

e_{i}(T) = \frac{1}{100}e^{(9.550426 - \frac{5723.265}{T} + 3.53068\cdot \log(T) - 0.00728332\cdot T)}

ove ei è la pressione di saturazione del vapore sul ghiaccio in [hPa] e T è la temperatura in [K].
A quasto punto, l' RH (umidità relativa) e la RHi (umidità relativa sul ghiaccio) che è definita in analogia all'umidità relativa normale come:

RH = \frac{e(T)}{e_{w}(T)}

RH_{i} = \frac{e(T)}{e_{i}(T)}

con, RH ed RHi appartenenti all'intervallo [0,1], e(T) la pressione parziale del vapore, ew(T) pressione di saturazione del vapore sull'acqua ed ei(T) pressione di saturazione del vapore sul ghiaccio.
A questo punto apparentemente non si sa come venirne fuori, perchè provando ad imporre la condizione RH > RHi questa è sempre banalmente falsa, dal momento che ei(T) < ew(T) per ogni T.
Ci viene in aiuto una grandezza che abbiamo già visto in precedenza e che ci consente di determinare direttamente la pressione parziale del vapore e(T) direttamente dai dati: il mixing ratio.
Ricordiamo che il mixing ratio r [g/Kg] è legato alla pressione p [hPa] (mb) e ad e(T) con la seguente formula:

r = 621.97\frac{e(T)}{p-e(T)}

che rigirata con un "ardito" passaggio matematico diventa:

\frac{r}{621.97}(p-e(T)) = e(T)

da cui segue che:

e(T)= \frac{rp}{621.97+r}

sostituendo e(T) così trovato nell'espressione per RHi, questa diventa:

RH_{i} = \frac{rp}{e_{i}(T)(621.97+r)}

A questo punto il gioco è fatto: ei è disponibile dal modello come funzione della temperatura, r è fornito dai dati della sonda, così come la temperatura e la pressione.
I dati con i quali le seguenti analisi vengono condotte sono relativi al giorno 11 settembre 2008 ore 12Z (+2 GMT) prelevati dalla stazione di osservazione di Pratica di Mare (RM).
Ecco un grafico che mostra l'evoluzione dell'umidità relativa non interpolata (punti blu) e RHi (punti viola):

Ecco invece cosa succede alla temperatura critica (blu) ed alla temperatura effettiva (viola):

Ora imponendo le condizioni (T < Tc) ed RHi > RH si evidenziano i seguenti intervalli di altezze:

Ora, l'area evidenziata in verde è relativa agli intervalli in cui l'umidità relativa è maggiore dell'umidità relativa sul ghiaccio (e pertanto se in questa regione c'è formazione di contrails, le contrails saranno persistenti e daranno origine a cirri), mentre l'area evidenziata in blu è relativa all'intervallo di altezze a cui si prevede la formazione di contrail persistenti. Come si può vedere c'è un'area abbastanza ampia (dai 12000 ai 18000 m circa) in cui si prevede formazione di contrail che daranno origine alle nuvole.
Anche oggi niente da fare per le contrail a 3000 m. :D

Analizzato brevemente questo completamento alla teoria di Appleman finora trattata, passiamo ad analizzare un'altro aspetto importante: la sensibilità alla variazione dei parametri. Ora, come spesso mi è stato fatto notare, le misurazioni delle sonde non sono del tutto affidabili. Non lo sono sia perchè la sonda non si muove solo lungo l'asse z ma in tutte e tre le dimensioni (anche presenti moti presenti di spin e nutazione, per il modello di una sonda meteorologica completo bisognerà aspettare un poco, ma mi impegnerò a firnire pure quello se fate i bravi ;) ), ma sopratutto perchè a temperature abbastanza basse la sonda commette degli errori rilevanti sul calcolo dell'umidità relativa, che insieme alla pressione è un'ingrediente fondamentale per il calcolo della temperatura critica. Ora vediamo un po' cosa succede se perturbiamo l'umidità relativa di 0.1 (10%). Normalmente si sarebbe fatta la derivata dell'espressione dispetto a tale parametro, ma dal momento che non vi è presente una formulazione esplicita ci si deve accontentare (almeno in prima battuta) di considerare semplicemente il ∆Tc che si viene a creare a causa di una variazione del 10% (incremento) dell'umidità relativa alla pressione di 1000 hPa:

Cosa vediamo dal grafico ? Abbastanza semplice: che finchè stiamo al di sotto del 60% di umidità relativa, tutto va relativamente bene. Ora se la sonda misura un'umidità relativa alta, in un punto a bassa temperature (dove si suppone appunto che le sonde si lascino fregare più facilmente) c'è da stare attenti, perchè la crescita dell'errore è praticamente esponenziale, con un picco massimo di ±4.7 C° ... insomma attenti a dove facciamo le misurazioni.
Circa gli errori sulla pressione la situazione è molto più tranquilla. Per esempio vediamo un po' cosa succede calcolando lo stesso errore con la stessa perturbazioe di prima per pressioni differenti:

Si evince una scarsa dipendenza dell'errore dalla pressione (rimane praticamente costante lungo l'asse della pressione) che invece viene dominato dall'andamento dell'umidità relativa.

Morale della favola: attenti a dove la sonda ha fatto le misurazioni (nello spazio (RH,T) temperature ed umidità), ma non preoccupatevi eccessivamente di errori sulla pressione, che tanto anche se sgarra un poco non cambia molto. Molto più pericolosa l'incertezza sull'umidità relativa.

Sono un particella di Bario ... (parte III)

"Iuuuu ?!? C'è Nessuno ?!?!"
[Particella di Bario in una scia chimica.]

Eccoci di nuovo qui. Arrivati finalmente alla terza ed ultima parte del trattamento dei modelli matematici per le particelle singole (nelle prossime considerazioni la particella di bario non sarà più sola, ma avrà un infinito non numerabile di amici/amiche).

Ordunque, cosa resta da trattare:
1) I diametri equivalenti;
2) La viscosità dell'aria in funzione della temperatura;
3) La densità dell'aria in funzione della temperatura;
4) L'equazione generale per la dinamica di una particella singola.

Quando si analizza il particolato, oltre a fornirne la distribuzione (per esempio dei diametri o delle masse in funzione del numero di particelle) occorre anche tenere in considerazione il raggio o il diametro equivalente. Insomma finora abbiamo visto che il raggio equivalente Rp compare quà e là nelle formule con una certa insistenza. Faccio notare per l'ennesima volta, che mai uno sciachimista mi abbia saputo fornire i numeri di Knudsen, la granulometria, le distribuzioni di qualsivoglia grandezza in funzione del numero di particelle, nè tantomeno la composizione. Mai inoltre ho visto delle considerazioni che possano essere chiamate "plausibili" (pretendere che siano rigorose è troppo). Quanto alla teoria lo scopo di questo blog è appunto informare: fornire gli strumenti necessari per fare delle considerazioni basilari, sviluppare di mano propria modelli semplici e iniziare a comprendere modelli più complicati che vengono spesso usati nei simulatori dell'evoluzione dei fenomeni atmosferici.
Dopo questa breve digressione però torniamo ai diametri (raggi) equivalenti.
Il diametro equivalente più semplice è senz'altro il diametro volumetrico equivalente, definito come:

D_{ve}=\sqrt[3]{\frac{6}{\pi}V_{p}}

che altro non è, che il diametro di una particella sferica avente lo stesso volume della particella non sferica presa in esame. Vp è ovviamente il volume della particella in esame. Ancora una volta, il discorso sulle semplici particelle sferiche è troppo approssimativo, quindi tanto per cambiare occorre introdurre un fattore di correzione alle varie formule: il fattore di forma dinamico (dynamic shape factor). Questo, altro non è che il rapporto tra la forza FD di Stokes misurata che agisce sulla particella e la forza FDve di Stokes che agirebbe su una particella avente diametro volumetrico equivalente:

\chi=\frac{F_{D}}{F_{D}^{ve}}

Ovviamente per le forme irregolari questo fattore è maggiore di 1.0 (ovviamente per piccoli numeri di Reynolds). Il numero di Reynolds è sostanzialmente un numerello che compare nelle equazioni di Navier-Stokes quando queste vengono adimensionalizzate, ed esprime il rapporto tra le forze di inerzia e le forze viscose (ovviamente rispetto all'aria in flusso laminare, e le particelle in esame i numeri di raynolds sono sufficientemente bassi). Guarda caso per particelle di quarzo (altro fondamentale ingrediente delle scie chimiche, noto anche come biossido di silicio, fondamentale componente degli aerosol desertici, volgarmente noto a tutti come "sabbia") il valore è stato misurato e controllato: χ=1.36, per le fibre invece (altro ingrediente fondamentale delle scie chimiche) che abbiano un rapporto lunghezza-diametro > 20, χ=1.06 (il fatto che sia più piccolo è dovuto al fatto che le fibre tendono ad orientarsi lungo il flusso del fluido), e per un cilindro (rigido) che abbia un rapporto lunghezza-diametro parti a 10, il fattore di forma dinamico è pari a χ=1.43.
Ora come introdurre questo fattore nelle formule? Nulla di più semplice: l'espressione della forza di Stokes che è stata usata nelle prime equazioni è:

F_{s}=\frac{6\pi \mu R_{p}}{C_{c}}\cdot u_{rel}

Ebbene, risolvendo la banale equazioncina per il fattore di correzione dinamico di forma possiamo ottenere la forza effettiva in funzione della forza teorica e del fattore di correzione. Quindi invece di usare il termine

\frac{6\pi \mu R_{p}}{C_{c}}\cdot (u_{z}-v_{z})

nell'equazione differenziale per il moto della particella useremo il termine:

\frac{6\pi \mu R_{ve}\chi}{C_{c}(R_{ve})}\cdot (u_{z}-v_{z})

dove Rve è il raggio volumetrico equivalente della particella.
Ora però manca una cosa, si è considerata la stessa velocità terminale (nell'aria quiescente), ma non la stessa densità. Quale sarebbe invece il diametro equivalente di una particella sferica che ha la stessa densità e la stessa velocità terminale della particella non sferica in esame? Il diametro di Stokes, definito come:

D_{st}=\sqrt{\frac{18 v_{t}\mu}{(\rho_{p}-\rho)g C_{c}(D_{st})}}

dove vt è la velocità terminale della particella. Si noti che a differenza della formulazione classica si è ritenuto più opportuno usare non la densità effettiva della particella, ma tenere in conto anche le proprietà di galleggiamento. Ora la fregatura di questa formula, è che per calcolare il diametro di Stokes, dobbiamo ogni volta risolvere un'equazione non lineare (data la brutta forma di Cc, a meno che qesto non risulti prossimo ad 1).
Il diametro di Stokes è legato al diametro volumetrico equivalente (raggio volumetrico equivalente moltiplicato per 2) dalla seguente relazione:

D_{st}=D_{ve}\sqrt{\frac{C_{c}(D_{ve})}{\chi C_{c}(D_{st})}}

Quello che in letteratura è chiamato Diametro Aerodinamico Equivalente altro non è che il diametro di Stokes, per il quale si suppone che la particella abbia la stessa densità dell'acqua.

Una volta trattati i diametri, passiamo a considerare la viscosità, che anche ricorre nelle formule abbastanza spesso. Fortunatamente la cosa è stata già studiata molto tempo fà, da un tale Sutherland, che ha tirato fuori una relazione per la viscosità dell'aria in funzione della temperatura:

\mu = 18.27\cdot 10^{-8}\frac{411.15}{T(z)+120}\sqrt{\left(\frac{T(z)}{291.15}\right)^{3}}

dove T(z) è la temperatura in Kelvin (ovviamente dipendente dall'altezza z), μ la viscosità dell'aria in hPa·s (millibar per secondo).

Per quanto riguarda la densità, ci sono sostanzialmente tre strade possibili:
1) Mettersi a combattere con l'equazione di stato di Lennard-Jones
2) Usare l'equazione di Van der Waals
3) Usare l'equazione di stato per i gas ideali

Quanto alla prima, la escludiamo perchè essendo il modello già approssimato non ha molto senso mettersi ad elaborare formulazioni super precise che poi comunque avrebbero effetti marginali
L'ultima è un po' troppo banale, ma comunque papabile, per la sua semplicità:

\rho = \frac{p(z)\cdot M}{R\cdot T}

dove p(z) è la pressione in Pa, M è la massa molare dell'atmosfera terrestre 0.0289644 [Kg/mol], T è la temperatura [K], R costante universale dei gas pari a 8.314472 [Pa·m³/(mol·K)], e ρ è la densità [Kg/m³].
Ovviamente se la pressione è p hPa, basterà mettere nella formula 100·p invece che p ed i conti tornano nuovamente.
La seconda possibilità (usare la formula di Van der Waals) è un buon compromesso tra complessità e precisione.
Consideriamo l'equazione di stato di Van der Waals:

\left(p(z)+a\frac{n^{2}}{V^{2}}\right)\cdot(V-n\cdot b)=nRT

ora dal momento che interessa la densità, bisogna ricorrere ad un'altra formuletta per poterla tirare fuori:

\frac{n}{V}=\frac{1}{V_{m}}=\frac{\rho}{M}

Usando questa espressione nell'equazione e mettendo in evidenza una V nella seconda parentesi del membro sinistro si perviene all'equazione:

\left(p(z)+a\frac{\rho^{2}}{M^{2}}\right)\cdot(1-\frac{\rho}{M}\cdot b)=\rho\frac{R}{M}T

A questo punto però siccome vogliamo modellare condizioni atmosferiche realistiche, dobbiamo considerare un paio di cosette: innanzitutto dobbiamo anche considerare il fatto che nell'aria non c'è solo aria, c'è anche il vapore ed in secondo luogo il fatto che la densità complessiva sarà data dalla somma di due contributi: ρa e ρv (rispettivamente densità dell'aria e del vapore), che soddisferanno due equazioni di Van der Waals leggermente differenti. Quanto alle pressioni, applicanto la legge di Dalton sull epressioni parziali, possiamo calcolarle tutte, dal momento che e(z,T) si può calcolare e p(z) è nota (fornita dalle misurazioni atmosferiche). L'equazione per l'aria:

\left(p(z)-e(z,T)+a\frac{\rho_{a}^{2}}{M^{2}}\right)\cdot(1-\frac{\rho_{a}}{M}\cdot b)=\rho_{a}\frac{R}{M}T

con b = 3.64·10^(-5)[m³/mol] ed a = 0.001358 [hPa (m^6)/mol²] e l'equazione per il vapore:

\left(e(z,T)+a\frac{\rho_{v}^{2}}{M^{2}}\right)\cdot(1-\frac{\rho_{v}}{M}\cdot b)=\rho_{v}\frac{R}{M}T

con b = 30.52·10^(-6)[m³/mol] ed a = 0.005536 [hPa (m^6)/mol²].
Ora tutte che tutte le costanti sono state definite, non resta che vedere cos'è la variabile e(z,T). Questa però altro non è che la pressione del vapore, definibile attraverso l'umidità relativa u (nell'intervallo [0,1]) e la pressione di saturazione del vapore attraverso la relazione:

e(T)=e_{s}(T)\cdot u(z)

ora dove es è la pressione di saturazione del vapore che abbiamo già incontrato nell'articolo intitolato Appleman v2.0 e per questa possiamo usare tranquillamente uno dei due modelli proposti.
Pertanto le due equazione per l'aria e per il vapore rispettivamente diventano:

\left(p(z)-e_{s}(T)\cdot u(z)+a\frac{\rho_{a}^{2}}{M^{2}}\right)\cdot(1-\frac{\rho_{a}}{M}\cdot b)=\rho_{a}\frac{R}{M}T

e

\left(e_{s}(T)\cdot u(z)+a\frac{\rho_{v}^{2}}{M^{2}}\right)\cdot(1-\frac{\rho_{v}}{M}\cdot b)=\rho_{v}\frac{R}{M}T

Per trovare la densità complessiva ρ(z,T) non resta che sostituire i numeretti e risolvere le due equazioni. Una volta ottenuti i risultati:

\rho(z, T) = \rho_{a}(z, T)+\rho_{v}(z,T)

Ora per il profilo di temperature, si possono benissimo prendere i dati di una sonda (come abbiamo già fatto, in "Appleman v2.0"), interpolarli ed usarli insieme a tutta la paccottiglia che abbiamo fin qui derivato nell'equazione differenziale. Stesso discorso vale per i dati relativi alla pressione ed all'umidità relativa.
Unica accortezza in questo caso è interpolarli in modo furbo, in modo tale da non creare discontinuità e punti di non derivabilità.

Ora però, siccome ho promesso di fornire l'equazione generale per il moto di una particella, eccovi accontentati:

m_{p}\frac{dv_{z}}{dt}=\frac{6\pi \mu(T(z)) R_{eq}\chi}{C_{c}(R_{eq})}\cdot (u_{z}-v_{z})+
V_{p}\rho(z,T(z))\frac{d u_{z}}{dt}+\frac{V_{p}}{2}\rho(z,T(z))\left(\frac{d u_{z}}{dt}-\frac{d v_{z}}{dt}\right)+

+6R_{eq}^{2}\sqrt{\pi\mu(T(z)) \rho(z,T(z))}\int\limits_{0}^{t}{\frac{\frac{du_{z}}{d\tau}-\frac{dv_{z}}{d\tau}}{\sqrt{t-\tau}}d\tau}+\sum\limits_{i}Fe_{i}

con Vp volume della particella, ρ somma delle soluzioni delle due equazioni di Van der Waals scritte prima, μ data dall'equazione di Sutherland ed il resto come prima.
Ora il secondo termine serve per tener conto del gradiente di pressione intorno alla particella dovuto all'accelerazione del fluido lungo il bordo della particella in questione. Il terzo termine è la forza necessaria per accelerare la massa apparente della particella relativamente al fluido. L'integrale brutto è l'integrale di Basset, che serve a tener conto delle forze che si creano a causa del fatto che il moto della particella perturba lo stato di quiete del fluido. Ed infine la sommatoria altro non è che la somma di tutte le forze esterne agenti sul sistema (campi elettrici, gravità, etc), nel caso si consideri solo la gravità, il termine diventa soltanto mp·g.
Ovviamente, il termine 2, 3 e 4 in genere vengono omessi per semplificare il modello.

Nei prossimi articoli, vedremo delle simulazioni con un modellino un po' più serio, ed inizieremo a formulare dei modelli continui basati sulle PDE, avvalendoci dei risultati ottenuti dall'analisi del moto di una particella singola.

Sono una particella di Bario ... (parte II)

"Iuuuu ?!? C'è Nessuno ?!?!"
[Particella di Bario in moto Browniano in una scia chimica.]

Una critica che spesso viene mossa ai modellini dagli sciachimisti (in quei pochi casi in cui riescano a muovere delle critiche che non siano "Il tuo modello matematico non mi piace!") consiste nell'asserire che siano sostanzialmente deterministici e che quindi non vadano bene per modellare fenomeni "caotici" (quando loro stessi non hanno la più pallida idea di come e perchè si definisca "caotico" un sistema). Ecco dunque che spinti da questa critica e dall'esigenza di calcolare il mean free path per una particella di aerosol/particolato dobbiamo iniziare a tenere conto anche di fenomeni aleatori come il moto Browniano.
Quello che finora ho sentito circa le dimensioni del particolato è più o meno classificabile come:
1) Piccolissimo
2) Al di sotto del millimetro (Dr.s Staninger)
3) Nanoparticelle
etc.
Nessuna indicazione precisa circa le dimensioni. Insomma una nanoparticella per essere definita tale ha un diametro di 1-100 nm (nanometri). A queste dimensioni occorre iniziare a tener conto dello Slip Correction Factor, perchè ci si avvicina al regime di transizione (Numero Kn di Knudsen pari a circa 1). Considerando che l'aria è composta prevalentemente di azoto e che questo ha un raggio di Van der Waals pari a circa 155 pm (picometri, 10^(-12) m), e che le particelle in esame hanno un raggio equivalente di circa 1 nm, c'è un rapporto di quasi uno a mille ed continuare ad usare l'approssimazione continua senza le dovute correzioni può essere rischioso.
Lo scopo che ci si prefigge nel condurre questa analisi è la "complicazione" della formula che è stata derivata nella prima parte, e la derivazione di una relazione per il calcolo del cammino libero medio e del coefficiente di diffusione da usare in formule più complicate, che modellino la diffusione con il laplaciano invece che con il moto Browniano e descrizioni dettagliate delle interazioni di ogni singola particella con le particelle vicine.

Per raggiungere i nostri malefici scopi (ovviamente, disinformare, confondere e oscurare, nonchè conquistare il mondo) ci serviamo di una semplicissima equazione:

m_{p}\frac{d\bold{v}}{dt}=-\frac{6\pi\mu R_{p}}{C_{c}}\bold{v}+m_{p}\cdot \bold{a}

con Rp il raggio equivalente della particella, mp massa della particella, a accelerazione Browniana, Cc lo slip correction factor, μ viscosità del mezzo in cui la particella è immersa ed ovviamente v la velocità della particella nello spazio.
L'equazione in esame inoltre è ben conosciuta nell'ambito della fisica statistica, tanto conosciuta da avere un nome proprio di equazione: Equazione di Langevin.
A prima vista sembra un'equazione tranquillissima, unica cosa è che le variabili in neretto (v ed a), sono in realtà delle variabili aleatorie, e l'equazione in questione in realtà un'equazione differenziale stocastica. Nulla di complicato (almeno concettualmente parlando) dal momento che l'accelerazione accelerazione Browniana è una variabile "casuale" a media nulla, che appunto contribuisce ad un moto discontinuo e casuale della particella in esame. Ora dal momento che essere spinti in una direzione è equiprobabile all'essere spinti nella direzione opposta, si suppone che vi sia la proprietà di isotropia del movimento, ossia che l'accelerazione non abbia una direzione preferita.
Forti di queste conoscenze, introduciamo una variabile ausiliaria r che ci dia la posizione della particella. Dal momento che il moto della particella in questione è isotropo < r > = 0 (anche r avrà media nulla, sempre dal discorso sul fatto che il moto Browniano non abbia una direzione preferita).
Ricordo che la media di una funzione f di una variabile aleatoria x continua (o equivalentemente il suo valore atteso) è definita come:

\langle f(x)\rangle=\int{f(x)\cdot P(x) dx}

Ora, dividendo l'equazione iniziale per la massa e calcolando il valore atteso delle variabili, si ottiene

\frac{d}{dt}\langle\bold{r}\cdot\bold{v}\rangle=-\frac{1}{\tau}\langle\bold{r}\cdot\bold{v}\rangle+\frac{3kT}{m_{p}}

Dove τ è definito come:

\tau =\frac{C_{c}m_{p}}{6\pi\mu R_{p}}

Integrando l'equazione normalmente (nella variabile < r·v >), si ottiene:

\langle\bold{r}\cdot\bold{v}\rangle=\frac{3kT\tau}{m_{p}}+\langle\bold{r}_{0}\cdot\bold{v}_{0}\rangle e^{-\frac{t}{\tau}}

ora notando che:

\langle\bold{r}\cdot\bold{v}\rangle=\langle\bold{r}\cdot\frac{\bold{r}}{dt}\rangle=\frac{1}{2}\frac{d}{dt}\langle r^{2}\rangle

l'equazione finale diventa:

\frac{1}{2}\frac{d}{dt}\langle r^{2}\rangle=\frac{3kT\tau}{m_{p}}+\langle\bold{r}_{0}\cdot\bold{v}_{0}\rangle e^{-\frac{t}{\tau}}

Asintoticamente (ma abbiamo visto che il tempo effettivo di settling delle particelle piccole alla velocità dinale è davvero molto breve) possiamo allegramente omettere il termine esponenziale ed integrare in dt, ottenendo:

\langle r^{2}\rangle=\frac{3kT\tau}{m_{p}}t=\frac{kTC_{c}}{\pi\mu R_{p}}t

Sempre per isotropia abbiamo che:

\frac{1}{3}\langle r^{2}\rangle=\langle x^{2}\rangle=\langle y^{2}\rangle=\langle z^{2}\rangle

e pertanto:

\langle x^{2}\rangle=\langle y^{2}\rangle=\langle z^{2}\rangle=\frac{kTC_{c}}{3\pi\mu R_{p}}t

Questo risultato, a cui Einstein è pervenuto per altra via, è stato confermato anche sperimentalmente. Ne consegue che la distanza attraversata quadratica media della particella è direttamente proporzionale al tempo per cui ha subito il moto browniano.
Ora però cosa centra tutto ciò con il coefficiente di diffusione ?
Proviamo a scrivere l'equazione di diffusione semplice in 3 dimensioni:

\frac{\partial N(x,y,z,t)}{\partial t}=D\cdot\nabla^{2}N(x,y,z,t)

dove l'operatore laplaciano di una funzione è definito come:

\nabla^{2}f(x,y,z,t) = \frac{\partial^{2} f(x,y,z,t)}{\partial x^{2}}+\frac{\partial^{2} f(x,y,z,t)}{\partial y^{2}}+\frac{\partial^{2} f(x,y,z,t)}{\partial z^{2}}

mentre la funzione N fornisce il numero di particelle che si sta muovendo di moto Browniano. Possiamo adesso considerare la diffusività dovuta al moto Browniano come un fenomeno macroscopico. Supponiamo ora di trovarci con N0 particelle nel piano (y,z) e supponiamo inoltre che N, non dipenda da y o z.
moltiplicando ambo i membri dell'equazione precedente per x² e facendone l'integrale in dx da -∞ a +∞, otteniamo:

\int\limits_{-\infty}^{+\infty}x^{2}\frac{\partial N}{\partial t}dx = \int\limits_{-\infty}^{+\infty}x^{2}D\frac{\partial^{2} N}{\partial x^{2}}dx

ora, l'integrale a primo membro dell'uguaglianza altro non è che praticamente il valore atteso di x² (una volta portato fuori il differenziale rispetto a t), dal momento che siamo partiti con N0 particelle, con N0 particelle rimaniamo, pertanto:

\int\limits_{-\infty}^{+\infty}x^{2}\frac{\partial N}{\partial t}dx = \frac{\partial}{\partial t}\int\limits_{-\infty}^{+\infty}x^{2}Ndx = N_{0}\frac{\partial \langle x^{2}\rangle}{\partial t}

Ora il secondo integrale è un po' più impicciato da risolvere, ma se ne viene vuori con relativa semplicità usando l'integrazione per parti, e ricordando che:

\int g(x)\frac{d^{2} f(x)}{d x^{2}}dx = g(x)\frac{df(x)}{dx}-\left(\frac{dg(x)}{dx}f(x)-\int f(x)\frac{d^{2} g(x)}{d x^{2}}dx \right)

Applicando questa semplice regoletta al nostro integrale vien fuori che:

\int\limits_{-\infty}^{+\infty}x^{2}D\frac{\partial^{2} N}{\partial x^{2}}dx = D\left(\left[x^{2}\frac{\partial{N}}{\partial x}\right]_{-\infty}^{+\infty} - \left[2xN\right]_{-\infty}^{+\infty}+\int\limits_{-\infty}^{+\infty}2Ndx\right)

i primi due termini sono nulli, mentre l'ultimo termine con l'integrale è pari proprio a 2N0.
Quindi ora avendo sviluppato gli integrali possiamo eguagliare i due risultati, ottenendo:

\frac{\partial{\langle x^{2}\rangle}}{\partial t} = 2D

integrando in dt, otteniamo:

\langle x^{2}\rangle = 2Dt

Uguagliando questa relazione con la relazione ottenuta dal modello iniziale con l'equazione differenziale stocastica, otteniamo quella che in letteratura si chiama relazione di Stokes-Einstein-Sutherland (a meno del fattore Cc):

D=\frac{kT C_{c}}{6\pi\mu R_{p}}

ove k è la costante di Boltzmann, Rp è il raggio equivalente della particella, Cc lo slip correction factor, T temperatura e μ viscosità del mezzo in cui è immersa la particella.
Ora finalmente abbiamo ottenuto un coefficiente di diffusione realistico da poter usare nei futuri modelli alle derivate parziali.
Ora tanto per avere un'idea delle forze che agiscono su particelle estremamente piccole, basta sostituire dei numeretti nelle relazioni precedenti e calcolare la distanza percorsa, per diffusione o per settling gravitazionale, ebbene per particelle con diametro equivalente di 0.01 μm lo spostamento dovuto alla diffusione è circa 1000 volte più grande dello spostamento dovuto alla forza di gravità. Tuttavia facendosi due conti, ci si rende conto che per particelle che non siano di dimensioni comparabili con le molecole del mezzo in cui sono immerse, la diffusione non dà certo un gran contributo al trasporto.
Supponendo che l'ipotesi sciachimista sia valida, notiamo subito che c'è qualcosa che potrebbe stonare: se troppo piccole le particelle saranno soggette di meno alla forza di gravità e più trasportabili per diffusione o per convezione, se troppo grandi, cadranno a terra. Ora la seconda ipotesi è in contrasto con l'ipotesi sciachimista per via della durata, la prima anche, perchè per essere di qualche utilità a terra le particelle non devono stare in quota troppo a lungo, altrimenti poi le andiamo a raccogliere in africa. Insomma le particelle dovrebbero stare in quota un tempo sufficientemente lungo (più di qualche ora) per essere il modello compatibile con l'ipotesi sciachimista, ma non troppo a lungo in modo da non ricadere a troppa distanza dal punto dove si voleva fare irrorazione e non finire in africa. Inoltre sviluppando un modello alle PDE ho notato che l'ipotesi potrebbe avere almeno altri due bachi pericolosissimi, che sostanzialmente smontano completamente la teoria (uno per esempio è il contrasto tra il rilevamento a terra di particelle nanometriche ed i processi di condensazione-nucleazione, l'altro è l'inconsistenza della durata delle scie, il rilevamento di particelle nanometriche a terra ed appunto i processi di condensazione; ma di questo parleremo nella prossima puntata).

Ora però due parole sul "mean free path". Il movimento di una particella è caratterizzato dalla velocità termica media:

\bar{c}_{p}=\sqrt{\frac{8kT}{\pi m_{p}}}

La teoria cinetica dei gas ci dice che questa può essere legata alla diffusività come segue:

D=\frac{1}{2}\bar{c}_{p}\lambda_{p}

usando la relazione di Stokes-Einstein-Sutherland si ottiene la seguente espressione per il cammino libero medio di una particella di particolato in un gas:

\lambda_{p} = \frac{C_{c}}{6\mu}\sqrt{\frac{2\rho k T R_{p}}{3}}

Con l'ovvio significato dei parametri.
Abbiamo quindi ricavato un modo per determinare il cammino libero medio ed il coefficiente di diffusione.

Nella terza parte, tratteremo più in dettaglio la viscosità, i diametri equivalenti e la relazione dinamica generale per una particella singola.

LHC (Large Hadron Collider) Avviato !

Eccoci finalmente al 10 settembre 2008. E cosa sarà mai successo di così importante oggi ? Ebbene, questa mattina alle 10.30 circa ora italiana, al CERN di Ginevra (Svizzera) nell' LHC finalmente sono stati fatti girare i primi fasci di particelle. Tanto per cambiare in barba a tutte le previsioni apocalittiche, non è successo nulla di strano, nessun buco nero, nessuna terra risucchiata (ma d'altronde non siamo ancora al 2012 ;) ), nessun universo distrutto. Capisco che i dati relativi all'LHC possano incutere un po' di riverenziale rispetto, ma non c'è nulla di cui preoccuparsi, nonostante sia una struttura sotterranea a forma di anello dalla lunghezza complessiva di 27 Km, composta prevalentemente di magneti superconduttori, con diverse strutture di accelerazione per aumentare l'energia delle particelle che vengono accelerate all'interno dell'anello.

In sostanza, dentro l'acceleratore due fasci di particelle viaggiano a velocità prossime a quelle della luce e con energie estremamente elevate, prima di essere fatti collidere l'uno con l'altro.
In via teorica con questo acceleratore sarà possibile arrivare ad energie di circa 7 TeV (Tera-Elettronvolt). Detto così non si capisce molto bene di cosa si sta parlando, ed sembra che ci sia qualcosa di strano se si considera che 7 TeV sono equivalenti a qualcosa come 1.1214 *10^(-6) J (Joule). Ora un Joule è pari a 1 N*m (Newton per metro), e praticamente è il lavoro che bisogna compiere contro la forza di gravità per sollevare una mela di circa 102 g di un metro. Detta così sembrerebbe un'energia ridicolamente bassa, ma quest'energia viene applicata ad una particella elementare di massa ancora più ridicola; ergo se si fa una semplice proporzione con un'ogetto delle dimensioni di un millimetro, la quantità di energia per unità di massa che percepirebbe l'ogetto in questione sarebbe un numero con un sacco di zeri. I due fasci presenti nell'acceleratore viaggiano ovviamente in direzioni opposte in due condutture separate, entrambe mantenute ad una condizione di vuoto ultra-elevato (circa 10^(-7) Pa). I fasci in questione vengono guidati lungo la struttura da un potente campo magnetico generato da magneti superconduttori: praticamente delle potentissime elettrocalamite costisuite da avvolgimenti di cavi speciali che cengono fatti operare in regime di superconduttività (quindi quasi senza dissipazione di energia ed una resistenza quasi nulla), a temperature che si aggirano intorno agli 0.15 K (-273 C°). Praticamente dopo l'antenna gravitazionale all'INFN di Frascati (RM) e quella del CERN (Ginevra, Svizzera) l'LHC è uno dei posti più freddi dell'intero universo. Per mantenere temperature così basse l'intero impianto è raffreddato ad elio liquido (eghm si, proprio quello che viene usato per gonfiare i palloncini alle feste, solo un tantinello più freddo), e richiede un'infrastruttura non indifferente.
Tanto per citare qualche dato: per il confinamento dei fasci, vengono usate 1232 magneti bipolari lunghi 15 m circa usati per piegare il fascio, e 392 magneti quadripolari lunghi 5-7 m usati per la focalizzazione dei fasci. Oltre a questi vengono usati anche altri magneti speciali per restringere il facio prima della collisione, insomma le particelle sono estremamente piccole e più il fascio è stretto più ci sono probabilità di avere delle collisioni, tanto per dare un'idea quello che viene fatto nell'acceleratore equivale a sparare degli aghi da due posizioni distanti 10 Km con una precisione tale da farli incontrare a metà strada.
Insomma un lavoro eccezionale, di importanza altrettanto eccezionale, perchè i risultati di queste ricerche potrebbero finalmente contribuire in modo significativo alla creazione di una teoria unificata (ricordo che la fisica quantistica fà letteralmente a botte con la Relatività Generale), ed eventualmente confermare o smentire la Teoria delle Stringhe.

Riferimenti:
Sito ufficiale dell' LHC

LHC News

WorldWide LHC Computing Grid

Sito personale di Marco Delmastro (fisico impegnato presso l' LHC)

Sono una particella di Bario ... (parte I)

"Iuuuu ?!? C'è Nessuno ?!?!"
[Particella di Bario in una scia chimica.]

Sono sicuro che molti di voi ricorderanno la famosa pubblicità di dell'acqua oligominerale a basso contenuto di sodio, dove una simpatica particella si aggirava per la bottiglia, cercando in vano suoi simili.
Cosa c'entra tutto questo con le scie chimiche ?
Partiamo dall'inizio del problema. Più volte è stato affermato che le scie chimiche contengano di tutto e di più (carciofini e tonno, criptonite verde, capperi, peli di castoro, nanomunghi™ etc ... ). Insomma sono piene di schifezze e questo più o meno lo avevamo capito tutti. Ebbene, tra le sostanze più gettonate ci sono alcuni metalli, Bario, Alluminio, Silicio etc.
Supponiamo per una volta che l'ipotesi sciachimista sia vera, ossia che le particelle in questione effettivamente siano presenti all'interno della scia e che la scia effettivamente contenga del Bario.
Iniziamo a creare un semplicissimo modellino LTI (lineare, tempo invariante) del comportamento si una particella nell'aria, ed in particolare l'evoluzione della sua velocità lungo l'asse z.
Innanzitutto, prendiamo l'asse z positivo verso il basso, ed iniziamo a scrivere l'equazione del comportamento della particella lungo l'asse z, avvalendoci dell'equilibrio delle forze:

m_{p}\frac{dv_{z}}{dt}=m_{p}\cdot g + a\cdot(u_{z}-v_{z})

il primo termine è dovuto alla gravità (F = m*a), il secondo alla differenza tra la velocità vz della particella e la velocità uz del fluido lungo l'asse z (velocità relativa).
Ora cos'è questo secondo termine misterioso ? Semplice: la forza di Stokes (forza di attrito viscoso), corredata da un fattore di correzione.
Tale forza ha espressione:

F_{s}=\frac{6\pi \mu R_{p}}{C_{c}}\cdot u_{rel}

dove μ è la viscosità dell'aria, Rp è il raggio equivalente (o fisico/effettivo) della particella, e Cc è lo "Slip Correction Factor". Lo slip correction factor viene applicato quando il numero di Knudsen diventa relativamente alto (per particelle molto piccole).
In particolare il numero di Knudsen Kn è definito come:

K_{n}=\frac{\lambda}{R_{p}}

ove λ è il "mean free path" (cammino libero medio), ossia la distanza media che la particella di particolato può percorrere in linea retta prima di beccare una particella del fluido in cui è immersa. Insomma per particelle grandi (per le quali cioè ha senso fare un'approssimazione continua del mezzo in cui sono immerse, e per cui la legge di Stokes predice bene il comportamento) il numero di Knudsen tende a 0. Per particelle piccole invece, le cui dimensioni sono paragonabili alle dimensioni delle particelle del fluido in cui sono immerse il numero di Knudsen tende ad 1. In quest'ultimo caso infatti la legge di Stokes (derivata da una legge di conservazione della massa e da un'equazione di conservazione dei momenti, appunto equazione di Navier-Stokes) necessita di una correzione per poter fornire dei risultati soddisfacenti. Da cui la necessità di introdurre lo slip correction factor nella formula di prima.
Cc risulta infatti essere definito come:

C_{c}=1+\frac{\lambda}{R_{p}}\left[1.257+0.4\cdot e^{-\frac{1.1R_{p}}{\lambda}}\right]

Pertanto ora che abbiamo praticamente definito tutto quello che c'è da definire, l'equazione differenziale lineare e stazionaria per la velocità che abbiamo considerato prima diventa:

m_{p}\frac{dv_{z}}{dt}=m_{p}\cdot g + \frac{6\pi\mu R_{p}}{C_{c}}\cdot(u_{z}-v_{z})

dividiamo tutto per la massa mp e otteniamo:

\frac{dv_{z}}{dt}=g + \frac{6\pi\mu R_{p}}{C_{c}\cdot m_{p}}\cdot(u_{z}-v_{z})

I più accorti a questo punto quantomeno protesteranno sul fatto che assumere costante la viscosità ed il mean free path, nonchè la massa ed il raggio sia un'eresia; infatti queste grandezze dipendono dalla temperatura ed un'altra serie di parametri che sostanzialmente possono essere ricondotti alla dipendenza dall'altezza a cui si trova la particella. Mi sento di rassicurare costoro sul fatto che le "dovute complicazioni" verranno introdotte più avanti.
La soluzione all'equazione precedente si ottiene in modo banale (se non ci fosse di mezzo un'integrale di convoluzione sarebbe ancora più facile) nel modo seguente:

v_{z}(t)=e^{-\frac{6\pi\mu R_{p}}{C_{c}\cdot m_{p}}\cdot t}\cdot v_{z}(0) + \int\limits_{0}^{t}
e^{-\frac{6\pi\mu R_{p}}{C_{c}\cdot m_{p}}\cdot (t-\tau)}\frac{6\pi\mu R_{p}}{C_{c}\cdot m_{p}}\cdot u_{z}(\tau)d\tau

Se si assume la velocità u, costante e non dipendente dal tempo, la soluzione si semplifica in:

v_{z}(t)=\frac{C_{c}g\cdot m_{p}+6\pi R_{p}\mu\cdot u_{z}+e^{-\frac{6\pi\mu R_{p}}{C_{c}\cdot m_{p}}\cdot t}\cdot (6\pi R_{p}\mu (v_{z}(0)-u_{z})-C_{c}g\cdot m_{p})}{6\pi R_{p} \mu}

A questo punto si può pure assumere che la particella in questione si sferica, e scrivere la massa mp tenendo conto delle proprietà di galleggiamento come:

m_{p}=\frac{4}{3}\pi R_{p}^{3}(\rho_{p}-\rho)

dove ρp è la densità della particella, e ρ è la densità del fluido in cui la particella è immersa (ovviamente anche questa densità dipenderà dalla temperatura e dall'altezza).
Espressione ancora un po' complicata, ma già si inizia a vedere cosa succede:
il termine esponenziale nella formula per la velocità per t tendente a +∞, tende a 0, annullando tutto l'ambaradan tra parentesi che lo stà a moltiplicare, guardacaso uccidendo anche la dipendenza dalle condizioni iniziali della velocità. D'altronde è noto che per sistemi lineari (e non solo) il punto di equilibrio non dipende dalle condizioni iniziali. L'espressione della velocità finale della particella, una volta fatte le dovute semplificazioni diventa pari a:

v_{z}^{*}=u_{z} + \frac{2 C_{c}R_{p}^{2}g\cdot (\rho_{p}-\rho)}{9\mu}

La velocità u, è da considerarsi con il segno -, se tende a riportare la particella in alto (ascensionale) e con il segno + se la componente lungo z del vento è diretta verso il basso. Per pervenire a questo semplice risultato asintotico non era necessario risolvere tutto quanto, bastava soltanto imporre:

\frac{dv_{z}}{dt}= 0 = g + \frac{6\pi\mu R_{p}}{C_{c}\cdot m_{p}}\cdot(u_{z}-v_{z}^{*})

e ricavare la velocità di equilibrio v* (questo giochetto lo si può fare se la velocità del vento è supposta essere costante).
L'ultimo termine dell'espressione per la velocità di equilibrio in particolare è certo piccolo ma non nullo. Cosa si evince dunque da questo semplice modellino ?
1) Se sufficientemente grande, una particella se ne viene giù a terra senza fare una grinza. Questo implica il fatto che se le correnti in quota sono tali da sostenere una particella pesante (per esempio di polvere di bario), a maggior ragione faranno restare in quota una particella di ghiaccio che ha densità minore. Pertanto non si può discriminare la composizione della scia di condensazione (o chimica) basandosi esclusivamente sui tempi di permanenza in quota delle particelle.
2) Se il tempo di permanenza è sufficientemente lungo, essendo le velocità terminali di equilibrio dipendenti dal raggio delle particelle (ed è ben diversa anche la velocità con la quale le particelle di peso diverso di assesteranno alla velocità terminale) le particelle viaggeranno a velocità differenti, e sostanzialmente si separeranno le une dalle altre. Insomma se ci fosse una scia composta da particolato di bario e goccioline/cristalli di acqua, è lecito supporre che dopo qualche ora, la scia sostanzialmente si separerà in due (o più) scie differenti, una con una maggior concentrazione di cristalli di ghiaccio, l'altra con una maggior concentrazione di particolato di bario. Non mi risulta che tali "doppie scie" siano mai state fotografate o riprese. Tale tecnica si usa nei laboratori (mediante centrifughe) per la separazione di vari composti.
Se in più ci aggiungiamo anche il fatto che le particelle hanno un'elevata igroscopicità formeranno dei fantastici nuclei di aggregazione ed accentueranno ancora di più gli effetti precedentemente descritti, per via dell'incremento della massa e del raggio.
Inutile dire che la granulometria dei presunti particolati non mi sia mai stata fornita, nonostante i sostenitori della teoria sciechimista in questione strepitino tanto sul fatto di aver fatto un'infinità di rilevazioni ed esperimenti, nè è stata fornita alcuna indicazione sui composti presumibilmente coinvolti.
Ora qualcosa circa la densità ρp e sul bario. Non può trattarsi del bario allo stato puro, perchè questo è estremamente reattivo e reagisce molto velocemente a contatto con l'aria. Ora i possibili candidati che ci si aspetterebbe di trovare in una miscela di irrorazione al bario è o l'ossido di bario (BaO) oppure il carbonato (BaCO3). Ora l'ossido di bario è prodotto per decomposizione termica del carbonato di bario a circa 1300 C°. Questa reazione inizia ad avvenire anche a 800 C°. L'ossido di bario però è sensibile alla presenza di H2O, e a contatto con questa prima si trasforma in idrossido di bario, poi continua a idratarsi, fino ad eventualmente diventare barite caustica. Unico particolare che i prodotti intermedi di questa reazione sono molto sensibili alla presenza di CO2 nell'ambiente circostante, ed a contatto con questa si produce ancora una volta il carbonato di bario. Pertanto sembra lecito supporre che è proprio questo il componente che ci si aspetterebbe di trovare nelle scie chimiche. Quindi sostituendo i valori opportuni nelle formule si può ottenere un modello più realistico, cercando di eliminare l'incertezza sul composto a partire da queste semplici considerazioni.
Il fatto però che i composti del bario in esame siano piuttosto igroscopici, accentua ulteriormente l'effetto di accrescimento della massa e della formazione dei nuclei di aggregazione che tenderanno a precipitare perchè appesantiti dall'acqua. Pertanto la lunga permanenza in quota di questi composti non sembra in prima analisi essere compatibile con le loro proprietà chimiche.
Ed a tutto questo si giunge supponendo per assurdo che l'ipotesi sciachimista sia vera.

Ulteriori analisi verranno condotte nella parte II.

Appleman v2.0

Dopo che nell'ultimo articolo è stato fatto vedere come fare delle previsioni "fai da te" sulla formazione delle contrail, con un metodo sostanzialmente "vecchio", è arrivato il momento di spiegare un po' più in dettaglio il processo e perchè no avvalersi di modelli un po' più recenti.

Innanzitutto, come vengono tracciate le curve delle temperature critiche ?
Abbastanza semplice, tutto ciò di cui si ha bisogno è un modello per la pressione di saturazione del vapore acqueo su una superficie di acqua superraffreddata. Avrei potuto anche dire che c'è bisogno di un generatore di campo tachionico trifase, ma in questa sede ne faremo a meno.

Battute a parte vediamo un po' cosa è cosa:
Pressione di saturazione del vapore - è la pressione di equilibrio dinamico della fase gassosa e della fase liquida (in questo caso dell'acqua), o più precisamente è la pressione alla quale, per una fissata temperatura T*, il tasso di evaporazione dell'acqua ed il tasso di condensazione del vapore sono uguali.
Per questi modelli si parla di superraffreddamento perchè il range di temperature per cui tali modelli sono validi sono al di sotto degli 0 C° (273.15 K).
Quindi finora nulla di particolarmente strano. Di modelli per la pressione di saturazione del vapore ce ne sono un'infinità, di seguito riporto ne riporto alcuni:

Goff-Gratch (1946):

\\log{(p_w)} =  -7.90298 (\frac{373.16}{T}-1)+ 5.02808\cdot \log{(\frac{373.16}{T})}- 1.3816\cdot10^{-7} (10^{11.344 (1-\frac{T}{373.16})}-1) + 8.1328\cdot 10^{-3} (10^{-3.49149 (\frac{373.16}{T}-1)}  -1)  + \log(1013.246)

ove pw è la pressione di saturazione del vapore sull'acqua in [hPa], T è la temperatura in [K].

Murphy-Koop (2005):

\log(100\cdot p_w) = 54.842763 - \frac{6763.22}{T} - 4.21\cdot\log(T)+ 0.000367T + \tanh{(0.0415(T - 218.8))}\cdot(53.878 - \frac{1331.22}{T} - 9.44523\cdot\log(T) + 0.014025 T)

ove pw è la pressione di saturazione del vapore sull'acqua in [hPa], T è la temperatura in [K],
log(x) è il logaritmo in base 10 di x, tanh(x) è la tangente iperbolica di x definita come:

\tanh(x)=\frac{\sinh(x)}{\cosh(x)}=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}

con e, base del logaritmo naturale, pari a 2.718281828...
Ora è banale invertire le due formule per ricavare pw in funzione di T: sia nel primo che nel secondo caso infatti si ottengono delle espressioni della forma:

p_{w}(T)=a\cdot e^{f(T)}

Ora, cosa centra questo con le temperature critiche, il contrail factor e tutto il resto ? Presto detto:
Le contrail, sono sostanzialmente nuvole che si formano, quando una miscela di aria fredda e gas di scarico di un motore jet raggiunge la saturazione rispetto all'acqua. In questo modo si formano prima delle goccioline liquide, poi dei cristalli, che derivano dal veloce congelamento di queste goccioline (per questo più avanti si userà la curva di saturazione dell'acqua invece che del ghiaccio, come proposto da qualche autore).
Supponiamo ora di considerare un volume di controllo interno alla scia, come influenza il motore jet le condizioni di quantità di vapore contenuto nel volume di controllo e la sua temperatura ? Abbastanza semplice: basta considerare la quantità

\frac{\Delta r}{\Delta T}

r è quello che si chiama "mixing ratio", ossia la massa di vapore (Kg) su massa d'aria secca. Pertanto quello che esprime la formula è sostanzialmente il comportamento del motore jet in esame: per ogni grado di aumento di temperatura, la variazione del contenuto di vapore nel volume di controllo. Questo fattore è appunto il Contrail factor, funzione delle caratteristiche del motore. Le varie formulazioni del contrail factor sono state derivate praticamente tutte in modo empirico. Ciò non toglie che si possa derivare una relazione più precisa usando un modello più preciso del motore e del processo di combustione.
Il modello proposto da Schumann per il contrail factor è:

\frac{\Delta r}{\Delta T} = CF = \frac{E_{I}C_{p}}{Q(1-\eta)}

ove EI è l'indice di emissione di vapore (dipende dal carburante usato), Cp è il calore specifico dell'aria, Q è la temperatura di combustione, η è l'efficienza del motore, ed è definita come:

\eta = \frac{F \cdot v}{Q\frac{dm}{dt}}

dove, v è la velocità del velivolo (o del motore in esame, dal momento che sarebbe più corretto definire l'efficienza di un motore singolarmente), F è la spinta del motore, Q è la temperatura di combustione, dm/dt è la quantità di combustibile per unità di tempo che viene immessa nel motore e che dovrà essere bruciata.
Come legare tutto questo al discorso sulla pressione di saturazione del vapore ?
La pressione di saturazione del vapore è legata al mixing ratio ed alla pressione atmosferica nel seguente modo:

r = 621.97\cdot\frac{e_s}{p-e_s}

dove, r è il mixing ratio, es è la pressione di saturazione del vapore sull'acqua, e p è la pressione atmosferica in hPa (mb).
Dal momento che non ci si accontenta di trovare soltanto un rapporto incrementale, si fa il passaggio al limite ed il rapporto incrementale (Δr/ΔT) diventa una derivata. Però a questo punto abbiamo tutte le relazioni per poter scrivere un'espressione analitica (per definizione il tutto è ancora uguale a CF):

\frac{dr}{dT}=621.97\cdot\frac{d}{dT}(\frac{e_s(T_c)}{p-e_s(T_c)})=CF

La temperatura che soddisfa la relazione precedente è appunto la temperatura critica per l'umidità relativa del 100% (dal momento che la la curva della pressione di saturazione del vapore sostanzialmente descrive le caratteristiche di un'aria satura e quindi con umidità relativa del 100%, ossia quelle condizioni in o un'abbassamento della temperatura o l'aggiunta di altro vapore provocherebbe la condensazione). Quindi la temperatura critica in qualche modo rappresenta delle condizioni limite per la formazione delle contrail (appunto si impone la tangenza alla curva di saturazione in un punto).
Sviluppando la derivata nella relazione precedente si ottiene:

621.97\cdot\frac{d}{dT}(\frac{e_s(T_c)}{p-e_s(T_c)})=621.97\cdot\frac{p\cdot \frac{de_s(T_c)}{dT}}{(p-e_s(T_c))^{2}}

eguagliando la relazione a CF, si ha:

\frac{\frac{de_s(Tc)}{dT}}{(p-e_s(T_c))^{2}}-\frac{CF}{p\cdot 621.97}=0

si scrive la relazione nella forma f(Tc)=0 per motivi di chiarezza concettuale. Ora data la complicata espressione di es(T) (i due modelli che sono stati mostrati all'inizio della discussione) non è possibile tirare fuori Tc in modo analitico:
nessun problema, per quanto l'espressione in questione possa essere complicata dal momento che le funzioni in questione sono tutte piuttosto regolari (sostanzialmente delle esponenziali) è possibile risolvere f(Tc)=0 in modo numerico-iterativo (ad esempio anche usando il metodo di Newton-Raphson non smorsato, non si hanno problemi di convergenza).
Ora indicando con Tc* la temperatura critica che soddisfa l'equazione precedente (temperatura critica per l'umidità relativa del 100%), la temperatura critica per un'umidità relativa arbitraria si ottiene risolvendo rispetto a Tc:

T_{c}-T_{c}^{*}+\frac{e_{s}(T_{c}^{*})-u\cdot e_{s}(T_{c})}{\frac{de_{s}(T_{c}^{*})}{dT}}=0

il termine al denominatore però è stato ricavato in precedenza, ed è quello che contiene il contrail factor, sostituendo questa relazione nell'ultima equazione si ottiene:

T_{c}-T_{c}^{*}+\frac{e_{s}(T_{c}^{*})-u\cdot e_{s}(T_{c})}{\frac{CF\cdot (p-e_{s}(T_{c}^{*})^{2})}{621.97\cdot p}}=0

Sostituendo ora nell'equazione di sopra i dati numerici relativi alla pressione atmosferica p misurata (hPa) e umidità relativa u (valori nell'intervallo [0,1]) si possono ottenere le relative temperature critiche.
Ecco quindi che si ottendono le curve che normalmente si vedono nei diagrammi di Appleman, per varie umidità relative (la curva più in basso è relativa all'umidità relativa dello 0%, quella più in alto del 100%, quelle intermedie corrispondono a rh incrementati del 10%), le curve sono state tracciate usando il modello Murphy-Koop ed un CF = 0.039:

E' interessante notare che il contrail factor usato da Appleman era ritenuto costante (a differenza della formulazione du Schumann), e pari a 0.0336, valori più aggiornati si aggirano intorno a 0.039 (Boeing 747).
Per esempio ecco cosa succede alla temperatura critica a 0% di umidità relativa usando due contrail factor differenti:

Sostanzialmente quindi alla fine si prevede la formazione di una contrail se la temperatura dell'ambiente è inferiore alla temperatura critica.
Ora cos'ha a che fare tutto ciò con le chemtrails ? Sono convinto che tutti più o meno si siano scontrati con la domanda sciachimista: "Ma prima non erano così, non ce n'erano così tante! Come me lo spieghi questo?". Ecco la risposta a questa domanda: aumento del traffico aereo (osservazione controllata empiricamente) e differente contrail factor (che guardacaso aumenta all'aumentare dell'efficienza, in modo descritto dalla formula di Schumann), maggiore di quello di qualche anno fa a causa di motori più potenti e più efficienti.
Ora avvalendoci ancora una volta di dati reperiti dalla stazione di Pratica di Mare (RM) su questo sito, ci accingiamo a fare delle previsioni per la formazione delle contrail. Stavolta invece di riportare i puntazzi su una jpg, facciamo una cosa un po' più figa:
La temperatura critica è funzione dell'umidità relativa e della pressione. La pressione e l'umidità relativa sono a loro volta funzioni dell'altezza, ed i loro valori sono disponibili dalle misurazioni. Non resta che interpolare (uso un'interpolazione cubica) i dati, mettere le funzioni così ottenute nella formula per il contrail factor ottenendo così una temperatura critica funzione dell'altezza, e confrontare i due grafici con il profilo di temperatura:

ed un ingrandimento della regione di interesse:

I dati sono stati prelevati alle ore 12Z (ore 14.00 Roma) di Domenica 7 Settembre 2008. In rosso è tracciato il profilo delle temperature critiche, in nero, il profilo di temperatura effettivo. Negli intervalli di altezze in cui la temperatura è minore della temperatura critica (curva nera inferiore alla curva rossa), si formeranno le contrail.
Anche oggi niente da fare per presunte contrail a 3000 m di altezza a quanto pare :D.

Riferimenti:
"Calculations of Aircraft Contrail Formation Critical Temperatures" di Mark L. Schrader (1996)

"Contrail Studies and Forecasts in the Subarctic Atmosphere Above FAirbanks,
Alaska" di Martin Stuefer e Gerd Wendler (2004)

"Improved Contrail Forecasting Techniques for the Subarctic Setting of Fairbanks, Alaska"
Gerd Wendler and Martin Stuefer (2002)

"On Conditions for Contrail Formation from Aircraft Exhausts" di U. Schumann (1996)

Un po' di teoria di base

Il fenomeno delle contrail era noto già dal 1921, e fu studiato in modo relativamente approfondito da un tale Appleman nel 1953, che per ovvi motivi tattici cercava di formulare un framework teorico in cui fosse possibile prevedere la formazione delle scie di condensazione.
Elaborò appunto un sistema che si basava sui "Grafici di Appleman" che permettevano di formulare una previsione del fenomeno.
Un Grafico di Appleman è sostanzialmente un oggetto così fatto:

Oppure così:

Insomma nulla di particolarmente strano o complicato finora. Per potersi cimentare nell'epica impresa di previsione della formazione delle contrail abbiamo bisogno di aggiungere un'altro ingrediente: Il profilo atmosferico di umudità e temperatura.
I dati atmosferici, per fortuna, sono liberamente disponibili in rete su questo fantastico sito.
Si sceglie la data, si scegli la regione, si clicca sulla stazione di interesse (Pratica di Mare (RM) per esempio), e si ottiene una bellissima tabella piena di numeretti, con tanto di orario, fatta più o meno così:

Insomma c'è di tutto e di più. Ora basta mettere un po' di dati sul grafico:
In particolare interessano le colonne: (PRES (Pressione), HGHT (Altezza), TEMP (Temperatura) RELH (Umidità relativa)).
Ora non resta che mettere in pratica i seguenti due banalissimi passaggi:
1) Dalla colonna TEMP e dalla colonna PRES si ottengono le informazioni circa il profilo di temperatura dell'atmosfera; ora è sufficiente riportare questi punti sul grafico nel punto di coordinate (TEMP,PRESS) (pallino rosso).
2) Dalla colonna RELH e dalla colonna PRES si ottengono le informazioni circa l'umidità relativa. Qui il gioco è lo stesso di prima solo che è un poco più complicato (ma non di molto):
Per segnare un punto relativo ad una pressione ed ad una umidità relativa, basta metterlo al punto di intersezione della curva della temepratura critica corrispondente all'umidità relativa che si sta considerando con il relativo valore della pressione.
Ecco cosa succede oggi su Pratica di Mare (RM):

Come si fa a stabilire se ad una certa quota si formeranno scie di condensazione ? Nulla di più semplice: basta che per quella particolare quota il pallino rosso sia ben a sinistra del pallino blu.
In pratica vi sono le condizioni per la formazione della scie nell'intervallo tra gli 11176 ed i 16540 metri di quota.

In pratica è quello che vedo da casa mia: neppure una scia persistente, con buona pace dei complottisti.

P.s.
Come fattomi gentilmente notare da brain_use, il modello applicato da Appleman è ormai datato e passato all'obsolescenza.

Delle Scie chimiche (Introduzione)

Navigando quà e là per la rete mi è capitato di imbattermi in persone che sostenevano a gran voce la cosidetta teoria delle scie chimiche. La principale "ipotesi" di questa teoria del complotto consiste più o meno nell'affermazione del fatto che le scie di condensazione persisteni lasciate dagli aerei che si vedono nel cielo, siano in realtà appunto delle scie "chimiche" che contengono i materiali ed i particolati più disparati:
Ormoni;
Nanofibre;
Particolati nanometrici;
Bario;
Alluminio;
Silicio;
Materiale biologico di varia origine e provenienza (magari non terrestre);
Insomma, di tutto e di più.
Quanto alle ipotesi sullo scopo presunto di queste scie, ci sono ancora più ipotesi che sui loro contenuti, si va dal controllo del clima, al controllo della mente, alle armi a microonde, allo sterminio della popolazione (circa 4 miliardi di persone nel progetto attuale)
L'ipotesi più fantasiosa è a mio avviso quella sulla sperimentazione di una "Bomba Gay" (Eghm, voi del NWO ne sapevate qualcosa ? ;) ), e qui non aggiungo altro. Ecco il link. Ricorre molto spesso anche il progetto HAARP (ecco il link a sito ufficiale, liberamente accessibile non solo dalla rete), insomma non ci sono credo di tutti i complottisti che concordino sui contenuti delle scie e sugli scopi.
Di seguito linko un po' di risorse:

Discussione sul forum di Focus - link
Discussione sul forum di Giovani.it - link
Blog di Hanmar - link
Blog di Straker (per par condicio) - link
Contrail Science - lo trovate tra i "Link Interessanti" nella colonnina di destra in questo blog.
Le risorse della rete sull'argomento sono davvero numerosissime, basta cercare su Google.

Le domande (e le asserzioni) fondamentali con cui molti sciachimisti sono soliti sostenere le loro teorie sono più o meno le seguenti:
1) Perchè le scie durano così tanto ?
1.a) Prima non duravano tanto a lungo e non si espandevano tanto, Perchè ora si espandono tanto e durano molto di più di quanto non durassero qualche anno fa?
2) Perchè sono irregolari ("intermittenti") ?
3) Perchè quando gli aerei girano e fanno manovre sulla mia testa fanno la scia strana e quando vanno via fanno la scia "normale" ?
4) Perchè sulla mia testa vedo reticoli di scie che si incrociano? Non sarà perchè stanno irrorando ?
5) Perchè gli aerei volano così bassi ?

Inutile dire che le foto dei presunti congegni per lo spruzzamento sono puntualmente smentite.