giovedì 11 settembre 2008

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.

2 commenti:

Cuorenero ha detto...

Ho ricontrollato le formule a plottato i grafici contronfanto con i risultati sperimentali ... molto vicini. Posterò i grfici nel prossimo articolo con tanto di simulazione.

blog ha detto...

Thank you
filmini izle