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.
