Brent Corrigan Free Stream

Simulation (Part I)

 
 Here we are again after a short pause (due examination of Control Systems (mental of the population) with multiple variables). In '
 

last article on modeling the behavior of a single particle at the end of everything we had worked out a complete model of the behavior of a particle immersed in an incompressible fluid with a laminar flow.
Now we have to put together the results obtained and presented in various articles, and insert them all in the original model to make use potersene during the simulation.

 recall for convenience that the model in question is defined by the following: 
 
 

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, C
c
Slip Correction Factor, R

 eq 
  il raggio equivalente della particella, m 
 p 
  massa della particella (considerando il galleggiamento), V 

p
volume della particella, u
z

  velocità del vento lungo l'asse z, v 
 z 
  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 particle is the sum of its results (ρ + ρ 
 to 
 
 
 v = ρ (z)). 

So far no problem, neither numerical nor analytical (explicit formulas for the roots of an equation of third degree even though there are some 'complicated). The only care for those who wish to use the explicit formulas: first, be sure to use the solution (the positive), according to all attention to the order you put the number of transactions (you may need to turn over a bit 'formula) trying to minimize the error propagation. In this case, just a method of Newton-Raphson no cushion to have a good and fast convergence.

Ora ultimo ma non per importanza il C

 c 
 . 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 C
c
sfido chiunque a risolversela a mano ed a tirare fuori un'espressione esplicita per C
c
. 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;