Fundamentos fluidos computacionales - MEC-223, Nicolas Thiers Tarea 1 - Soluci´on on num´erica eric a de la ecuaci´ ecua ci´on on de Navier-Stokes Universi Uni versidad dad T´ ecnica ecni ca Federico eder ico Santa Mar´ıa ıa Departamento de Mec´anica anica Fernan Gonz´ alez alez 18 de octubre de 2016
1.
Intr Introd oduc ucci ci´ on o ´n
Se resolver´a la ecuaci´on on de Navier-Stokes (N-S) para el caso de un flujo encerrado entre dos placas paralelas, el cual se mueve debido a la gravedad.
Figura 1: Flujo empujado por la gravedad entre dos placas paralelas Para implementar el c´odigo odigo a la soluci´on on de este problema se utilizara util izara el m´ etodo etodo de proyecci´on on correcci´on, on, por la cual se separara la ecuaci´on on general. ∂V ∇ p + (V · ∇)V = − + ν ∇2 V + g ∂t ρ
(1)
El cual consiste en utilizar una velocidad auxiliar V , que representa rep resenta una velocidad velocid ad intermedia intermed ia cuyas caracter´ car acter´ısticas ıstica s son equivalentes a las del vector velocidad normal V , para calcular en primer lugar ∗
V − V n = −(V n · ∇)V n + ν ∇2 V n t ∗
(2)
La cual corresponde corresponde a proyec proyecci´ ci´ on on y luego la ecuaci´on: on: +1 V n +1 − V ∇ pn =− ρ t ∗
1
(3)
Universidad T´ ecnica Federico Santa Mar´ ıa Departamento de Mec´anica
Que corresponde a la correcci´on. Si a esta ecuaci´on se le aplica un operador ∇ y se re-ordenan t´erminos se encuentra que ∇ · V n +1 el cual representa la divergencia de la velocidad cuyo valor es nulo debido a la conservaci´on de masa y si consideramos al flujo como incompresible, entonces se obtiene: 2 n+1
∇
=
p
ρ t
(4)
∗
∇ · V
Es decir, el objetivo es encontrar el vector V con la ecuaci´on (2) con informaci´on de un vector V n , el cual en una primera iteraci´on temporal sera nulo. Luego al tener el vector V se proceder´a a encontrar la presi´on p con la ecuaci´on (4) y finalmente con la ecuaci´on (3) se obtiene el vector V n +1 que corresponde a la velocidad a un paso temporal posterior. Las EDP necesitan de condiciones de borde para ser resueltas, en este problema se considerara la condici´on de adherencia para las paredes que fija como nula a las velocidades. Para la presi´on tambi´en es necesaria la condici´on de borde, se considera la ecuaci´on (3) y la evaluamos en el borde aplicando el producto punto con el vector normal n a la superficie. ∗
∗
+1 V n +1 · n − V · n ∇ pn n =− · ρ t ∗
(5)
Por la condici´on de no deslizamiento en las paredes se tiene que V · n = 0 lo que implica la aparici´o n de una condici´ on de borde para la presi´on, de la forma: ∂p n+1 =0 ∂n
(6)
Ademas se considerara condiciones de periodicidad para la velocidad y la presi´on en la parte superior e inferior del flujo.
2.
Esquemas num´ ericos
Se analizara un flujo bidimensional cuya velocidad sera denotada V = (u, v ), de la ecuaci´on (2) entonces se puede escribir n u − un ∂u n ∂ 2 un ∂ 2 un n ∂u = −un + ν + −v ∂x ∂y ∂x 2 ∂y 2 t ∗
n v − vn ∂v n ∂ 2 v n ∂ 2 v n n ∂v = −un + ν + −v ∂x ∂y ∂x 2 ∂y 2 t ∗
(7)
−g
u = 0, v = 0 en x = 0, x = H u ,v periodicos en y = 0, y = L ∗
∗
∗
∗
Se procede entonces a discretizar estas ecuaciones con diferencias finitas atrasadas para la primera derivada y diferencias centradas para la segunda derivada, al re-ordenar para dejar todo en funci´on de V n se obtiene
ui,j = u n i,j + t ∗
n vi,j = v i,j + t ∗
n − ui,j
n − ui,j
n un i,j − ui
1,j
−
x
n vi,j − vin
1,j
−
x
n − vi,j
n − vi,j
n un i,j − ui,j
y
n n vi,j − vi,j
1
−
y
+ + 1
−
ν
n n un i+1,j − 2ui,j + u i
1,j
−
(x)2 n vin+1,j − 2vi,j + vin ν (x)2
1,j
−
+
n n un i,j +1 − 2ui,j + ui,j
(y)2 n n n vi,j +1 − 2vi,j + v i,j + (y )2
1
−
1
−
−g
(8)
Luego se procede con la ecuaci´on (4) de forma diferencial ∂ 2 pn+1 ∂ 2 pn+1 ρ ∂u ∂ v + = + 2 2 ∂x ∂y ∂y t ∂x
LATEX
∗
∗
(9) 2
Universidad T´ ecnica Federico Santa Mar´ ıa Departamento de Mec´anica
∂p n+1 = 0 en x = 0, y = L ∂x pn+1 peri´ odico en y = 0, y = L
Se puede notar que esta ecuaci´on al lado derecho corresponde a una EDP el´ıptica, y el lado derecho corresponde a la divergencia de la velocidad intermedia V . Esta EDP se resuelve llev´andola a un sistema matricial A · p = f (x) + b, donde A corresponde a una matriz de dimensiones nx − 2 × ny − 2 con nx y ny los nodos de la malla que se esta utilizando. p corresponde a un vector de tama˜no n x − 2 · ny − 2 , f (x) contiene informaci´on de las fuentes que inciden en el comportamiento de p, en este caso a la divergencia de V y b es un vector que contiene informaci´on tanto de las condiciones de borde como de las esquinas de la malla. Este es un sistema de ecuaciones denominado Poisson, con el cual se obtendr´a un vector p (ver referencia 3) ∗
∗
Luego de obtener el vector de presi´on el algoritmo utiliza una funci´on que transforma este vector en una matriz y le aplica las condiciones de borde correspondientes. De la ecuaci´on de correcci´on (3) cuya forma diferencial es: un+1 − u 1 ∂p n+1 = ρ ∂x t n+1 v 1 ∂p n+1 −v = ρ ∂y t ∗
(10)
∗
un+1 = 0, v n+1 = 0 en x = 0 y x = H un+1 , v n+1 peri´ odico en y = 0, y = L
Se procede a discretizar la presi´on con diferencias atrasadas, obteniendo finalmente
+1 un = u i,j + i,j ∗
n+1 vi,j
= v i,j + ∗
t x · ρ t y · ρ
+1 n+1 pn i,j − pi 1,j −
+1 n+1 pn i,j − pi,j −1
(11)
El algoritmo consiste en copiar el vector V n +1 obtenido con la ecuaci´on anterior, para calcular V y repetir el proceso en n t iteraciones temporales. Cabe destacar que la para el m´ etodo de P oisson la matriz A se mantiene constante en toda iteraci´on temporal por lo que no es necesario re-calcularla, lo que permite ahorrar tiempo. Ademas esta matriz contiene una gran cantidad de valores nulos, por lo que se utiliza un formato de compresi´on llamado CSR (compressed storage row ) y una librer´ıa de phyton llamada scipy que permite resolver sistemas matriciales con estas matrices CSR de forma mas eficiente (ver referencia 4). Entonces en resumen el algoritmo realiza las siguientes operaciones: ∗
1. Se encuentra u y v con (8) para un V i nicial = 0 ∗
∗
2. Se calcula la divergencia de V
∗
3. Se resuelve la ecuaci´ on de Poisson para obtener la presi´on 4. Se aplican las condiciones de borde para la presi´ on 5. Se calcula V n +1 con (11) 6. Se copia el vector V n +1 como V inicial y se repite el proceso hasta un numero nt de iteraciones 7. Se crea un gr´afico donde el fondo representa la presi´on y los vectores representan la velocidad del flujo LATEX
3
Universidad T´ ecnica Federico Santa Mar´ ıa Departamento de Mec´anica
3.
Soluci´ on anal´ ıtica
Para poder comparar los resultados obtenidos por el m´ etodo num´erico se necesita de la soluci´on anal´ıtica de la ecuaci´ on de N-S para las condiciones del flujo dado. Se tiene entonces la ecuaci´on en su forma diferencial para el eje x e y ademas de la ecuaci´on de conservaci´on de masa ∂u ∂u ∂u ∂ 2 u ∂ 2 u 1 ∂p + u + v = + ν + 2 ∂t ∂x ∂y ρ ∂x ∂x 2 ∂y ∂v ∂v ∂v 1 ∂p + u + v = ∂t ∂x ∂y ρ ∂y
+
∂ 2 v ∂ 2 v ν + 2 ∂x 2 ∂y ∂u ∂ v + =0 ∂x ∂y
(12)
Esta ecuaci´on se resolver´a para un flujo estacionario, es decir cuando el perfil de velocidad se estabilice y no existan
∂ = 0. Como el fluido es empujado solo por la gravedad la cual apunta ∂t en el eje −y , considerando ademas que el flujo es laminar se puede decir que u(x, y ) = 0, lo que inmediatamente implica que todas sus derivadas sean nulas tambi´ en y que ademas no exista una diferencia de presi´ on en el eje x ∂p por lo tanto =0 ∂x ∂v De la conservaci´on de masa entonces se encuentra que = 0. Con todas estas consideraci´on se encuentra que la ∂y presi´ on solamente depende del eje y y que la EDP que modela el problema es:
cambios con respecto al tiempo, por lo tanto
∂p ∂ 2 v + ρ · g = µ 2 ∂y ∂x
Las condiciones de borde son v (x = 0) = 0 y v (x = L ) = 0. Ademas se sabe que la presi´on es del orden 10 ∂p que permite realizar la aproximaci´on de que = 0. As´ı la soluci´ on finalmente queda: ∂y ρg 2 v ( x) = (x − xL) 2µ
LATEX
(13) 13
−
, lo
(14)
4
Universidad T´ ecnica Federico Santa Mar´ ıa Departamento de Mec´anica
4.
Soluci´ on num´ erica
Al implementar el algoritmo en python se obtienen los siguientes resultados para distintos pasos de tiempo
(a) nt=10
(b) nt=50
Figura 2: Soluci´ on num´erica para 10 y 50 iteraciones, vectores de velocidad(flechas) y presi´on(color de fondo)
(a) nt=100
(b) nt=1000
Figura 3: Soluci´ on num´erica para 100 y 1000 iteraciones, vectores de velocidad(flechas) y presi´on(color de fondo)
LATEX
5
Universidad T´ ecnica Federico Santa Mar´ ıa Departamento de Mec´anica
Se puede notar f´acilmente que el flujo se estabiliza entre nt = 100 y nt = 1000, por que los vectores de velocidad cualitativamente son muy parecido. Para analizar un tiempo superior, se realizo una simulaci´on para n t = 5000
Figura 4: Soluci´ on num´erica para nt = 5000 Cuantitativamente se utilizara la norma 2 para obtener el error de la simulaci´on num´erica con respecto a las distintas iteraciones
Figura 5: Error porcentual en funci´on de las iteraciones temporales Se puede notar que el error aumenta a medida que el tiempo aumenta, se puede concluir que el flujo se estabiliza cuando el error es menor a un 10 %.
4.1.
Efecto del numero del Reynolds
El numero de Reynolds se expresa como: LATEX
6
Universidad T´ ecnica Federico Santa Mar´ ıa Departamento de Mec´anica
Re =
V prom · D ν
(15)
Donde la velocidad promedio se puede obtener por medio del caudal Q que fluye a trav´ es de las placas con un ancho L x = 2[m] y una profundidad unitaria, ademas D es una dimensi´on caracter´ıstica del problema, en este caso D = L x ρg Q = V prom · A = 2µ
L
x
(x2 − xLx ) dx · 1
(16)
0
Lo que se obtiene finalmente es V prom =
ρgL 2x 12µ
(17)
Te´ oricamente el numero de Reynolds sirve para poder ver la transici´on entre el flujo laminar y un flujo turbulento. El valor del numero de Reynolds que marca la transici´on puede tomarse como Re = 2300, cuando un flujo tiene un bajo Reynolds generalmente es unidimensional y a medida que la velocidad aumenta el fluido no es capaz de retener la cantidad de movimiento que lleva el flujo, visto de otra forma, la energ´ıa que tiene el fluido no es capaz de mantenerse en una dimensi´on, entonces si el Reynolds aumenta el fluido tender´ıa a ser bidimensional y posteriormente tridimensional.
Figura 6: Forma de la velocidad para flujo laminar y el flujo turbulento La soluci´ on anal´ıtica (14) es para un flujo estacionario y laminar, sin embargo los efectos de la turbulencia vuelven al flujo un problema transitorio, donde las aproximaciones echas para su obtenci´on no son aplicables. Se espera entonces que al aumentar el numero de Reynolds el error relativo a la soluci´on anal´ıtica aumente. Cabe destacar que la ecuaci´on (17) solo sirve para calcular el Reynolds de la soluci´on anal´ıtica, para el caso de la soluci´ on num´erica se calculara un promedio de la velocidad para realizar el an´alisis.
4.2.
Estabilidad
La ecuaci´on de N-S y la forma que se esta estudiando en este caso tiene la caracter´ıstica de ser un problema de conveccion-difusi´ on, debido a que el termino ν · ∇2 V corresponde a la difusi´on de la cantidad de movimiento del fluido y el termino (V · ∇)V corresponde a al conveccion del fluido. Cuando se realiza un esquema num´ erico de este tipo de EDP aparecen dos n´umeros que afectan a la estabilidad y convergencia, la primera es la que aparece en la ecuaci´on (8) al lado izquierdo y que corresponde al termino
LATEX
7
Universidad T´ ecnica Federico Santa Mar´ ıa Departamento de Mec´anica
V · t 1 que se conoce como numero de Courant el cual te´oricamente debe ser C r < . El segundo que tambi´en 2 x t aparece en la ecuaci´on (8) que corresponde a ν 2 . x C r =
5.
Referencias
[ 1 ] Chorin, A. J. The numerical solution of the Navier-Stokes equations for an incompressible fluid. Bull. Am. Math. Soc. 73 928-931 (1967) [ 2 ] Chorin, A. J. Numerical Solution of the Navier-Stokes Equations. Math. Comp. 22 745-762 (1968) doi:10.1090/s00255718-1968-0242392-2 [ 3 ] Laboratorio 6 Fund. fluidos computacionales http://ramos.mec.utfsm.cl/mod/resource/view.php?id=35589 [ 4 ] Recomendaciones de python Fun. fluidos http://ramos.mec.utfsm.cl/file.php?file= %2F372 %2FTareas %2Frecomendaciones computacionales
LATEX
8