Processi di advezione-diffusione, Guide, Progetti e Ricerche di Idraulica

Scarica Processi di advezione-diffusione e più Guide, Progetti e Ricerche in PDF di Idraulica solo su Docsity! DIPARTIMENTO DI INGEGNERIA E ARCHITETTURA CORSO DI LAUREA MAGISTRALE DI INGENGERIA CIVILE CURRICULUM: PROTEZIONE DAL RISCHIO SISMICO, IDRAULICO E AMBIENTALE RELAZIONE TECNICA CORSO DI IDRAULICA AMBIENTALE METODI NUMERICI DI SIMULAZIONE DEI PROCESSI DI TRASPORTO ADVETTIVO – DIFFUSIVO: MATLAB Docente: Prof. Vincenzo Armenio Studenti: Giovanni Coan Sebastiano Troian Data di prima stesura: 18/05/2021 1 Giovanni Coan Sebastiano Troian SOMMARIO 1 Abstract ................................................................................................................................... 2 2 Introduzione alla teoria sulla dispersione di inquinanti ............................................................. 4 3 Impostazione metodo di risoluzione numerica .......................................................................... 7 3.1 Processo puramente diffusivo con rilascio istantaneo e puntuale ....................................... 7 3.2 Processo puramente diffusivo con rilascio istantaneo e distribuito su una lunghezza ....... 13 3.3 Processo puramente diffusivo con rilascio istantaneo e puntuale tra due pareti ................ 16 3.4 Processo puramente diffusivo con rilascio istantaneo e puntuale con decadimento .......... 19 3.5 Problema di diffusione anisotropa in dominio 3D ............................................................ 21 3.6 Processo di advezione – diffusione con rilascio istantaneo e puntuale ............................. 24 3.7 Processo di advezione – diffusione con rilascio istantaneo e puntuale con decadimento .. 26 3.8 Processo advettivo – diffusivo: Problema bidimensionale ............................................... 28 4 Giovanni Coan Sebastiano Troian 2 INTRODUZIONE ALLA TEORIA SULLA DISPERSIONE DI INQUINANTI Lo studio dell’emissione di sostanze all’interno del Atmospheric Boundary Layer (ABL)1 ha come scopo l’individuazione dell’andamento della concentrazione nello spazio e nel tempo. In questo caso ci limiteremo a studiare il trasporto della sostanza, senza considerare senza considerare gli effetti della stratificazione della colonna di fluido (sia riferita all’atmosfera che ad un bacino idrico). L’analisi del sistema può essere effettuata applicando due approcci diversi:  Approccio Euleriano in cui la dinamica del trasporto delle sostanze è descritta prendendo come riferimento un volume di controllo fisso nello spazio, ma aperto, in quanto consente lo scambio sia di massa che di energia con il resto del fluido;  Approccio Lagrangiano dove si prende a riferimento un volume materiale in moto con il fluido, che costituisce un sistema chiuso in quanto non scambia massa, ma solo energia con il resto del fluido. Nello studio dei processi e nell’applicazione degli algoritmi numerici per ottenere i risultati si adotta l’approccio Euleriano del volume di controllo. Quest’ultimo permette la scrittura delle equazioni di bilancio, riferite appunto al volume di controllo. Per la definizione del modello numerico faremo, però, riferimento alle leggi di conservazione scritte in forma indefinita, cioè valide in qualsiasi punto del dominio di moto, e ottenute da quelle riferite al volume di controllo attraverso un’analisi differenziale. Di seguito riportiamo l’espressione finale delle leggi di conservazione.  Equazione di continuità: Nel caso di fluido comprimibile abbiamo: 𝜕𝜌 𝜕𝑡 + ∇ ∙ (𝜌𝐯) = 0 Nel caso di incomprimibilità del fluido diventa: div(𝒗) = 0  Equazione di conservazione della quantità di moto: 𝜌 𝑑𝐯 𝑑𝑡 = 𝜌𝐠 + ∇ ∙ 𝜎𝑖𝑗 1 “L’ Atmospheric Boundary Layer è quella porzione di troposfera direttamente influenzata dalla superficie terrestre, che risponde all’immissione di energia da essa proveniente con scale temporali dell’ordine dell’ora.” (Stull, 1989). 5 Giovanni Coan Sebastiano Troian Il tensore degli sforzi 𝜎𝑖𝑗 può essere separato, grazie alle equazioni costitutive, in una parte isotropa di modulo pari alla pressione idrostatica locale p, e una parte deviatorica che tiene conto degli sforzi tangenziali viscosi. 𝜎𝑖𝑗 = −𝑝𝑖𝑖 + 𝜏𝑖𝑗 Prendendo in considerazioni fluidi newtoniani, incomprimibili e con variazioni locali di temperatura trascurabili, quindi con viscosità dinamica e cinematica costanti, il deviatore degli sforzi può essere espresso come: 𝜏𝑖𝑗 = 2𝜇𝜀𝑖𝑗 Quindi ottengo: 𝜌 𝑑𝐯 𝑑𝑡 = −∇𝑝 + 𝜇∇2𝐯 + 𝜌𝐠  Equazione di conservazione della massa di una concentrazione disciolta in un corpo fluido: 𝜕𝐶 𝜕𝑡 + 𝜕𝐶 𝜕𝑥𝑖 𝑣𝑖 = 𝑘 𝜕2𝐶 𝜕𝑥𝑖𝜕𝑥𝑖 + 𝑆 Concentriamo la nostra attenzione su quest’ultima equazione. Si possono distinguere due tipi di trasporto della massa di un fluido, quello advettivo e quello diffusivo. 𝜕𝐶 𝜕𝑡 + 𝜕𝐶 𝜕𝑥𝑖 𝑣𝑖 = 𝑘 𝜕2𝐶 𝜕𝑥𝑖𝜕𝑥𝑖 + 𝑆 Trasporto advettivo: è il flusso della sostanza che viene trasportata in un campo di moto con velocità v. Il flusso advettivo di una sostanza attraverso una superficie dA è proporzionale alla velocità del fluido in direzione normale a dA: 𝑭𝑎 = 𝐶𝐯𝑛𝑑𝐴 Trasporto diffusivo: è il termine di trasporto dovuto alle interazioni molecolari. Questo termine è proporzionale al gradiente di concentrazione, attraverso il coefficiente di diffusione 𝑘, e si attua nel senso delle concentrazioni decrescenti (down gradient). La legge che rappresenta questo fenomeno è la legge di Fick: 𝑭𝑑 = −𝜌𝑘𝛁𝐶 6 Giovanni Coan Sebastiano Troian Sempre con riferimento alla stressa equazione, consideriamo un moto di tipo turbolento in cui il flusso è caratterizzato da elementi di casualità nella distribuzione spaziale e temporale delle proprietà macroscopiche. Il comportamento turbolento è puramente caotico e può essere trattato in termini statistici. A questo fenomeno si cerca di trovare una soluzione “media” attraverso l’approccio delle equazioni mediate di Reynolds2. La proprietà statistica della turbolenza porta alla riscrittura delle equazioni di Navier – Stokes in termini di moto medio considerando le componenti di velocità e pressione come somma di una parte stazionaria media dipendente solo dallo spazio e una fluttuazione dipendente sia dallo spazio che dal tempo ma con media nulla. L’equazione che si ottiene è la seguente: 𝜕𝑈?̅? 𝑈?̅? 𝜕𝑥𝑗 = − 1 𝜌 𝜕?̅? 𝜕𝑥𝑖 + 1 𝜌 𝜕 𝜕𝑥𝑗 ( 𝜇 ( 𝜕𝑈?̅? 𝜕𝑥𝑗 + 𝜕𝑈?̅? 𝜕𝑥𝑖 ) − 𝜌𝑢𝑖𝑢𝑗̅̅ ̅̅ ̅) I termini aggiuntivi 𝜌𝑢𝑖𝑢𝑗̅̅ ̅̅ ̅ sono detti sforzi di Reynolds, che sono incogniti. In questo modo otteniamo un sistema di cinque equazioni in 14 incognite, non determinato. Per risolvere il problema della chiusura, Boussinesq introdusse il concetto di viscosità turbolenta: 𝜏𝑖𝑗 = 2𝜌 (𝜈𝑇 + 𝜈)𝑆𝑖𝑗 ̅̅̅̅ = 𝜌 (𝜈 + 𝜈𝑇) ( 𝜕𝑈?̅? 𝜕𝑥𝑗 + 𝜕𝑈?̅? 𝜕𝑥𝑖 ) Allo stesso modo posso riscrivere il termine di trasporto diffusivo dell’equazione di conservazione della massa di una concentrazione introducendo il termine di diffusività turbolenta 𝑘𝑇. Per concludere, l’equazione di conservazione della massa di una concentrazione disciolta in un corpo fluido, che sta alla base dei processi da noi studiati, può essere riscritta come: 𝜕𝐶 𝜕𝑡 + 𝜕𝐶 𝜕𝑥𝑖 𝑣𝑖 = (𝑘 + 𝑘𝑇) 𝜕2𝐶 𝜕𝑥𝑖𝜕𝑥𝑖 + 𝑆 Questo nell’ipotesi semplificativa che (𝑘 + 𝑘𝑇) sia costante nello spazio3. Infine, dato che 𝑘𝑇 ≫ 𝑘, possiamo approssimare il coefficiente di diffusività a quello di diffusività turbolenta: 𝐷~𝑘𝑇 𝜕𝐶 𝜕𝑡 + 𝜕𝐶 𝜕𝑥𝑖 𝑣𝑖 = 𝐷 𝜕2𝐶 𝜕𝑥𝑖𝜕𝑥𝑖 + 𝑆 2 RANS: Reynolds Averaged Navier – Stockes 3 In realtà abbiamo che 𝑘𝑇(𝐱; t). 9 Giovanni Coan Sebastiano Troian Figura 3.1: Condizioni iniziali processo puramente diffusivo con rilascio puntuale e istantaneo Partendo da queste condizioni iniziali, viene applicato alla soluzione discreta un loop temporale e spaziale, con le relative condizioni al contorno del dominio. Quest’ultime sono individuate come principio di parità dei flussi alle estremità del dominio spaziale. Di seguito viene riportata la gaussiana che rappresenta l’andamento della soluzione nel tempo e nello spazio. Si può notare come, l’andamento della concentrazione sia simmetrico rispetto al punto di rilascio. Trattando il caso di immissione istantanea la concentrazione al punto di rilascio decadrà gradualmente nel tempo, mentre, negli altri punti del dominio, dato che la massa rilasciata si conserva (significa che area sottesa alle curve gaussiane è costante nel tempo), la concentrazione prima aumenta per poi decrescere. soluzione numerica condizioni al contorno di parità dei flussi 10 Giovanni Coan Sebastiano Troian Figura 3.2: Profilo di concentrazione al variare del tempo nel caso di processo puramente diffusivo con rilascio istantaneo e puntuale Per una migliore visualizzazione dell’andamento gaussiano della soluzione riportiamo di seguito i valori della concentrazione nel dominio stampati in un minor numero di istanti temporali. Figura 3.3: Profilo di concentrazione al variare del tempo nel caso di processo puramente diffusivo con rilascio istantaneo e puntuale Per valutare la bontà e il livello di approssimazione del metodo di simulazione numerica confronto i risultati ottenuti con quest’ultimo, con i risultati ottenuti applicando la soluzione analitica. 11 Giovanni Coan Sebastiano Troian  Confronto andamento della soluzione analitica (arancione) e numerica (blu) al punto di rilascio al procedere del tempo: Figura 3.4: Confronto andamento della soluzione analitica e numerica al punto di rilascio in funzione del tempo Notiamo come l’andamento nel tempio sia decrescente e, rispetto alla soluzione analitica, i valori ottenuti inizino a discostarsi a valori del tempo pari a 1.5 × 104𝑠, cioè all’incirca quattro ore.  Confronto andamento della soluzione analitica e numerica nello spazio al tempo generico 𝑇 = 14000𝑠 ≈ 4ℎ: Figura 3.5: Confronto soluzione analitica e soluzione numerica al tempo T=4h nel caso di processo puramente diffusivo con rilascio istantaneo e puntuale nel dominio spaziale Discretizzazione spaziale C o n ce n tr az io n e Tempo C o n ce n tr az io n e 14 Giovanni Coan Sebastiano Troian Figura 3.8: Andamento temporale della concentrazione nel caso di processo puramente diffusivo con rilascio istantaneo e distribuito Confrontiamo i risultati ottenuti con il metodo numerico e quelli ottenuti con la soluzione analitica in corrispondenza della mezzeria del dominio, dimostrando che la concentrazione nel tempo tende a dimezzarsi. La soluzione analitica del problema diffusivo, con rilascio istantaneo e distribuito su una lunghezza è: 𝐶(𝑥; 𝑡) = 𝐶0 2 [1 + erf ( 𝑥 √4𝐷𝑡 ) ] Condizione al contorno di parità dei flussi e valore iniziale di C all’estremo superiore del dominio 15 Giovanni Coan Sebastiano Troian Figura 3.9: Confronto soluzione analitica e numerica nella sezione centrale del dominio nel caso di processo puramente diffusivo e rilascio istantaneo e distribuito Si può notare come, sia la soluzione analitica (blu) che numerica (arancione), procedano da un valore iniziale di 7.3 × 10−3 𝑘𝑔 𝑚3⁄ al valore di concentrazione dimezzato, proprio come volevasi dimostrare. 16 Giovanni Coan Sebastiano Troian 3.3 PROCESSO PURAMENTE DIFFUSIVO CON RILASCIO ISTANTANEO E PUNTUALE TRA DUE PARETI Considerando il problema puramente diffusivo con rilascio istantaneo e puntuale tra due pareti, facciamo riferimento al caso pratico del rilascio di benzene da parte di una cisterna in un canale navigabile. In questo caso, però, studiamo la dispersione dell’inquinante nella sezione trasversale al canale. In questo caso i parametri fisici e le variabili cambiano, in quanto la diffusività trasversale è completamente diversa da quella in direzione longitudinale e il dominio è delimitato da due superfici, quella del pelo libero del canale (trascuriamo eventuali evaporazioni) e quella del fondale (trascuriamo eventuali infiltrazioni). Dopo aver fissato le condizioni iniziali, per rappresentare la presenza di una parete “impermeabile” poniamo le condizioni al contorno di flusso nullo 𝐷 𝜕𝐶 𝜕𝑥 = 0. condizioni al contorno di flusso nullo sulle pareti laterali 19 Giovanni Coan Sebastiano Troian 3.4 PROCESSO PURAMENTE DIFFUSIVO CON RILASCIO ISTANTANEO E PUNTUALE CON DECADIMENTO Consideriamo ora il caso di un contaminante rilasciato istantaneamente, si diffonde e decade nel tempo. L’equazione generale che regola il fenomeno è la seguente: 𝜕𝐶 𝜕𝑡 = 𝐷 𝜕2𝐶 𝜕𝑥2 − 𝐾𝐶 con: 𝐾: coefficiente di decadimento della concentrazione 1 [𝑠]⁄ ; Nel caso specificato, con rilascio centrato nel dominio, la soluzione analitica è del tipo: 𝐶(𝑥, 𝑡) = 𝑀 √4𝜋𝐷𝑡 exp (− 𝑥2 4𝐷𝑡 − 𝐾𝑡) dove, rispetto al caso iniziale si è aggiunto solamente il termine di decadimento. Allo stesso modo, anche la soluzione numerica rimane inalterata a meno del termine di decadimento. Le condizioni iniziali, definite come soluzione analitica al tempo 𝑡𝑖𝑛 = 10 ∗ 𝑑𝑡 variano, mentre le condizioni al contorno di parità dei flussi rimangono inalterate. Di seguito riportiamo i risultati ottenuti. Figura 3.13: Condizioni iniziali processo diffusivo con rilascio istantaneo e puntuale e con decadimento 20 Giovanni Coan Sebastiano Troian Figura 3.14: Andamento concentrazione processo diffusivo con rilascio istantaneo e puntuale e con decadimento 21 Giovanni Coan Sebastiano Troian 3.5 PROBLEMA DI DIFFUSIONE ANISOTROPA IN DOMINIO 3D Consideriamo il caso di un processo di pura diffusione nelle tre direzioni spaziali e, come nella maggior parte dei sistemi ambientali, consideriamo la diffusione anisotropa, in quanto non procede alle stesse velocità nelle direzioni orizzontali e verticali. In quest’ultima direzione, infatti, le scale del moto sono minori e quando si ha stratificazione stabile la turbolenza viene inibita. Potremmo considerare i coefficienti di diffusione diversi in tutte le direzioni, ma per semplificare la rappresentazione del problema su Mathlab consideriamo 𝐷𝑥 = 𝐷𝑦. Avremmo quindi: 𝑞𝑥 = −𝐷 𝜕𝐶 𝜕𝑥 𝑞𝑦 = −𝐷 𝜕𝐶 𝜕𝑦 𝑞𝑧 = −𝐷𝑧 𝜕𝐶 𝜕𝑧 Consideriamo il problema di un rilascio istantaneo e puntuale vicino ad una parete, per esempio il caso di una esplosione che avviene ad una altezza H dal suolo. Lavorando con il metodo delle immagini otteniamo la seguente soluzione analitica: Considerando la diffusività isotropa sul paino x-y la distribuzione della concentrazione su di esso avrà andamento circolare. Su Mathlab abbiamo semplificato la rappresentazione della soluzione considerando il piano x-z, in quanto sarà simmetrica rispetto al piano y-z. La condizione di stabilità dell’algoritmo, in questo caso, sarà quella più restrittiva tra la direzione x e la direzione z, per questo motivo è stato scelto il valore di dt più piccolo. Di seguito riporto le condizioni al contorno e il ciclo spaziale e temporale utilizzato, con le istantanee dell’andamento della distribuzione a diversi istanti di tempo. 24 Giovanni Coan Sebastiano Troian 3.6 PROCESSO DI ADVEZIONE – DIFFUSIONE CON RILASCIO ISTANTANEO E PUNTUALE Passiamo a considerare il caso in cui il fluido all’interno del quale avviene la diffusione si muova lungo una direzione preferenziale con una certa velocità 𝑢. Come si può notare dall’immagine riportata qui sopra, per quanto riguarda la definizione dei parametri, rispetto al caso precedente dobbiamo considerare la velocità media del fluido. Invece, la novità riguarda l’introduzione di una nuova condizione di stabilità: 𝐶𝐹𝐿 = 𝑈 ∆𝑡 ∆𝑥 < 1 Il valore di 𝑑𝑡 finale sarà il minimo tra quello relativo al limite di stabilità diffusivo (𝑃𝑒 ≪ 1) e al limite advettivo (𝑃𝑒 ≫ 1). Per quanto riguarda le condizioni iniziali, considerando un rilascio puntuale e istantaneo, si assumono pari alla soluzione analitica calcolata in mezzeria al dominio. Come è già stato effettuato nel caso speculare di pura diffusione. Figura 3.18: Condizioni iniziali caso processo advettivo-diffusivo con rilascio istantaneo e puntuale Condizione di stabilità diffusiva Condizione di stabilità advettiva 25 Giovanni Coan Sebastiano Troian Il risultato grafico, stampato ad istanti temporali ridotti rispetto a quelli di calcolo per consentirne una migliore visualizzazione, è il seguente: Figura 3.19: Andamento concentrazione nel tempo, nel caso di un processo advettivo-diffusivo con rilascio puntuale e istantaneo Per determinare la bontà della soluzione trovata confrontiamo i risultati della simulazione numerica con quelli ottenuti attraverso la soluzione analitica: 𝐶(𝑥, 𝑡) = 𝑀 √4𝜋𝐷𝑡 exp (− (𝑥 − 𝑢𝑡)2 4𝐷𝑡 ) 26 Giovanni Coan Sebastiano Troian 3.7 PROCESSO DI ADVEZIONE – DIFFUSIONE CON RILASCIO ISTANTANEO E PUNTUALE CON DECADIMENTO Vediamo ora la situazione di advezione – diffusione, come il caso precedente, ma con la presenza anche del termine di decadimento. Le variazioni nella scrittura del codice, rispetto al caso presedente, sono le seguenti:  termine di decadimento nelle condizioni iniziali  termine di decadimento nell’algoritmo numerico Il risultato, stampato ad istanti temporali ridotti rispetto a quelli di calcolo, per consentirne una miglior visualizzazione, è il seguente: Figura 3.20: Distribuzione concentrazione processo advettivo-diffusivo con rilascio istantaneo e puntuale e termine di decadimento 29 Giovanni Coan Sebastiano Troian La simulazione ha quindi fornito i seguenti elaborati grafici (dinamici poiché in funzione del tempo) e per semplicità vengono riportate alcune istantanee del processo. 30 Giovanni Coan Sebastiano Troian Figura 3.22: Istantanee processo di advezione - diffusione bidimensionale 31 Giovanni Coan Sebastiano Troian