giovedì 16 ottobre 2008

Simulazione (parte II)



Dopo le ultime simulazioni di un semplice modellino lineare della particella, eccoci qui alle prese con l'equazione completa della dinamica della particella. Si suppongono numeri di Reynolds relativamente bassi e numeri di Knudsen prossimi allo 0. Ricordo che il numero di Reynolds non è altro che il rapporto tra le forze di inerzia e le forze viscose. Ovviamente la scelta di usare approssimazioni per bassi numeri di Reynolds è perfettamente giustificata, dal momento che la massa delle particelle prese in esame è estremamente ridotta. I numeri di Knudsen (il rapporto tra il ''mean free path'' ed il raggio equivalente) esprimono sostanzialmente quanto sia valida l'approssimazione continua del mezzo rispetto alle dimensioni della particella, che nel nostro caso nonostante sia (piccola, piccolissima, minuscola, microscopica, nanodimensionata) è comunque significativamente più grande delle molecole che compongono l'aria. Lo slip correction factor quindi dovrebbe essere circa pari a 1, ma per eccesso di zelo e per gusto sadico del fare del male alla CPU lo includiamo e lo rendiamo perfino dipendente dall'altezza a cui si trova la perticella. Tutte le particelle che vengono considerate nelle successive simulazioni inoltre sono non-Strakeriane. Ricordo a chi non lo sapesse che una particella Strakeriana è una particella di cui sono ignote (e non possono essere misurate in modo alcuno con nessuno strumento) le dimensioni, che non ubbidisce ad alcuna legge fisica nota o ignota, la cui dinamica non è esprimibile con nessun formalismo presente, passato o futuro e la cui composizione varia a piacimento.

Ecco dunque la brutta bestia di equazione che ci toccherà simulare (in questi momenti ringrazio Turing e Von Neumann per aver inventato i PC):

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


+\frac{V_{p}}{2}\rho(z)\left(\frac{d u_{z}}{dt}-\frac{d v_{z}}{dt}\right)
+6R_{eq}^{2}\sqrt{\pi\mu(z) \rho(z)}\int\limits_{0}^{t}
\frac{\displaystyle \left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=\tau}}{\sqrt{t-\tau}}d\tau-m_{p}g

I parametri presenti nell'equazione sono stati già accuratamente descritti nell'articolo precedente, pertanto non si ritiene necessario ripetere tutta la pappardella ancora una volta.
Quanto ai termini che compaiono nell'equazione:

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

non è altro che la forza di Stokes che agisce sulla particella in moto nel fluido.

\frac{9}{4}\pi \rho(z) \left( \frac{R_{eq}\chi}{C_{c}(z)} \right)^2 (u_{z}-v_{z})|u_{z}-v_{z}|

E' la correzione di Ossen all'equazione di Stokes, che serve per tener conto di quello che succede anche lontano dalla particella.

V_{p}\rho(z)\frac{d u_{z}}{dt}

E' un termine che viene introdotto dal gradiente di pressione che si viene a creare nel fluido a causa del moto della particella, in quanto il fluido in questione, in prossimità delle pareti della particella tende ad accelerare.

\frac{V_{p}}{2}\rho(z)\left(\frac{d u_{z}}{dt}-\frac{d v_{z}}{dt}\right)

Esprime l'effetto del flusso potenziale intorno alla particella in esame, e serve a tener conto anche dell'inerzia del fluido che circonda la particella.

+6R_{eq}^{2}\sqrt{\pi\mu(z) \rho(z)}\int\limits_{0}^{t}
\frac{\displaystyle \left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=\tau}}{\sqrt{t-\tau}}d\tau

Rappresenta la forza di Basset e tiene conto della storia dell'accelerazione della particella, accoppiandola alla densità del mezzo. In pratica è il termine che tiene conto della vorticità del fluido che si viene a creare a causa del moto della particella.

m_{p}g

E' la cara vecchia forza di gravità (come detto all'inizio, vengono considerate soltanto particelle non-Strakeriane).

Una volta che abbiamo snocciolato tutti i termini e tutte le costanti, non possiamo ancora fare alcun chè per la risoluzione dell'equazione, in quanto c'è una singolarità che si annida nell'integrale. Mi spiego meglio: il termine con la radice al denominatore dell'integrale è discontinuo in τ =t. Dal momento che quel che c'è sopra è praticamente una costante, eccoci che nel fare i conti abbiamo un termine infinito come integrando (costante/0). Nessuna paura però, con un'ardito passaggio matematico che consiste nell'aggiungere e togliere la stessa cosa a una certa quantità uccidiamo la singolarità brutta, e ne creiamo una che numericamente è molto più tranquilla (0/0):

\int\limits_{0}^{t}\frac{\displaystyle \left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=\tau}}{\sqrt{t-\tau}}d\tau
\underbrace{-\int\limits_{0}^{t}\frac{\displaystyle \left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=t}}{\sqrt{t-\tau}}d\tau
+\int\limits_{0}^{t}\frac{\displaystyle \left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=t}}{\sqrt{t-\tau}}d\tau}_{=0}

Ovviamente non abbiamo fatto altro che aggiungere 0 all'integrale originario (e chi lo ha detto che le addizioni di zeri sono
futili ???).
Ora però possiamo (in virtù della linearità dell'operatore integrale) fare dei primi due integrali uno solo:

\int\limits_{0}^{t} \left[\left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=\tau}-\left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)
\right|_{t=t}\right] \frac{\displaystyle 1}{\sqrt{t-\tau}}d\tau
+\int\limits_{0}^{t}\frac{\displaystyle \left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=t}}{\sqrt{t-\tau}}d\tau

La cosa carina in questo passaggio è che il secondo integrale può essere tranquillamente risolto, perchè le funzioni al numeratore non dipendono più da τ e possono essere portate fuori dall'integrale:

\int\limits_{0}^{t}\frac{\displaystyle \left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=t}}{\sqrt{t-\tau}}d\tau =
\left(\frac{d}{dt}(u_{z}-v_{z})\right)\cdot\int\limits_{0}^{t}\frac{\displaystyle 1}{\sqrt{t-\tau}}d\tau=
\left(\frac{du_{z}}{dt}-\frac{dv_{z}}{dt}\right)\cdot 2\sqrt{t}

Sostituendo il tutto nell'equazione otteniamo:

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


+6R_{eq}^{2}\sqrt{\pi\mu(z) \rho(z)}\left( \int\limits_{0}^{t}\left[\left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=\tau}-\left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=t}\right] \frac{\displaystyle 1}{\sqrt{t-\tau}}d\tau\right) +


+6R_{eq}^{2}\sqrt{\pi\mu(z) \rho(z)}\left(\left(\frac{du_{z}}{dt}-\frac{dv_{z}}{dt}\right)\cdot 2\sqrt{t} \right)-m_{p}g



Ora, dividiamo tutto per mp tanto per ricondurci a qualcosa di più o meno leggibile.
In linea di massima abbiamo una cosa del tipo:

\frac{dv_{z}}{dt}=
l(z)(u_{z}-v_{z})
+g(z)(u_{z}-v_{z})|u_{z}-v_{z}|
+(3 h(z)+2\sqrt{t}f(z))\frac{d u_{z}}{dt}


-(h(z)+2\sqrt{t}f(z)) \frac{d v_{z}}{dt}
+f(z) \int\limits_{0}^{t}\Psi\left(\tau,t,v_{z}(\tau),\frac{dv_{z}}{dt}(\tau)\right) d\tau -g

con:

\begin{array}{l}
f(z)=6\frac{R_{eq}^{2}\sqrt{\pi\mu(z) \rho(z)}}{m_{p}}\\
\\
\Psi\left(\tau,t,v_{z}(\tau),\frac{dv_{z}}{dt}(\tau)\right)=
\left[\left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=\tau}-\left.\left(\frac{d}{dt}(u_{z}-v_{z})\right)\right|_{t=t}\right]
\frac{\displaystyle 1}{\sqrt{t-\tau}}\\
\\
h(z)=\frac{V_{p}}{2 m_{p}}\rho(z)\\
\\
l(z)=\frac{6\pi \mu(z) R_{eq}\chi}{C_{c}(z)m_{p}}\\
\\
g(z)=\frac{9}{4 m_{p}}\pi \rho(z) \left( \frac{R_{eq}\chi}{C_{c}(z)} \right)^2
\end{array}

Siccome però non abbiamo intenzione di stare a combattere con equazioni integrali-algebrico-differenziali, con un ardito passaggio matematico otteniamo:

(1+h(z)+2\sqrt{t}f(z))\frac{dv_{z}}{dt}=
l(z)(u_{z}-v_{z})
+g(z)(u_{z}-v_{z})|u_{z}-v_{z}|
+(3 h(z)+2\sqrt{t}f(z))\frac{d u_{z}}{dt}


+f(z) \int\limits_{0}^{t}\Psi\left(\tau,t,v_{z}(\tau),\frac{dv_{z}}{dt}(\tau)\right) d\tau -g

che con un passaggio matematico altrettanto ardito diventa:

\frac{dv_{z}}{dt}=
M(z,t)(u_{z}-v_{z})
+G(z,t)(u_{z}-v_{z})|u_{z}-v_{z}|
+H(z,t)\frac{d u_{z}}{dt}


+F(z,t)\int\limits_{0}^{t}\Psi\left(\tau,t,v_{z}(\tau),\frac{dv_{z}}{dt}(\tau)\right) d\tau -P(z,t)

Dove:

\begin{array}{l}
M(z,t)=\frac{l(z)}{1+h(z)+2\sqrt{t}f(z)}\\
\\
G(z,t)=\frac{g(z)}{1+h(z)+2\sqrt{t}f(z)}\\
\\
H(z,t)=\frac{3 h(z)+2\sqrt{t}f(z)}{1+h(z)+2\sqrt{t}f(z)}\\
\\
F(z,t)=\frac{f(z)}{1+h(z)+2\sqrt{t}f(z)}\\
\\
P(z,t)=\frac{g}{1+h(z)+2\sqrt{t}f(z)}
\end{array}

Ora però è il momento di iniziare a pensare alla risoluzione del problema vero e proprio. In questo caso si procederà con una tecnica di risoluzione approssimata analitica, invece che puramente numerica dal momento che per gestire la memoria introdotta dal termine integrale nell'equazione sarebbe stato necessario scrivere codice da zero ed implementare a mano metodi numerici ad hoc per l'equazione. Quindi almeno in questa sere per la simulazione si preferisce illustrare un approccio che non sia la classica variazione su tema dei vari Runge-Kutta e che non si veda esattamente tutti i giorni. L'approccio in questione è appunto il "Variational Iteration Method". Vediamo però un po' come funziona: costruiamo a partire dall'ultima relazione che abbiamo ricavato una formulazione variazionale del problema, ed un funzionale di correzione. In sostanza a questo punto lo scopo è cercare di ricavare una formula iterativa a partire dall'equazione originaria, che consenta per passi successivi di arrivare alla soluzione. Ad ogni passo la soluzione verrà migliorata, e man mano che si procede, ci si muove verso un certo punto di minimo: la soluzione esatta. Si spera inoltre che il metodo in questione sia convergente e stazionario. Molto più semplice da scrivere che da spiegare. Ora riscriviamo l'equazione in modo da eguagliarla tutta a 0.

\frac{dv_{z}}{dt}-
M(z,t)(u_{z}-v_{z})
-G(z,t)(u_{z}-v_{z})|u_{z}-v_{z}|
-H(z,t)\frac{d u_{z}}{dt}


-F(z,t)\int\limits_{0}^{t}\Psi\left(\tau,t,v_{z}(\tau),\frac{dv_{z}}{dt}(\tau)\right) d\tau +P(z,t) = 0

Per maggiore chiarezza per adesso scriviamo semplicemente

v' - V(t,z,v)-F(z,t)\int\limits_{0}^{t}\Psi(\tau,t,v,v')d\tau = 0

Ora però dobbiamo risolvere per la posizione lungo l'asse z, per cui si sostituisce v con z':

z'' - V(t,z,z')-F(z,t)\int\limits_{0}^{t}\Psi(\tau,t,z',z'')d\tau = 0

A questo punto però dal momento che non è strettamente necessario ricondurci ad un sistema del primo ordine per ora lasciamo le cose come stanno senza aggiungere equazioni ausiliarie.Una volta scritta in forma così generale ed elementare l'equazione, scriviamo adesso il funzionale di correzione che adopereremo per trovare la soluzione approssimata:

\xi_{n+1}(t)=\xi_{n}(t)+\int\limits_{0}^{t}\lambda(s)\left[\xi_{n}''(s) - V(s,\xi_{n},\xi_{n}')-
F(\xi_{n},s)\int\limits_{0}^{s}\Psi(\tau,s,\tilde{\xi}_{n}',\tilde{\xi}_{n}'')d\tau \right]ds

Ove i λ(s) sono i moltiplicatori generalizzati di Lagrange (per adesso non noti). Chi ha fatto degli esami (o letto qualcosa) di controllo ottimo o ottimizzazione e programmazione dinamica sà bene di cosa sto parlando. Questi verranno scelti in modo ''ottimo'' rispetto alle variazioni di ξn, cioè: δξn.
Ora se scriviamo la variazione del funzionale di correzione si ottiene banalmente che:

\delta\xi_{n+1}(t)=\delta\xi_{n}(t)+\delta\int\limits_{0}^{t}\lambda(s)\left[\xi_{n}''(s) - V(s,\xi_{n},\xi_{n}')-
F(\xi_{n},s)\int\limits_{0}^{s}\Psi(\tau,s,\tilde{\xi}_{n}',\tilde{\xi}_{n}'')d\tau \right]ds

Ora dal momento che stiamo considerando variazioni limitate, il termine che contiene Ψ(.) sparisce e della funzione di cui sopra, e resta soltanto:

\delta\xi_{n+1}(t)=\delta\xi_{n}(t)+\delta\int\limits_{0}^{t}\lambda(s)\left[\xi_{n}''(s) - V(s,\xi_{n},\xi_{n}')\right]ds

Ora spezziamo l'integrale in due pezzi:

\delta\xi_{n+1}(t)=\delta\xi_{n}(t)+
\delta\int\limits_{0}^{t}\lambda(s)\xi_{n}''(s)ds-
\delta\int\limits_{0}^{t}\lambda(s)V(s,\xi_{n},\xi_{n}')ds

Ora sviluppiamo con l'integrazione per parti il primo termine integrale per due volte:

\delta\xi_{n+1}(t)=\delta\xi_{n}(t)+
\left.\lambda(s)\right|_{s=t}\delta\xi_{n}'(s)-\left.\lambda'(s)\right|_{s=t}\delta\xi_{n}(s)+
\int\limits_{0}^{t}\lambda''(s)\delta\xi_{n}(s)ds
-\delta\int\limits_{0}^{t}\lambda(s)V(s,\xi_{n},\xi_{n}')ds

riscrivendo un po' i primi termini mettendo in evidenza δξn otteniamo:

\delta\xi_{n+1}(t)=\left(1-\left.\lambda'(s)\right|_{s=t}\right)\delta\xi_{n}(t)+
\left.\lambda(s)\right|_{s=t}\delta\xi_{n}'(s)+
\int\limits_{0}^{t}\lambda''(s)\delta\xi_{n}(s)ds
-\delta\int\limits_{0}^{t}\lambda(s)V(s,\xi_{n},\xi_{n}')ds

Ora procediamo a definire le condizioni di stazionarietà della variazione del funzionale di correzione (questo è dovuto al fatto che rispetto ai λ stiamo cercando la variazioni minima, quindi imponiamo la variazione nulla). La cosa consiste nel considerare trascurabile il termine integrale con la funzione V, e quindi imporre come terna di condizioni semplicemente:

\left\{
\begin{array}{l}
\left(1-\left.\lambda'(s)\right|_{s=t}\right) = 0\\
\\
\left.\lambda(s)\right|_{s=t} = 0\\
\\
\lambda''(s)=0
\end{array}
\right.

che risolte, danno luogo al moltiplicatore di Lagrange:

\lambda(s)=(s-t)

ed al funzionale di correzione:

\xi_{n+1}(t)=\xi_{n}(t)+\int\limits_{0}^{t}(s-t)\left[\xi_{n}''(s) - V(s,\xi_{n},\xi_{n}')-
F(\xi_{n},s)\int\limits_{0}^{s}\Psi(\tau,s,\xi_{n}',\xi_{n}'')d\tau \right]ds

Teoricamente si potrebbe andare a smucinare nel termine integrale con V(.), ma ci si rende velocemente conto che non ha molto senso perchè non ha termini lineari in ξn ed i coefficienti sono per lo più funzioni implicite, pertanto in ogni caso non ne uscirebbero fuori dei termini che vadano ad influire sull'impostazione delle condizioni.

Ogni metodo iterativo che si rispetti però ha bisogno di qualcosa da cui partire. Il nostro metodo non fa certo eccezione. Guarda caso abbiamo già a disposizione un po' di dati per poter formulare un'approssimazione della scelta iniziale della soluzione approssimata. Se ci ricordiamo delle vecchie simulazioni con il modello lineare, e della soluzione esplicita dell'equazione differenziale associata, ci rendiamo presto conto che il gioco è praticamente fatto: la velocità è stabile e tende a quella del vento molto velocemente, la posizione è praticamente definita da una legge lineare.
Ecco quindi che praticamente abbiamo che:

\xi_{0}(t) = z_{0}+v_{z}(0)t

Arrivati a questo punto però volendo si può usare un piccolo trucco, dal momento che è immediatamente evidente che la funzione

\xi_{0}(t) = z_{0}+u_{z}t

nel caso della velocità del vento costante è quella che fà si che si raggiunga immediatamente l'equilibrio e che le iterazioni successive non modifichino in alcun modo il risultato.

Dopo numerosi test cercando di preservare maggiormente l'equazione però si giunge alla conclusione che un'equazione così complicata non possa essere gestita in modo conveniente da un metodo analitico, sopratutto a causa delle funzioni implicite presenti.
Per la soluzione numerica si procederà nell'articolo successivo con il metodo della Collocazione.

giovedì 18 settembre 2008

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.

giovedì 11 settembre 2008

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.

mercoledì 10 settembre 2008

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)

lunedì 8 settembre 2008

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.