TFG-3708TIERRAGÓMEZ,VICENTE

Vista previa del material en texto

Trabajo Fin de Grado
Grado en Ingeniería Aeroespacial
Rotura Autoinducida de Chorros Capilares
Autor: Vicente Tierra Gómez
Tutores: Pedro Ángel Vázquez González
Heliodoro González García
Dpto. de Física Aplicada III
Escuela Técnica Superior de Ingeniería
Universidad de Sevilla
Sevilla, 2021
Trabajo Fin de Grado
Grado en Ingeniería Aeroespacial
Rotura Autoinducida de Chorros Capilares
Autor:
Vicente Tierra Gómez
Tutores:
Pedro Ángel Vázquez González
Profesor Titular
Heliodoro González García
Profesor Titular
Dpto. de Física Aplicada III
Escuela Técnica Superior de Ingeniería
Universidad de Sevilla
Sevilla, 2021
Trabajo Fin de Grado: Rotura Autoinducida de Chorros Capilares
Autor: Vicente Tierra Gómez
Tutores: Pedro Ángel Vázquez González
Heliodoro González García
El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profe-
sores:
Presidente:
Vocal/es:
Secretario:
acuerdan otorgarle la calificación de:
El Secretario del Tribunal
Fecha:
A mi familia,
a Sevilla
y todos los amigos que me ha brindado.
I
Resumen
El presente estudio tiene como objetivo la identificación de regímenes cuasi-periódicos enla generación de gotas de un chorro capilar sin estimulaciones externas. Para ello, se ha
hecho uso del software Basilisk con la intención de crear un módelo numérico capaz de simular
el comportamiento del fluido. Mediante múltiples barridos, se han podido esclarecer resultados
cualititavos y cuantitativos de cuáles son las variables que consiguen el efecto deseado. Asimismo,
la similitud que pueda existir entre estos resultados y las predicciones teóricas, confirmará la validez
del modelo numérico y las hipótesis inicialmente planteadas.
II
Abstract
The study that lies in this document, has for goal the identification of single droplets under aquasi periodic regimen in a capillary jet without any external stimulation. For this purpose,
Basilisk software has been used in order to build a numerical model able to simulate the fluid
behaviour. Through multiple sweeps, qualitative and quantitative results had been clarified as mean
to obtain the variables that achieve the desired effect. Likewise, the similarity that may exist between
these results and the theoretical predictions, will confirm the validity of the numerical model and
the hypotheses initially made.
III
Índice Abreviado
Resumen II
Abstract III
Índice Abreviado IV
Notación VII
Índice de Figuras IX
1 Introducción 1
1.1 Antecedentes históricos 1
1.2 Aplicaciones actuales 2
1.3 Objetivo 3
1.4 Software Basilisk 3
1.5 Software MATLAB 4
2 Descripción del problema 6
2.1 Descripción física 6
2.2 Condiciones iniciales y de contorno 8
3 Rutinas Empleadas 12
3.1 Basilisk 12
3.2 MATLAB 23
4 Análisis de los resultados 29
4.1 Para distintos números de Weber y número de Ohnesorge 0.003 39
5 Análisis de la influencia de la viscosidad 45
6 Conclusiones y Trabajos Futuros 55
Apéndice A Códigos Basilisk 58
A.1 Código Basilisk Inicial 58
A.2 Código Basilisk Definitivo 65
A.3 Código Basilisk Modificado 73
Índice de Tablas 82
Bibliografía 83
IV
Índice
Resumen II
Abstract III
Índice Abreviado IV
Notación VII
Índice de Figuras IX
1 Introducción 1
1.1 Antecedentes históricos 1
1.2 Aplicaciones actuales 2
1.3 Objetivo 3
1.4 Software Basilisk 3
1.5 Software MATLAB 4
2 Descripción del problema 6
2.1 Descripción física 6
2.2 Condiciones iniciales y de contorno 8
2.2.1 Condiciones iniciales 8
2.2.2 Condiciones de contorno 9
2.2.3 Ecuaciones adimensionales 10
3 Rutinas Empleadas 12
3.1 Basilisk 12
3.1.1 Rutina Principal 12
3.1.2 Rutinas Secundarias 20
Rutina de Snapshots 21
Rutina de Facets 21
Rutina de Movie y MoviePPM 22
3.2 MATLAB 23
3.2.1 Rutina de representación gráfica 23
3.2.2 Rutina de Transformada de Fourier 26
4 Análisis de los resultados 29
Evolución de la longitud del chorro para distintos números de Weber 30
Comparativa de los valores de ruptura media para distintos
números de Weber 30
V
Índice VI
Para Weber 3 con altura de celda igual a 1.09375 ·10−1 32
Para Weber 3 con altura de celda igual a 9.375 ·10−2 33
Para Weber 4 con altura de celda igual a 9.375 ·10−2 34
Transformadas de Fourier para Weber 2 con altura de celda igual
a 9.375 ·10−2 36
Evolución de xFront para Weber 2 con altura de celda igual a
9.375 ·10−2 37
Evolución de xFront para Weber 3 38
Evolución de xFront para Weber 4 38
4.1 Para distintos números de Weber y número de Ohnesorge 0.003 39
Transformadas de Fourier para Weber 5 y Ohnesorge 0.003 39
Transformadas de Fourier para Weber 6 y Ohnesorge 0.003 40
Transformadas de Fourier para Weber 7 y Ohnesorge 0.003 41
Transformadas de Fourier para Weber 8 y Ohnesorge 0.003 42
5 Análisis de la influencia de la viscosidad 45
Evolución de la longitud del chorro para Weber 4 y distintos
Ohnesorge 46
Comparativa de los valores de ruptura media para distintos
números de Ohnesorge 46
Evolución de la longitud del chorro para Weber 4 y Ohnesorge 0.001 47
5.0.1 Para Weber 4 y número de Ohnesorge 0.001 48
Transformadas de Fourier para Weber 4 y Ohnesorge 0.001 48
5.0.2 Para Weber 4 y número de Ohnesorge 0.003 49
Transformadas de Fourier para Weber 4 y Ohnesorge 0.003 49
5.0.3 Para Weber 4 y número de Ohnesorge 0.005 50
Transformadas de Fourier para Weber 4 y Ohnesorge 0.005 50
5.0.4 Para Weber 4 y número de Ohnesorge 0.009 51
Transformadas de Fourier para Weber 4 y Ohnesorge 0.009 51
5.0.5 Para Weber 4 y número de Ohnesorge 0.013 52
Transformadas de Fourier para Weber 4 y Ohnesorge 0.013 52
6 Conclusiones y Trabajos Futuros 55
Apéndice A Códigos Basilisk 58
A.1 Código Basilisk Inicial 58
A.2 Código Basilisk Definitivo 65
A.3 Código Basilisk Modificado 73
Índice de Tablas 82
Bibliografía 83
Notación
∇ Gradiente
∂ Derivada parcial
© Copyright
~ex Vector unitario en la dirección axial
H Altura de celda
~n Vector normal a la interfase del chorro
Oh Número de Ohnesorge
P Presión
R Radio del chorro
r Coordenada de dirección radial
t Tiempo
t0 Tiempo adimensional
~u Velocidad del chorro
ux Velocidad axial del chorro
uy Velocidad radial del chorro
U0 Velocidad adimensional
We Número de Weber
x Coordenada de dirección axial
δs Función delta de Dirac
κ Curvatura de fase
µ Viscosidad dinámica
µaire Viscosidad dinámica del medio externo
µagua Viscosidad dinámica del fluido bajo estudio
µ̄ Viscosidad dinámica adimensional
ρ Densidad
ρaire Densidad del medio externo
ρagua Densidad del fluido bajo estudio
ρ̄ Densidad adimensional
σ Tensión superficial
φ Variable binaria
VII
VIII
Acrónimos
ETSI Escuela Técnica Superior de Ingeniería de Sevilla
DOD Drop-on-Demand
CIJ Continuous Inkjet
CFD Computational Fluid Dynamics
FFT Fast Fourier Transform
CFL Número de Courant-Friedrichs-Levy
Índice de Figuras
1.1 Métodos de Impresión por Inyección [1] 2
1.2 Longitud de ruptura en condiciones de microgravedad [2] 3
1.3 Logo MATLAB © 1994-2021 The MathWorks, Inc 4
2.1 Régimen de jetting y dripping [3] 6
2.2 Esquema de la geometría del sistema 7
2.3 Zonas del dominio donde se imponen condiciones de contorno 9
3.1 Comando para iniciar la rutina SelfExcitement.tst 12
3.2 Líneas 10-27 de la rutina principal 14
3.3 Líneas 34-44 de la rutina principal 14
3.4 Líneas 47-65 de la rutina principal 15
3.5 Líneas 115-132 de la rutina principal 16
3.6 Líneas 138-166 de la rutina principal 16
3.7 Líneas 168-202 de la rutina principal 17
3.8 Líneas 205-224 de la rutina principal 18
3.9 Ficheros data- generados por la rutina 18
3.10 Líneas 231-250 de la rutina principal 19
3.11 Líneas 254-270 de la rutina principal 19
3.12 Líneas 272-301 de la rutina principal 20
3.13 Líneas 303-324 de la rutina principal 20
3.14 Líneas 326-331 de la rutina principal 21
3.15 Líneas 353-367 de la rutina principal 21
3.16 Ejemplo de archivo Ux.png 21
3.17 Líneas 338-351 de la rutina principal 22
3.18 Líneas 369-387 de la rutina principal 22
3.19 Captura de imagen del archivo .mpg generado 22
3.20 Fichero apto para su lectura en MATLAB 23
3.21 Transformación de vector tras lecturaen matriz deseada 24
3.22 Construcción de vectores y matriz con las variables de interés 24
3.23 Elaboración de gráficas para la evolución de la velocidad 25
3.24 Elaboración de gráfica para la evolución espacial del chorro 25
3.25 Rutina de transformada de Fourier para la evolución de la velocidad 26
3.26 Obtención de velocidad promedio 26
3.27 Rutina de transformada de Fourier para la evolución del chorro 27
4.1 Evolución de la longitud del chorro para distintos números de Weber 30
IX
Índice de Figuras X
4.2 Comparativa de los valores de ruptura media para distintos números de
Weber 30
4.3 Weber 3 con altura de celda igual a 1.09375 ·10−1 32
4.4 Weber 3 con altura de celda igual a 9.375 ·10−2 33
4.5 Weber 4 con altura de celda igual a 9.375 ·10−2 35
4.6 Transformadas de Fourier con Weber 2 y altura de celda igual a 9.375 ·10−2 36
4.7 Evolución de xFront con Weber 2 y altura de celda igual a 9.375 ·10−2 37
4.8 Evolución de xFront con Weber 3 38
4.9 Evolución de xFront con Weber 4 38
4.10 Transformadas de Fourier con Weber 5 y Ohnesorge 0.003 39
4.11 Transformadas de Fourier con Weber 6 y Ohnesorge 0.003 40
4.12 Transformadas de Fourier con Weber 7 y Ohnesorge 0.003 41
4.13 Transformadas de Fourier con Weber 8 y Ohnesorge 0.003 42
5.1 Evolución de la longitud del chorro para Weber 4 y distintos Ohnesorge 46
5.2 Comparativa de los valores de ruptura media para distintos números de
Ohnesorge 46
5.3 Evolución de la longitud del chorro para Weber 4 y Ohnesorge 0.001 47
5.4 Transformadas de Fourier con Weber 4 y Ohnesorge 0.001 48
5.5 Transformadas de Fourier con Weber 4 y Ohnesorge 0.003 49
5.6 Transformadas de Fourier con Weber 4 y Ohnesorge 0.005 50
5.7 Transformadas de Fourier con Weber 4 y Ohnesorge 0.009 51
5.8 Transformadas de Fourier con Weber 4 y Ohnesorge 0.013 52
1 Introducción
“La mente que se abre a una nueva idea, jamás volverá a
su tamaño original”
Albert Einstein
1.1 Antecedentes históricos
El estudio de la física que rodea a los fenómenos relacionados con la tensión superficial, ha
sido objeto de interés durante siglos. Pero no fue hasta principios del siglo XIX cuando comenzó a
formalizarse su estudio, gracias a la labor pionera de Thomas Young [4] y Pierre-Simon Laplace [5].
Pero debido a la escasez de conocimientos y medios para ser investigada, en sus comienzos el
estudio se limitó a los fenómenos estáticos con la aportación de físicos como Eötvös [6], Lenard [7]
y Bohr [8].
Tras esta primera aproximación, de manera totalmente independiente, Savart y Plateau abrieron
las puertas a futuros estudios relacionados con la formación de gotas. Savart logró, a pesar de las
dificultades que suponía la observación de la ruptura de gotas con el ojo desnudo, extraer retratos
del proceso con una precisión asombrosa [9] . Mientras que por el otro lado, Plateau, apoyándose
en técnicas experimentales desarrolladas por él mismo, fue capaz de obtener la relación existente
entre la inestabilidad del chorro y la tensión superficial [10] [11].
τ =
√
r3ρ
σ
. (1.1)
Siguiendo la línea de estudio ofrecida por Plateau, Rayleigh desarrolló el primer método de
predicción de formación de gotas mediante un análisis lineal de estabilidad [12]. Introduciendo la
dinámica de fluidos al proceso de ruptura, definió el “tiempo capilar”, parámetro descrito con la
expresión 1.1, que permite establecer la escala temporal en la que los fenómenos capilares tienen
lugar. Aunque constituyó un gran avance en su contexto histórico y es de gran utilidad para describir
el crecimiento de inestabilidades, estemétodo falla a la hora de predecir la formación de gotas satélite.
Durante años, la herramienta diseñada por Rayleigh era lo único de lo que se disponía a la hora de
analizar el comportamiento de roturas de chorros. Pero esta no permitía comprender la formación de
gotas satélites. Esta preocupación, fue la que propició el desarrollo del "drop weight method" [13];
en el que se define la tensión superficial con el peso de la gota. Además, Edgerton et al. [14],
1
1.2 Aplicaciones actuales 2
profundizaron años más tarde en el estudio realizado por Lenard, redescubriendo y ampliando lo
analizado por el físico.
Tras el artículo presentado por Edgerton et al., no fue hasta finales del siglo pasado que se
consiguió arrojar algo de luz en el estudio de los aspectos no lineales de la ruptura de gotas. Primero
Peregrine [15], con una detallada secuencia de la ruptura y nuevas ideas acerca de las limitaciones
del chorro; y después Kowalewski [16], con imágenes experimentales increíblemente precisas
haciendo uso de la técnica estroboscópica, fueron lo más relevante que experimentó este campo de
estudio el siglo pasado.
Actualmente, con el auge de aplicaciones en las que los fenómenos de ruptura de chorros resultan
de gran utilidad, este estudio se ha convertido en una labor más pragmática que teórica. Dentro de
esta tendencia, no querría desaprovechar la oportunidad de hacer mención a los estudios elaborados
por el Departamento de Física Aplicada III [17] [18] y los proyectos desarrollados por alumnos de
esta escuela bajo el tutelaje del mismo [19] [20].
1.2 Aplicaciones actuales
Tal y como se expuso al final del apartado anterior, han surgido una amplia variedad de aplicaciones
técnicas para las cuales, estos fenómenos de ruptura, resultan de un inestimable provecho. Es por
ello, que aporta una enorme riqueza no perder la perspectiva práctica del estudio realizado en este
documento. Entre ellas, se va a destacar la impresión por inyección, que ha experimentado un auge
en los últimos años.
Esta técnica consiste en la impresión sobre cualquier superficie física mediante gotas generadas a
través de la ruptura del chorro. Dentro de esta aplicación, podemos distinguir dos métodos, el DOD
(de sus siglas en inglés “drop-on-demand”) y el CIJ (acrónimo de“continuous inkjet”) [21].
Figura 1.1 Métodos de Impresión por Inyección [1].
En el DOD, la ruptura del chorro se alcanza mediante la aplicación de un pulso de presión
en el inyector, lo que da lugar a la formación de gotas en la salida de este. Mientras que en el
CIJ, se tiene un chorro continuo que es estimulado tras la salida del inyector a traves de elementos
piezoeléctricos, lo que provoca su ruptura. Este último método, incorpora un canal tras este elemento
1.4 Software Basilisk 3
piezoeléctrico, dando la posibilidad de reciclar el fluido y mejorando así la producción y el control
de las perturbaciones inducidas. Ambos métodos se recogen en la imagen 1.1.
Cabe mencionar, que al tratarse de un campo de estudio aún en una fase muy temprana de su
desarrollo, muchas de las futuras aplicaciones que podría tener son todavía innovadoras ideas
esperando a prosperar. Aunque no sería muy descabellado pensar en futuras aplicaciones dentro del
ámbito de la salud, campo donde la fluidomecánica ha permitido entender numerosos fenómenos
cardíacos [22], entre otros.
1.3 Objetivo
Esta investigación tiene como objetivo la identificación de un régimen de goteo periódico de un
chorro capilar. A diferencia de algunos trabajos previamente realizados en esta materia [19], el
motor de este régimen se trata del chorro, sin necesidad de que intervenga una perturbación externa.
Es por ello, que se trata de seguir la línea de investigación aportada por el artículo “Self-destabilising
loop of a low-speed water jet emanating from an orifice in microgravity” [2], que ha sido de gran
utilidad a la hora de encontrar respaldo experimental a los resultados numéricos hallados.
Figura 1.2 Longitud de ruptura en condiciones de microgravedad [2].
Partiendo de esta premisa, y empleando el software Basilisk, cuya utilización será posteriormente
justificada en mayor detalle, se dispuso a realizar un amplio barrido para identificar bajo qué valores
del número de Weber -parámetro indispensable en los fénomenos de tensión superficial que cuanti-
fica el peso relativo de los efectos inerciales frente a su tensión superficial- se lograba este régimen.
Una vez hallada esta horquilla de valores,se llevaron a cabo diversas simulaciones bajo supuestos
diferentes, con la finalidad de poder discernir si los resultados que se mostrarán a continuación,
son fruto de la física intrínseca del problema planteado o se deben a la casuística de la simulación.
Además de poder establecer relaciones cualitativas acerca de cómo afectan estos parámetros al
desarrollo del fenómeno.
1.4 Software Basilisk
Basilisk es un software de código abierto englobado dentro de lo que se conoce como dinámica de
fluidos computacional, más comunmente denotado como CFD. Trata de mejorar y ser sucesor del
1.5 Software MATLAB 4
conocido software Gerris, y se encuentra desarrollado por los mismos autores [23] [24].
En lo que se refiere a este documento, se ha decidido trabajar con este software por los siguientes
motivos.
• Alta exactitud en los resultados. En particular, se trata de un eficiente software numérico en
lo que se refiere a los problemas de tensión superficial, indispensable en los fenómenos de
chorros capilares.
• Gran flexibilidad en la realización de cálculos. Basilisk cuenta con la posibilidad de delimitar
el estudio a unas cajas, “boxes”, cuyas dimensiones son establecidas por el usuario. Además,
se permite destinar los cálculos de cada caja a un procesador distinto; lo que permite agilizar
en gran medida todo el proceso numérico, además de la gran libertad que le otorga al usuario
para centrar su estudio a una región de interés.
• Libre selección demallado ymallado adaptativo. Durante las fases preliminares del estudio, los
barridos iniciales se realizaron con mallas muy gruesas gracias a esta alternativa, reduciendo
así el tiempo de cálculo. Además, el propio software es capaz de identificar los lugares donde
coexisten ambas fases para poder refinar el mallado donde sea estrictamente necesario, y no
caer en grandes divergencias entre el modelo numérico y el fenómeno físico sin comprometer
los tiempos de simulación.
• Animaciones. El software incorpora diferentes herramientas que generan imágenes y vídeos
de la simulación, permitiendo así visualizar la evolución del chorro.
1.5 Software MATLAB
MATLAB en cambio, es un software de cómputo numérico especializado en la manipulación de
matrices y cálculo algebraico, así como en la representación de datos y funciones [25]. Es esta
última aplicación, junto con la experiencia que he desarollado durante todo mi recorrido académico,
la que ha determinado su uso.
Además de esto, cuenta con una amplia gama de funciones ya preinstaladas que son de gran
utilidad para el estudio en cuestión. Aunque serán analizadas en profundidad en posteriores páginas,
entre estas funciones destaca la FFT, tansformada rápida de Fourier; capaz de llevar al dominio de
la frecuencia las evoluciones temporales de las variables bajo investigación.
Figura 1.3 Logo MATLAB © 1994-2021 The MathWorks, Inc.
1.5 Software MATLAB 5
2 Descripción del problema
“La religión es una cultura de fe, la ciencia es una cultu-
ra de duda”
Richard P. Feynman
En lo que se refiere a este apartado, se va a proceder a realizar una explicación de la física que
se esconde detrás de este estudio, enunciando las hipótesis realizadas, con el fin de establecer un
problema asumible desde el punto de vista del cálculo numérico. Por último, se establecerán tanto
las condiciones iniciales como las de contorno, tratando de que se ajusten lo máximo posible al
problema que tratamos de simular.
2.1 Descripción física
El modelo trata de replicar la evolución de un chorro de radio de salida igual a R, proveniente de un
depósito carente de interés en el estudio, por lo que se ha supuesto infinito para poder despreciar el
ruido que pudiese provocar en los resultados. El motor de este movimiento es la inercia, ya que al
inicio del chorro (conocido como “inlet” en inglés), se le ha impuesto la velocidad. En oposición a
esta fuerza impuesta, el chorro experimentará los efectos de la tensión superficial al avanzar a través
del aire, decelerando consecuentemente el fluido bajo análisis.
Figura 2.1 Régimen de jetting y dripping [3].
6
2.1 Descripción física 7
Es esta competición entre inercia y capilaridad, la que determina la longitud intacta de chorro,
puesto que la tendencia a formar gotas, que disminuyen la superficie líquida, debe luchar contra la
resistencia a los cambios de velocidad que ello requiere. Bajo estos dos efectos toma lugar nuestro
análisis. En el ámbito académico a este régimen se le conoce como “ jetting′′, que a diferencia del
“dripping′′, forma gotas tras la ruptura de una columna fluida. Un ejemplo gráfico de estos dos
regimenes se aprecia en la figura 2.1
Tal y como expuso Rayleigh en su estudio [12], el carácter axisimétrico del problema permite
que las interacciones que existen entre chorro y aire puedan extrapolarse al estudio del problema
bidimensional. Estableciendo el sistema de referencia inercial en la salida del chorro, con el eje x
en el eje de simetría del chorro y con el eje r perpendicular a este, queda definido un sistema como
el que se muestra en la imagen 2.2.
Figura 2.2 Esquema de la geometría del sistema.
La geometría y los parámetros que definen el problema, se ven reflejados también en la imagen
2.2. El radio del chorro R, con su correspondiente velocidad, denotada por~u. Para caracterizar cada
uno de los fluidos, son necesarias las densidades y viscosidades dinámicas de estos, ρ y µ con
sus respectivos subíndices. El vector~n es el vector normal a la interfaz entre el chorro y el medio
externo con dirección al aire.
Previo a la introducción de las ecuaciones que definen el problema bajo análisis, es necesario
enumerar las hipótesis que se han tenido en cuenta a la hora de hacer frente su estudio.
• En primer lugar, se desprecia el efecto de la gravedad de todo el estudio capilar. Para satisfacer
esta hipótesis, se asume que el estudio se da en ambientes de microgravedad.
• El chorro capilar se ha considerado un fluido newtoniano e incompresible que fluye sin
interacción a través de un medio establecido como aire en condiciones ambiente.
• Se asume la ausencia de gradientes térmicos y de concentración.
• La densidad y viscosidad dinámica del fluido en el chorro se han considerado constantes a
efectos de las dos hipótesis anteriores.
2.2 Condiciones iniciales y de contorno 8
∇ ·~u = 0, (2.1)
ρ
∂~u
∂ t
+ρ~u ·∇~u =−∇P+µ∇2~u+σκδs~n, (2.2)
donde la curvatura, κ , definida por κ = ∇ ·~n, P es la presión del fluido y σ su tensión superficial.
Como todo problema fluidomecánico, las ecuaciones que gobiernan el flujo son las conocidas
como ecuaciones de Navier-Stokes, definidas en 2.1 y 2.2. Nótese como en la ecuación 2.2, se ha
hecho uso de la función delta de Dirac, δs. Esta se define de tal manera que toma valor unidad
exclusivamente en la interfase entre ambos, lugar donde los efectos de la tensión superficial tienen
lugar.
∂φ
∂ t
+~u ·∇φ = 0. (2.3)
Para definir por completo el problema, recurrimos a una variable binaria, que se denota por φ ,
llamada trazador. Esta nace de la necesidad de discernir entre las dos fases existentes en el estudio.
Valiendo la unidad cuando se trata del chorro, y cero cuando se localice en el medio exterior, se
define la ecuación 2.3, capaz de caracterizar la evolución de la interfase.
2.2 Condiciones iniciales y de contorno
Para finalizar esta sección del documento, y así cerrar todo lo que concierne a la descripción física
del problema, solo resta plantear las condiciones iniciales y de contorno.
2.2.1 Condiciones iniciales
ux(t = 0) = u0, (2.4)
En primer lugar, la condición inicial exigida para el presente estudio se muestra en 2.4. En esta,
se impone que la velocidad axial inicial sea igual al módulo de la velocidad en la salida del chorro.
Como consecuencia de la búsqueda de un fenómeno axisimétrico, el valor de la velocidad transversal
debe ser nulo.
2.2 Condiciones iniciales y de contorno 9
2.2.2 Condiciones de contorno
Las condiciones de contorno, tomando como referencia la figura 2.3 son
Entrada:
ux=
{
u0 si r ≤ R,
0 si r > R.
ur = 0,
Salida:
∂ux
∂x
=
∂ur
∂x
= 0,
P = 0,
Arriba:
ux = ur = 0.
(2.5)
La frontera inferior corresponde al eje axisimétrico. Las condiciones de contorno a la salida son las
típicas condiciones de salida de flujo.
Figura 2.3 Zonas del dominio donde se imponen condiciones de contorno.
2.2 Condiciones iniciales y de contorno 10
2.2.3 Ecuaciones adimensionales
Finalmente, a la hora de implementar estas ecuaciones en el software numérico, se ha procedido a
tomar variables adimensionales con vistas a generalizar el cálculo. Al ser la capilaridad la que estable-
ce la fuerza dominante, se decide tomar esta como la base de la adimensionalización de las variables.
U0 =
√
σ
ρR
, (2.6)
P0 =
σ
R
, (2.7)
t0 =
R
U0
=
√
ρR3
σ
. (2.8)
De este modo, la velocidad, la presión y el tiempo pasan a seguir las relaciones 2.6, 2.7 y 2.8
respectivamente.
ρ̄(φ) = φ +(1−φ) ρaire
ρagua
, (2.9)
µ̄(φ) = φ +(1−φ) µaire
µagua
. (2.10)
De la misma manera, la densidad y la viscosidad dinámica; dependientes del fluido en cuestión,
pasan a ser definidas como 2.9 y 2.10. Con estas variables, las ecuaciones resultantes se muestran a
continuación.
Las ecuaciones de Navier-Stokes adimensionales son
∇ ·~u′ = 0, (2.11)
ρ̄
∂~u′
∂ t ′
+ ρ̄~u′ ·∇~u′ =−∇P′+Ohµ̄∇2~u′+ κ̄δ ′s~n, (2.12)
La ecuación adimensional para el trazador es
∂φ
∂ t ′
+~u′ ·∇φ = 0. (2.13)
Las condiciones iniciales adimensionales son
u′x(t ′ = 0) =
√
We = β ,
φ =
{
1 r < 1 y x < Li,
0 resto.
(2.14)
El programa elimina en cada paso de cálculo las gotas que se forman y se queda con el volumen mas
grande de líquido. Para asegurarse de que este volumen corresponde siempre al chorro, al inicio se
fija un volumen de líquido en forma de rectándulo de longitud Li. En las simlaciones típicamente
Li = 5.
2.2 Condiciones iniciales y de contorno 11
Las condiciones de contorno adimensionales son
Entrada:
u′x =
{
β si r ≤ 1,
0 si r > 1.
u′y = 0,
φ =
{
1 si r ≤ 1,
0 si r > 1.
Salida :
∂u′x
∂x′
=
∂u′r
∂x′
= 0,
∂φ
∂x′
= 0,
P′ = 0,
Arriba :
u′x = u′r = 0,
∂φ
∂ r′
= 0.
(2.15)
Los parámetros adimensionales que aparecen en las expresiones previas son
We =
√
ρaguaU20 R
σ
, Oh =
µagua√
ρaguaσR
. (2.16)
El número de Weber, We, representa el balance entre la inercia y la tensión superficial. Para We
altos la inercia domina y para We bajos lo hace la tensión superficial. El número de Ohnesorge
representa el efecto de la viscosidad: valores altos corresponden a viscosidad elevada y valores
bajos a flujos poco viscosos. En (2.14) se define el parámetro β como la raíz cuadrada de We.
3 Rutinas Empleadas
Toda la física es o bien imposible o trivial. Es imposible
hasta que la entiendes, y luego se vuelve trivial
Ernest Rutherford
Tal y como se introdujo en el apartado 1 del documento, todos los resultados numéricos han
sido obtenidos gracias al software Basilisk. Mientras que las representaciones gráficas y cálculos
ajenos al problema fluidomecánico se han hallado a través de MATLAB. En esta sección, previa a
la exposición de los resultados obtenidos, se enuncian y hacen comprensibles las rutinas empleadas
en dichos softwares.
3.1 Basilisk
En lo que concierne al código realizado en Basilisk, es necesario mencionar previamente que se
han utilizado dos rutinas. La primera de ellas corresponde al estudio preliminar, en el que se realizó
un barrido para poder identificar el abanico de valores que se encuentra bajo el régimen que ocupa
este análisis. Es por ello, que esta primera rutina no será incorporada aquí, ya que se trata de una
versión más rudimentaria que la que finalmente se ha desarrollado, en la que se hacía uso de un
mallado adaptativo que no presentaba buenos resultados. De todos modos, esta primera versión se
encuentra en el anexo A.1.
3.1.1 Rutina Principal
Dejando de lado esta primera aproximación, se va a proceder a realizar una minuciosa explicación,
línea a línea, de la rutina principal empleada. Esta se encuentra adjuntada en el anexo A.2.
Figura 3.1 Comando para iniciar la rutina SelfExcitement.tst.
A la hora de empezar a correr la simulación, se recurre a la orden que aparece en 3.1. Para poner
en marcha la orden make SelfExcitement.tst, es necesario definir de antemano las variables del
problema que se encuentran en 3.2. Estas se corresponden con,
12
3.1 Basilisk 13
• MPI: Define el número de procesadores que se emplearán para simular la evolución del
chorro.
• NX: Se trata del número de “boxes” que posee la simulación. El código esta elaborado de tal
manera que cada caja tenga asociada un único procesador. Es por esto que MPI y NX deben
de ser iguales.
• DISPLAY: Es una variable relacionada con la visualización del desarrollo del chorro en
la simulación, que toma exclusivamente tres valores. Si se le introduce el valor 1, la rutina
empieza a correr con la posibilidad de que en algún momento posterior puedas recurrir a esta
visualización. Para 0, la simulación comienza sin generar ningun tipo de representación gráfica.
Y para -1, la rutina no se inicia hasta que no se abre la herramienta para la visualización.
• LEVEL: Establece la altura de la caja mínima en el mallado. esta sigue la ley exponencial
mostrada a continuación,
H = NX ·LBOX/2LEV EL siendo H el lado de la celda. (3.1)
• WEBER: Determina el número de Weber que experimenta el chorro bajo estudio.
• TEND: Corresponde con el tiempo total, en segundos, que se quiere simular dentro del
problema.
• DTMAX: Fija el paso en segundos que existe entre cada evento de la simulación.
• DTOUT: Establece el paso también en segundos, de las capturas y los videos de la simulación.
• JACOBI: Si se define como 1, utiliza el método iterativo de Jacobi para resolver el problema.
• LBOX: Impone la longitud de la simulación.
• SNAPSHOTS: Es una variable binaria, definida de tal manera que si toma el valor 1, muestra
en cada instante definido por DTOUT la representación gráfica del chorro.
• MOVIE: Se trata de otra variable binaria, que si se le otorga el valor 1, genera un archivo
.ogv con la evolución del chorro.
• DUMP: Permite guardar los archivos de volcado de memoria si se le da el valor 1.
• FACETS: Si se le da el valor 1, guardará el perfil del chorro para determinados instantes
temporales en un archivo independiente.
• MOVIEPPM: Análoga a la variable MOVIE, solo que en esta ocasión el formato del vídeo
es .mpg, con la finalidad de que cualquier dispositivo pueda reproducirlo.
3.1 Basilisk 14
Figura 3.2 Líneas 10-27 de la rutina principal.
Nótese como no todas estas variables han sido definidas al ejecutar la orden 3.1. Esto se debe a
que algunas de estas han sido ya incorporadas con valores predeterminados, y solo es necesario
recurrir a ellas en el caso de que se quiera modificarlas.
Figura 3.3 Líneas 34-44 de la rutina principal.
Entrando en lo que al propio código se refiere, lo primero que se encuentra es la llamada de las
funciones necesarias para definir el estudio, junto con la igualdad que se expuso en 2.9 y 2.10. Estas
líneas del código se corresponden con la imagen 3.3. Las funciones empleadas en esta rutina son,
• La función grid/multigrid.h permite establecer el mallado que discretiza las variables a
resolver en las ecuaciones diferenciales. La opción de multigrid, nos permite además, dis-
criminar zonas del dominio de estudio donde no es necesario incluir un mallado muy fino,
optimizando así el proceso de resolución numérico.
• La línea de comando donde se define axi.h se emplea en los problemas donde existe simetría
de revolución alrededor de un eje. Es por esto, que nuestra simulación hace uso de esta
función, donde se impone la utilización de un sistema de coordenadas cilíndricas.
• La órden donde se incluye navier-stokes/centered.h es utilizada para la resolución iterativa
de las ecuaciones de Navier-Stokes para fluidos incompresibles. El método utilizado, parecido
al incluido en Gerris, hace uso de un ciclo temporal genérico, un paso temporal límitado por
el número de CFL [26], el esquemade advección Bell-Colella-Glaz [27] para la resolución
de ecuaciones diferenciales y de un solver de viscosidad.
3.1 Basilisk 15
• La orden donde se hace referencia al archivo two-phase.h, se corresponde con la función
empleada en simulaciones donde se tienen dos fluidos inmiscibles. En ella, se define una
variable binaria f que correspondería con la vista en ??.
• Para la línea donde se hace alusión a tension.h, se incluye con el fin de modelar los fénomenos
relacionados con la tensión superficial. Junto con el archivo anterior, los fenómenos en la
interfase quedan totalmente descritos.
• La utilización del archivo view.h viene motivada por la necesidad de obtener representaciones
gráficas de la simulación en cuestión. Todas las renderizaciones de estas representaciones se
realizan con OpenGL.
• La función tag.h permite vincular cada celda dentro de un vecindario con una “etiqueta”, de
este modo se pueden discriminar los distintos volúmenes fluidos presentes en la simulación.
Su implementación responde a la posibilidad de poder eliminar las distintas gotas que se
formasen ajenas al volumen fluido principal y así reducir el ruido que pudiesen provocar estas
perturbaciones.
Tras definir las funciones que se emplean en la rutina, se ha decidido otorgar el valor predetermi-
nado a los parámetros que caracterizan la simulación. Dicha orden se realiza con vistas a agilizar el
proceso de ejecución del código y no impedir el análisis si no se declaran al comenzar la rutina
todas la variables que lo definen. Una ejemplificación se puede ver en la imagen 3.4.
Figura 3.4 Líneas 47-65 de la rutina principal.
Una vez hecho esto, se definen las propiedades físicas de los fluidos junto con las variables
globales del mismo, donde se incluyen la longitud del chorro y las dimensiones del canal a analizar.
Estas líneas del código se corresponden con 3.5. Entre las propiedades de los fluidos se encuentran,
por este mismo orden: la densidades y las viscosidades dinámicas de ambos, la longitud de la
salida del chorro en metros, la tensión superficial del agua, una variable definida como contador
y el número adimesional de Ohnesorge, fundamental en los fenómenos capilares que relaciona la
magnitud que poseen las fuerzas viscosas frente a las de tensión superficial.
3.1 Basilisk 16
Figura 3.5 Líneas 115-132 de la rutina principal.
A continuación de las líneas de código anteriores, se han definido las condiciones de contorno
que gobiernan el problema, como se puede apreciar en la figura 3.6. En primer lugar, se establece un
bucle del tipo “if” en el que se impone la velocidad en el inlet a los puntos que cumplan la condición
de tener menor coordenada transversal que 2 ·R, es decir, que formen parte de lo que se ha definido
como chorro. En caso de que la condición se incumpla, se impone el valor nulo de la velocidad
transversal.
Figura 3.6 Líneas 138-166 de la rutina principal.
Los siguientes tres conjuntos de órdenes, se corresponden con las condiciones de contorno ??,
?? y ??, expuestas en la sección anterior. Las variables u.n, u.t, p y f se corresponden con la
velocidad axial y transversal, la presión y un parámetro que recoge la evolución del perfil del chorro
respectivamente.
Las condiciones de contorno han sido impuestas mediante el uso de las funciones dirichlet y
neumann. estas definen las ecuaciones que su propio nombre indican, la condición de frontera de
Neumann 3.2 y la condición de frontera de Dirichlet 3.3.
{ dy
dx(0) = α1
dy
dx(1) = α2
(3.2)
3.1 Basilisk 17
{
y(0) = α1
y(1) = α2
(3.3)
Una vez establecidas todas las variables necesarias para definir el problema desde un punto de
vista físico, se prosigue con lo que corresponde a la propia estructura del programa que realiza la si-
mulación, que se encuentra en la figura 3.7. Para comenzar, se recurre a la línea habitual de cualquier
rutina que emplee lenguaje en C++ y se inicia el mallado en las líneas contiguas, entre la 177 y la 181.
Figura 3.7 Líneas 168-202 de la rutina principal.
Entre las líneas 183 y 187 se define el valor de las variables físicas del fluido, haciéndose dis-
tinción entre agua y aire mediante la variable binaria. La ultima línea, hace alusión a los flujos de
Stokes, para los cuales los efectos de la inercia son insignificantes frente a la viscosidad. Por ello, al
dar la orden stokes = false, se impide que se omita el término de advección en las ecuaciones de
Navier-Stokes.
Para concluir con la figura 3.7, restan las líneas de la 190 a la 202. Tratan de un bucle con
una condición “if”, en el que se establece que si el valor de la variable identificatoria del proceso
(“pid” de Process ID en inglés) es igual a 0 -que la rutina se encuentre en la primera iteración de
la simulación- se mostrará en pantalla toda la información de color púrpura entrecomillada con el
comando fprintf. De este modo, se crea un encabezado para todos los archivos que se generen de la
simulación.
Para la representación de los valores que adquieren estas variables, es necesario aclarar que
dependiendo de la variable a mostrar, es imprescindible indicar si se trata de un decimal con signo
de un entero, cuya representación sería la de %d\n, de una conversión a coma flotante con signo en
notación científica, que se denota con %e\n, o de una conversión a coma flotante haciendo uso de
la notación que requiera menor espacio, en cuyo caso se usa %g\n.
Realizando una estructura similar para guardar los datos en ficheros, se hace uso de lo mostrado
en la imagen 3.8. En estas líneas, se genera un fichero por cada núcleo que se tiene dentro del
proceso y que recopila los parámetros de interés en el estudio. La orden sprintf, le otorga el nombre
“data-0nºdenúcleo” y presenta el formato que se ve en 3.9.
3.1 Basilisk 18
Figura 3.8 Líneas 205-224 de la rutina principal.
Estas líneas, referentes a la figura 3.8, representan el final del planteamiento inicial, y lo que resta
de rutina se corresponde con los eventos generados para la evolución del chorro en cada instante
temporal.
Figura 3.9 Ficheros data- generados por la rutina.
Las líneas asociadas a los eventos son las que se pueden ver en la imagen 3.10. Se define el tiempo
inicial, igual a 0, y las respectivas variables que se corresponden con las ecuaciones 2.4. Se hace
uso además de la orden boundary para no encontrar excepciones a la hora de realizar ecuaciones
diferenciales. Esto se debe a que, al encontrarse dividido el dominio en celdas, cuando se inicia
un proceso, el valor de las celdas dentro de un vecindario se actualiza, mientras que las que se
encuentran fuera del dominio computacional, denominadas “ghost-cells”, no. Dicha orden arregla
este problema, dando valor a estas ghost cells una vez se ha iniciado el proceso.
Una vez impuestas las condiciones iniciales, se recurre a la función remove droplets, que tal y
como su nombre indica, elimina la aparición de gotas excesivamente pequeñas como para tenerlas en
cuenta en este estudio, evitando así el ruido generado. Así mismo, la función incorpora la posibilidad
de eliminar las burbujas que surjan de estos procesos mediante la orden bubbles=true.
3.1 Basilisk 19
Figura 3.10 Líneas 231-250 de la rutina principal.
Tras esto, se realiza la localización del punto de máxima longitud del chorro, denominado de aquí
en adelante como xfront. Para ello, se emplean las líneas de código que se encuentran en la figura
3.11. En ellas, se realiza un filtro inicial donde se discriminan los valores de f, parámetro entre 1 y
0 que representa el perfil del chorro, menores que la cota THR. Con los valores almacenados de
f donde se tiene el volumen fluido bajo estudio, se halla el máximo y se almacena en la variable
xfront en cada paso temporal.
Figura 3.11 Líneas 254-270 de la rutina principal.
A continuación, se procede a guardar las actualizaciones de las variables de interés. Esto se
corresponde con lo que se puede visualizar en la imagen 3.12. De manera análoga a la vista en 3.7 y
3.8, se establece el cabecero de los ficheros para cuandoel iterante, variable i, se encuentra en el
comienzo.
Posteriormente, se dispone a completar el fichero con las variables que definen el problema. Para
su extracción se han definido 5 puntos. Estos se encuentran a una distancia de uno, dos, tres, cinco
y diez radios respectivamente. Las variables guardadas en el archivo, son el tiempo, el incremento
de tiempo, el número de iteraciones realizadas, la evolución espacial del chorro antes definida, el
número de iteraciones realizadas en cada uno de los tres primeros puntos definidos y sus respectivas
velocidades. Nótese como para el cálculo de las velocidades, se ha hecho uso de la orden interpolate,
que tal y como su propio nombre indica, interpola la función en primer lugar, u.x en esta caso, en
las coordenadas que se indiquen.
3.1 Basilisk 20
Figura 3.12 Líneas 272-301 de la rutina principal.
Para dar por finalizado el análisis de la rutina principal, quedan por definir las líneas correspon-
dientes a 3.13. En estas, se crea un archivo en el que se recogen todas las velocidades axiales del
chorro para toda su longitud. Véase como para cada núcleo se genera un archivo .dat con la orden
sprintf.
Para obtener estas magnitudes, se recurre a un bucle “while”, y se definen las cuatro variables
necesarias para realizar esta función. En primer lugar, x0 se corresponde con el valor inicial de la
celda en cuestión, x1 con el límite de esta, xx con la variable espacial que va a aumentar en cada
iteración y dx con el paso con el que aumenta xx. Descritas estas variables, se introduce la condición
de que xx no supere a x1 y que mientras esto suceda, se guarden los valores de la velocidad y presión
mediante la orden interpolate.
Figura 3.13 Líneas 303-324 de la rutina principal.
3.1.2 Rutinas Secundarias
Lo analizado anteriormente constituye la estructura fundamental de la rutina, en lo que a obtención
de los resultados se refiere. Pero además de lo expuesto, se han incluido rutinas de postprocesado
que merecen ser mencionadas y analizadas. Estas subrutinas son llamadas siempre que a la hora de
iniciar el código se les establezca su valor correspondiente, tal y como se indicó en 3.1. Además,
se va a seguir una estructura similar a la sección anterior, en la que se realizará una meticulosa
explicación de las funciones y órdenes empleadas en la obtención de resultados auxiliares.
3.1 Basilisk 21
Rutina de Snapshots
En lo que se refiere a la subrutina utilizada para la obtención de representaciones gráficas del chorro,
se han empleado las líneas de código que se encuentran en 3.14 y en 3.15. En primer lugar, se
recurre a la función snapshots incorporada en el conjunto view.h antes mencionado. De este modo,
se crea un evento que comienza en el instante inicial y que aumenta a cada paso temporal dentro de
la simulación. Su valor límite se establece al comenzar la rutina.
Figura 3.14 Líneas 326-331 de la rutina principal.
Figura 3.15 Líneas 353-367 de la rutina principal.
En 3.15 se hace uso de la función draw-vof, capaz de representar el volumen fluido, junto con
la orden squares, que permite generar un mapa de color tras establecer dos valores límites. En
particular, estos son los valores máximo y mínimo de la velocidad axial, hallados en las líneas 359
y 360. A continuación, se genera un archivo .png que se actualiza a cada paso temporal estipulado
dentro del código.
En la línea 364, se actualizan estas variables a cada paso de tiempo y la función draw-string
aporta la leyenda a la imagen. Un ejemplo de estos archivos generados es el que se ve en 3.16.
Figura 3.16 Ejemplo de archivo Ux.png.
Rutina de Facets
En las líneas 3.17, se busca almacenar el contorno del chorro mediante la variable f. Para ello, se
crea un archivo con el nombre “facets-nºPID”, que se corresponde con el antes mencionado Process
ID. Después, se genera un escalar, el ff, que se corresponde con un filtro de alta frecuencia de la
variable f, en el que se le otorga el valor nulo a todos aquellos valores que no superen el umbral de
1 ·10−6. De este modo, la salida queda mucho más depurada.
3.1 Basilisk 22
Figura 3.17 Líneas 338-351 de la rutina principal.
Rutina de Movie y MoviePPM
Con la intención de generar algunos archivos de vídeo que hagan acopio de lo simulado, se escriben
las líneas 3.18. En ellas, se ha utilizado dos formatos de video distintos con la finalidad de poder ser
visualizados sin la necesidad de tener instalado el paquete ffmpeg, imprescindible en el visionado
de archivos .ogv.
Figura 3.18 Líneas 369-387 de la rutina principal.
Las líneas que abarcan de las 369 a la 373, se corresponden con las que generan el archivo .ogv,
función viene incomporada dentro del archivo view.h.
Las líneas hasta ahora no analizadas, hacen alusión a la generación del archivo .mpg, que es
reproducible desde el software instalado por defecto en Windows. Para ello, se crea un evento que se
actualiza en cada paso temporal dentro de la simulación hasta el límite impuesto. En él, se programa
el almacenaje de todas las imagenes .ppm al archivo final en formato video. El resultado obtenido
es el que se muestra en la imagen 3.19.
Figura 3.19 Captura de imagen del archivo .mpg generado.
3.2 MATLAB 23
3.2 MATLAB
3.2.1 Rutina de representación gráfica
Una vez se han generado los ficheros con las evoluciones de la velocidad en los puntos antes
mencionados, es necesario representar estos valores frente al tiempo. La estructura que se aprecia
en 3.9 corresponde con los archivos generados por la rutina. En ellos va incorporado, además de los
valores de la velocidad; el valor de la evolución del chorro, xfront, y el valor del tiempo dentro de la
simulación, entre otros.
MATLAB, a la hora de realizar la lectura de un fichero, devuelve un valor nulo si encuentra
caracteres. Como consecuencia de esta incapacidad de discriminar los valores a representar del
encabezado del archivo, es necesario corregir el formato en el que se presentan los resultados. De
este modo, nuestros ficheros deben tomar la forma mostrada en 3.20.
Figura 3.20 Fichero apto para su lectura en MATLAB.
Al introducir este fichero, se realiza la lectura del mismo e inserta los valores en un vector columna,
nombrado en esta ocasión como A, cuyas componentes son todas las variables separadas por un espa-
cio. Para poder tratar todas estas variables con mayor soltura y que se tenga una estructura análoga a
la ofrecida por el fichero inicial, se transforma este vector en una matriz donde cada columna repre-
sente la evolución de todas las variables. Para ello, se emplean las líneas de código mostradas en 3.21.
En ellas, se define un parámetro numérico, denotado como parn, que fija el número de filas
que posee cada columna. Tras obtener este valor, se crea una matriz auxiliar, llamada “aux”, de
coeficientes aleatorios con las dimensiones deseadas, es decir, 12 columnas y parn filas. Con
la ayuda de un bucle del tipo “for”, que recorre el vector A componente a componente, se van
sustituyendo los valores aleatorios en el vector aux por los correspondientes en A. Para finalizar, se
iguala la matriz auxiliar obtenida por A, con el fin de trabajar con la misma notación.
3.2 MATLAB 24
Figura 3.21 Transformación de vector tras lectura en matriz deseada.
Estas líneas de código son ejecutadas en tantas ocasiones como ficheros “data” generados por
la rutina de Basilisk resulten de interés. Esto se debe a que, como se mencionó anteriormente,
los resultados que se obtienen en cada “box” dentro de la simulación corresponden a un fichero
independiente. Así, para este caso particular -con una longitud de caja igual a 4 ·R- los tres primeros
puntos de estudio se encuentran dentro del primer fichero (ya que el punto P3 se encuentra a una
distancia de 3 ·R), el cuarto punto en el segundo y el último en el tercero.
Figura 3.22 Construcción de vectores y matriz con las variables de interés.
El siguiente paso para la representación consiste en la extracción de los datos de interés que
tenemos en la matriz. De esta forma, tal y comose puede apreciar en 3.22, se generan 2 vectores y
una matriz. El vector t se compone de la evolución del tiempo dentro de la simulación, mientras que
el vector x f ront almacena la evolución del chorro en el eje x. Por último, en la matriz v se guardan
las componentes de la velocidad de tal modo que cada columna se corresponde con un punto de
estudio. Cabe mencionar que en la figura 3.22, B y C hacen alusión a las dos matrices asociadas a
los ficheros data para la segunda y tercera caja respectivamente.
Para concluir con la rutina asociada a la representación gráfica, se define el parámetro Beta que
se corresponde con la raíz cuadrada del número de Weber, al igual que se estableció en la rutina de
Basilisk. Este valor se representará junto a las evoluciones de la velocidad para corroborar que los
resultados son coherentes. Haciendo uso de la función plot de MATLAB, se obtienen finalmente
las representaciones deseadas.
3.2 MATLAB 25
Figura 3.23 Elaboración de gráficas para la evolución de la velocidad.
Dado que se han representado los valores en más de una gráfica -como se ve en la imagen 3.23-
se utiliza la función f igure en dos ocasiones. Esta función, junto con los comandos hold on y hold
off, permite la representación de todo lo que se encuentre entre estas dos órdenes en una misma
gráfica. Además, se han incorporado órdenes que facilitan una mejor interpretación de la imagen
resultante.
En primer lugar, la función legend incorpora una leyenda a la representación gráfica. Las órdenes
xlabel e ylabel le añaden una etiqueta al eje en cuestión, otorgando la capacidad de nombrar qué
variables estamos representando. El comando ylim tiene la función de establecer los límites en el
eje de ordenadas a representar. Para concluir, las funciones grid on y grid minor aportan un fino
mallado con el fin de poder localizar mejor en qué orden de valores se encuentran los resultados.
Procediendo de manera análoga para la evolución espacial del chorro -correspondiente a 3.24-
se da por finalizada esta rutina. Veáse cómo en ambas líneas de código, para la representación se
han descartado los primeros valores obtenidos. Esto se debe a que se ha tratado de eliminar las
imposiciones iniciales que se tienen en la simulación. Además de mejorar los resultados de las
velocidades en puntos muy alejados de la salida del chorro en los instantes iniciales.
Figura 3.24 Elaboración de gráfica para la evolución espacial del chorro.
3.2 MATLAB 26
3.2.2 Rutina de Transformada de Fourier
Tanto para la evolución temporal de la velocidad, como para la distancia que recorre el chorro, se
han representado las tranformadas de Fourier con el fin de observar las frecuencias dominantes;
con el fin de encontrar periodicidad en los resultados.
Con vistas a realizar dicha labor para las velocidades, se ha confeccionado el código que se ve en
la figura 3.25. En primer lugar, se definen las variables correspondientes al periodo de la muestra,
que viene impuesto en el código elaborado en Basilisk, y la longitud de la muestra, que será igual a
la variable parn mencionada.
Figura 3.25 Rutina de transformada de Fourier para la evolución de la velocidad.
Tras esto, se elabora un bucle para los cinco distintos puntos donde se realiza el análisis de la ve-
locidad axial. La primera función empleada en este bucle es la fft, capaz de calcular la transformada
discreta de Fourier del vector velocidad usando un algoritmo de transformada rápida de Fourier
implementado en el software.
Nótese que la matriz introducida es velprom y no la matriz vel antes definida. Esta matriz viene
dada por las líneas de código 3.26, en donde se suprime la componente promedio de la velocidad
impuesta por las condiciones iniciales. A pesar de saber que esta componente es igual a Beta, para
eliminar errores de cálculo numérico se ha optado por calcular la media a todo el vector en cuestión y
restar dicha media. Esta operación se ha llevado a cabo con la función mean, que ha sido comparada
en todas las ocasiones con el valor de Beta; confirmando que es el mismo valor con un error de un
orden despreciable.
Figura 3.26 Obtención de velocidad promedio.
3.2 MATLAB 27
Después de realizar el algoritmo de transformada rápida de Fourier en la rutina 3.25, se cálcula
el espectro bilateral, denotado como P2, y el espectro unilateral basado en el espectro bilateral y la
longitud de la entrada, llamado P1 en este código.
Para poder representar la transformada, se define el dominio de la frecuencia f tal y como se
muestra en la línea 98. Una vez determinado, de nuevo se recurre al comando f igure para elaborar
una gráfica para los cinco distintos puntos de estudio. Y con la orden plot, se representa el valor de
la transformada rápida de Fourier de la velocidad en el dominio de la frecuencia. Es remarcable
el hecho de que se han vuelto a utilizar las mismas funciones de diseño de gráficas, con el fin de
aportar mayor clarividencia a la hora de visualizar las representaciones.
Extrapolando esta misma rutina elaborada para la velocidad a la evolución espacial del chorro,
obtenemos el código 3.27. Como diferencia entre ambas rutinas, se debe hacer mención a lo que se
encuentra en la línea 110. En ella, se puede apreciar como la media se toma para todos los valores
del xfront menos para el inicial. Esto se debe a que el primer valor que devuelve la rutina de Basilisk
carece de interés para el estudio.
Figura 3.27 Rutina de transformada de Fourier para la evolución del chorro.
3.2 MATLAB 28
4 Análisis de los resultados
No se puede enseñar nada a un hombre, sólo se le puede
ayudar a descubrirse a sí mismo
Galileo Galilei
La estrategia seguida para la verificación del régimen de autoexcitado se fudamenta en una
hipótesis. Se ha admitido que un indicio de este régimen es la aparición de frecuencias asociadas
a la presencia de un régimen periódico en puntos intermedios del chorro. Para ello, se mide en
distintos puntos la evolución temporal de la velocidad axial en el eje.
En primer lugar, se realizó un descarte preliminar gracias a la retroalimentación con lo que se iba
simulando. En esta sección del apartado, se hizo uso de la rutina A.1, extremadamente simplista, ya
que lo que primaba era la cadencia con la que se obtenían los resultados.
Tras este primer acercamiento hacia el problema bajo estudio, se determinó en una primera
instancia, que la horquilla de valores para los que se apreciaba cierto comportamiento oscilatorio
iba de Weber 2 hasta Weber 4. Para velocidades mayores, no se apreciaron periodicidades en las
magnitudes evaluadas.
29
30
Evolución de la longitud del chorro para distintos números de Weber
Figura 4.1 Evolución de la longitud del chorro para distintos números de Weber.
Comparativa de los valores de ruptura media para distintos números de Weber
Figura 4.2 Comparativa de los valores de ruptura media para distintos números de Weber.
Una muestra de este efecto, se puede apreciar en las representaciones 4.1 y 4.2, realizadas en
una etapa final del estudio, pero colocadas aquí para reafirmar esta criba inicial. En ella, se puede
observar como a partir de números de Weber mayores a 4, el fenómeno de autoexcitado no tiene
lugar y no se provoca la ruptura prematura del chorro. Sin embargo, la comparación entre 1.2 y 4.2
muestra discrepancias entre lo observado por Umemura [2] y lo simulado.
Tras estos primeros análisis, y debido a que esta elección se debía a la mera observación de las
gráficas obtenidas, se dispuso a comprobar qué frecuencias eran las dominantes en la evolución
31
temporal de las velocidades; en el caso de que nuestra suposición de que existía cierta perodiciodad
fuese cierta. Esto, no hizo más que corroborar las suposiciones iniciales.
A la vista de que esto suponía un buen punto de partida para analizar en mayor profundidad la
naturaleza de estos resultados, se optó por aumentar la finura del mallado y el tiempo de exposición, y
de este modo, verificar si el efecto observado se debía a la incapacidad de la simulacióna reproducir
la realidad o si por el contrario, es inherente a la física del problema.
Una muestra de esto, se ve en las figuras 4.3 y 4.4. En ellas, se ha representado la evolución
temporal de la velocidad axial en dos puntos con sus respectivas transformadas de Fourier. Cabe
aclarar que la elección de estos dos puntos no es arbitraria, y se han tomado de manera que el punto
P4 se encontrase a medio camino entre el orificio de salida y la ruptura del chorro, con el fin de
tener unos resultados independientes a las perturbaciones que pudiesen provocar estos; y el punto
P5 cerca de donde se produce dicha ruptura. En el estudio llevado acabo por el físico nipón Akira
Umemura [2], este valor es del orden de 20, como se ve en 1.2.
K = 2πFrec./β . (4.1)
A la vista de las representaciones obtenidas, queda clara la predominancia de una frecuencia
alrededor de 0.25 Hz, que se mantiene con independencia del punto de estudio y el tamaño de celda
utilizado en el mallado. Resulta de gran interés, comparar estos picos con el máximo de Rayleigh,
K, situado en torno a 0.7. Para convertir el valor en el dominio de la frecuencia a la longitud de
onda adimensional, se hace uso de la fórmula en 4.1. Estos valores han sido tabulados en 4.1, en los
que se puede apreciar la cercanía con el valor hallado por Rayleigh.
Tabla 4.1 Valores de la longitud de onda adimensional en P4 y P5.
Frec(P4) K(P4) Frec(P5) K(P5)
Weber 3 H = 1.09375 ·10−1 0.2385 0.8652 0.2385 0.8652
Weber 3 H = 9.375 ·10−2 0.237 0.8597 0.237 0.8597
Weber 4 H = 1.09375 ·10−1 0.3045 0.9566 0.2774 0.8718
Weber 4 H = 9.375 ·10−2 0.312 0.9802 0.285 0.8954
32
Para Weber 3 con altura de celda igual a 1.09375 ·10−1
Figura 4.3 Weber 3 con altura de celda igual a 1.09375 ·10−1.
33
Para Weber 3 con altura de celda igual a 9.375 ·10−2
Figura 4.4 Weber 3 con altura de celda igual a 9.375 ·10−2.
34
Estos mismos resultados para Weber igual a 3, se reproducían para Weber igual a 4, con la
salvedad de que la evolución de la velocidad se encontraba altamente amortiguada y este efecto era
mucho más tenue. Con la finalidad de ejemplificar esta afirmación, se han adjuntado las figuras 4.5,
que de manera análoga a lo representado para Weber 3, se han tomado los puntos P3, localizado
a una distancia de 3 radios del inlet, P4 y P5 para representar la evolución de la velocidad y sus
transfromadas en el dominio de la frecuencia.
Además de esto, se ha incorporado una vista aumentada de dicha evolución de la velocidad para
poder visualizar el comportamiento oscilatorio que presenta. De nuevo, las representaciones de
estas en el dominio de la frecuencia presentan un máximo que es invariable con el lugar donde se
analice, aunque su presencia se mitiga a medida que el análisis se realiza cerca del punto de ruptura
del chorro. Los valores de los máximos de la frecuencia están tabulados en 4.1.
Para Weber 4 con altura de celda igual a 9.375 ·10−2
35
Figura 4.5 Weber 4 con altura de celda igual a 9.375 ·10−2.
36
Sin embargo, a pesar de que para Weber 2 las tendencias se mantenían, a medida que se disminuía
el tamaño de las celdas del mallado se introducía gran cantidad de ruido a las visualizaciones de la
transformada de Fourier de las evoluciones de la velocidad, como se puede comprobar en 4.6.
En un primer momento, se hipotetizó acerca de la posibilidad de que la simulación no caracterizase
correctamente el fenómeno bajo estudio. Pero tras comprobar que los valores de la evolución del
xfront, mostrados en 4.7, estaban respaldados por el estudio de Umemura [2], se teorizó acerca de
la posibilidad de que se tratase del caso límite entre “jetting” y “dripping”. Ante esta hipótesis, se
optó por realizar simulaciones en video y comprobar qué sucedía. Y en efecto, el régimen tratado
se encontraba en este límite. Es por esto, que se decidió descartar esta valor del número de Weber
de este estudio.
Transformadas de Fourier para Weber 2 con altura de celda igual a 9.375 ·10−2
Figura 4.6 Transformadas de Fourier con Weber 2 y altura de celda igual a 9.375 ·10−2.
37
Evolución de xFront para Weber 2 con altura de celda igual a 9.375 ·10−2
Figura 4.7 Evolución de xFront con Weber 2 y altura de celda igual a 9.375 ·10−2.
En este punto del estudio, se encuentra que la simulación es incapaz de aumentar el dominio
computacional sin comprometer el tiempo de cálculo o la altura de las celdas que componían el
mallado. Este obstáculo, era el causante de que no se lograse visualizar la evolución del xfront para
valores mayores que 3, ya que la ruptura se producía fuera del dominio. Con vistas a solventar este
problema se modifica la rutina principal.
La particularidad de esta, es la asignación de celdas del mallado a núcleos independientes. De
esta manera, se puede ampliar la longitud del chorro sin comprometer el tiempo de cálculo. La
modificación a la que hace referencia esta parte se encuentra en A.3. Cabe resaltar que para esta
parte se hizo uso de la computadora del departamento, ya que disponía de 16 núcleos y agilizaba en
gran medida el proceso.
En lo que respecta a Weber 3, estos resultados no arrojan mucha luz al estudio previo, ya que
la altura de la malla es similar a la más fina y el anterior código no imponía ninguna limitación
en el dominio. Sin embargo, este nuevo enfoque permite confirmar lo relativo a la evolución de la
longitud del chorro, apartado que se ha refinado en esta sección. Es por esto, que en este caso, sí
se puede afirmar que los resultados concuerdan con lo hallado experimentalmente [2], tal y como
observamos comparando 1.2 y 4.8.
38
Evolución de xFront para Weber 3
Figura 4.8 Evolución de xFront con Weber 3.
Sin embargo para Weber igual a 4, correspondiente con la imagen 4.9, estos valores se salen de
lo estimado por Umemura [2], ya que deberían rondar los 40, como se ve en 1.2. Este hecho es el
que ha motivado el estudio complementario que se recoge en 5. En él, se trata de evaluar cuál es la
longitud media de ruptura del chorro frente a un parámetro de gran interés en estos fenómenos, el
número de Ohnesorge.
A pesar de esto, los valores de las evoluciones de la velocidad sí que se mantienen dentro de lo
visto hasta ahora. Esto, junto con el análisis que viene a continuación, tratará de justificar la validez
del estudio.
Evolución de xFront para Weber 4
Figura 4.9 Evolución de xFront con Weber 4.
4.1 Para distintos números de Weber y número de Ohnesorge 0.003 39
4.1 Para distintos números de Weber y número de Ohnesorge 0.003
En esta sección del documento, que ha permitido la obtención de las representaciones 4.1 y 4.2,
se va a corroborar que la horquilla de valores tomada inicialmente era correcta y que no se ha
descartado de la criba preliminar ningún resultado que pudiese ser de interés. Para ello, se va
a mostrar la transformada al dominio de la frecuencia de la velocidad axial en dos puntos, uno
previo a la longitud de ruptura y otro tras este, para distintos valores del número de Weber. La
elección de estos dos puntos, tiene como fin comprobar que los efectos del autoexcitado son visibles
exclusivamente en estos puntos previos a la rotura.
Transformadas de Fourier para Weber 5 y Ohnesorge 0.003
Figura 4.10 Transformadas de Fourier con Weber 5 y Ohnesorge 0.003.
4.1 Para distintos números de Weber y número de Ohnesorge 0.003 40
Transformadas de Fourier para Weber 6 y Ohnesorge 0.003
Figura 4.11 Transformadas de Fourier con Weber 6 y Ohnesorge 0.003.
4.1 Para distintos números de Weber y número de Ohnesorge 0.003 41
Transformadas de Fourier para Weber 7 y Ohnesorge 0.003
Figura 4.12 Transformadas de Fourier con Weber 7 y Ohnesorge 0.003.
4.1 Para distintos números de Weber y número de Ohnesorge 0.003 42
Transformadas de Fourier para Weber 8 y Ohnesorge 0.003
Figura 4.13 Transformadas de Fourier con Weber 8 y Ohnesorge 0.003.
A la vista de estos resultados, se puede inferir lo siguiente. Para todos estos números de Weber,
en los puntos previos a la longitud de ruptura media, se encuentran picos alrededorde la frecuencia
que se observa en Weber 4. Sin embargo, para los puntos tras la ruptura, los picos han desaparecido.
Tabla 4.2 Valores de F(vel) y longitud de onda adimensional para distintos números de Weber.
Número de Weber F(vel) Amplitud de la transformada
5 0.000439 0.969548
6 0.000135 1.028715
7 6.2651·10−5 0.517468
8 4.04291 ·10−5 0.521841
En la tabla 4.2, se han recogido los valores de la función en el dominio de la frecuencia y la
longitud de onda adimensional en los puntos previos a la ruptura. Es por esto que, tomando los
valores de esta tabla como apoyo, se concluye que, a pesar de que sí existen frecuencias dominantes
4.1 Para distintos números de Weber y número de Ohnesorge 0.003 43
para estos números de Weber, el valor que toman estas son tan pequeñas que impiden la aparición
del régimen autoexcitado.
4.1 Para distintos números de Weber y número de Ohnesorge 0.003 44
5 Análisis de la influencia de la
viscosidad
Nada en la vida es para ser temido, es sólo para ser
comprendido. Ahora es el momento de entender más, de
modo que podamos temer menos
Marie Curie
En esta sección del estudio, con vistas a dotar de mayor amplitud al análisis, se va comprobar si
al aumentar la viscosidad el comportamiento autoexcitado desaparece. Ya que la viscosidad atenúa
las perturbaciones que se propagan aguas arriba desde la ruptura del chorro.
Por ello, se ha optado por representar para Weber igual a 4 la evolución temporal del punto de
máxima longitud del chorro para distintos valores del número de Ohnesorge. Y de este modo, poder
determinar la consecuencia de los anteriores resultados y comprobar que existe un límite a partir
del cual el fluido es tan viscoso que se erradican los efectos de la autoestimulación.
Resaltar que para todas las simulaciones previas, el valor del número de Ohnesorge era igual a
0.00333105, resultado de tomar el parámetro característico definido como ROD en 3.5, igual a 0.001.
Asimismo, junto a las representaciones mencionadas, se incorporan las respectivas transformadas
al dominio de la frecuencia de la evolución de las velocidades en dos puntos, uno localizado antes
de la longitud de ruptura media y otro después, con el fin de corroborar que este efecto desaparece
si analizamos puntos a la derecha de este valor. La motivación de estas gráficas nace de la necesidad
de poseer mayores argumentos en favor de la veracidad de los resultados hallados.
Con vistas a mejorar la visualización, y consecuentemente, las conclusiones inferidas, se ha
elaborado una imagen que agrupa las evoluciones del xFront para los distintos números de Ohnesorge,
5.1, y la tabla 5.2, en la que se comparan los valores medios de ruptura. La línea vertical denota
la desviación típica de los valores hallados, haciendo así una estimación de lo esperado para los
valores intermedios.
45
46
Evolución de la longitud del chorro para Weber 4 y distintos Ohnesorge
Figura 5.1 Evolución de la longitud del chorro para Weber 4 y distintos Ohnesorge.
Comparativa de los valores de ruptura media para distintos números de Ohnesorge
Figura 5.2 Comparativa de los valores de ruptura media para distintos números de Ohnesorge.
Gracias a estas representaciones, se puede apreciar un comportamiento inicial algo lineal, que
finalmente se amortigua, pudiendo inferir cierta saturación a partir de Ohnesorge 0.007, indicando
así el límite de la viscosidad para la cual se produce el fenómeno autosostenido.
Cabe mencionar, que para valores de Ohnesorge menores a 0.003, el fluido simulado es tan poco
viscoso que se tienen grandes errores en lo que concierne al cálculo numérico. Un ejemplo de esto
se encuentra en la imagen 5.3, correspondiente a la evolución temporal de la longitud del chorro
para Ohnesorge igual a 0.001, valor para el que se agrava el efecto.
47
Esto es consecuencia de que la rutina elimina las gotas que se producen en la ruptura, reteniendo
un único volumen de agua. En la gráfica se observa como el chorro no se rompe y llega hasta el
final del dominio fluido, donde desaparece por completo y vuelve a iniciarse en la salida.
Evolución de la longitud del chorro para Weber 4 y Ohnesorge 0.001
Figura 5.3 Evolución de la longitud del chorro para Weber 4 y Ohnesorge 0.001.
48
5.0.1 Para Weber 4 y número de Ohnesorge 0.001
Dado que para estos valores, el fluido es muy poco viscoso, los resultados reflejan la presencia
de excesivo ruido de perturbaciones. A pesar de tener una señal tan poco filtrada, esta frecuencia
dominante se encuentra presente en el punto previo a la ruptura y desaparece después, tal y como se
puede apreciar en la figura 5.4.
Transformadas de Fourier para Weber 4 y Ohnesorge 0.001
Figura 5.4 Transformadas de Fourier con Weber 4 y Ohnesorge 0.001.
49
5.0.2 Para Weber 4 y número de Ohnesorge 0.003
La principal diferencia encontrada para Ohnesorge 0.003, imagen 5.5, es la eliminación de gran
parte del ruido ocasionado por esta falta de viscosidad. Esto provoca que el pico de frecuencias se
visualice de manera mucho más clara para el punto previo a su ruptura. De nuevo, si centramos el
análisis en puntos tras la longitud de ruptura, este efecto no existe.
Transformadas de Fourier para Weber 4 y Ohnesorge 0.003
Figura 5.5 Transformadas de Fourier con Weber 4 y Ohnesorge 0.003.
50
5.0.3 Para Weber 4 y número de Ohnesorge 0.005
Para esta pareja de representaciones, 5.6, asociadas a Ohnesorge igual a 0.005, que no supera el
límite de saturación mostrado en 5.1, la presencia de este frecuencia dominante en el punto previo al
lugar de ruptura se hace bastante notable. En cambio, para el punto posterior a la ruptura, este efecto
vuelve a desaparecer por completo. Se puede comprobar también, que el efecto de las perturbaciones
producido por la falta de viscosidad se ha erradicado.
Transformadas de Fourier para Weber 4 y Ohnesorge 0.005
Figura 5.6 Transformadas de Fourier con Weber 4 y Ohnesorge 0.005.
51
5.0.4 Para Weber 4 y número de Ohnesorge 0.009
Para el caso de Ohnesorge igual a 0.009, representación 5.7, el valor otorga tanta viscosidad
que provoca que el pico de frecuencias, aunque se localice en el mismo lugar, se encuentra muy
disminuido.
Transformadas de Fourier para Weber 4 y Ohnesorge 0.009
Figura 5.7 Transformadas de Fourier con Weber 4 y Ohnesorge 0.009.
52
5.0.5 Para Weber 4 y número de Ohnesorge 0.013
Al igual que sucedía para el caso anterior, para Ohnesorge 0.013 reflejado en 5.8, el pico en las
frecuencias es significativamente más pequeño que para los casos de escasa viscosidad.
Transformadas de Fourier para Weber 4 y Ohnesorge 0.013
Figura 5.8 Transformadas de Fourier con Weber 4 y Ohnesorge 0.013.
En la tabla 5.1 se recogen los valores de F(vel) de los picos y la longitud de onda adimensional.
Tras todas estas representaciones, se puede concluir que el pico de las frecuencias localizado en
torno a 0.25 Hz, es independiente del valor del Ohnesorge que se emplee. Además, si se supera
el límite observado en 5.1 de 0.007, la viscosidad es tan elevada que el pico en el dominio de la
frecuencia, y consecuentemente el efecto de autoexcitado, se atenúa notablemente.
53
Tabla 5.1 Valores de F(v) para los distintos números de Ohnesorge..
Número Ohnesorge F(vel) Amplitud de la transformada
0.001 0.049683 0.546656
0.003 0.169623 0.766561
0.005 0.000616 0.885951
0.009 0.003047 0.926926
0.013 0.002595 0.923785
Por último, recalcar que la hipótesis de que este efecto desaparecía para puntos localizados tras
el valor de la longitud de ruptura media se ha visto reforzada con estas simulaciones.
54
6 Conclusiones y Trabajos Futuros
“La ciencia aumenta nuestro poder en la medida en que
reduce nuestra soberbia.”
Herbert Spencer
Tras exponer los resultados de manera cronológicamente coherente, y tomando el testigo del
trabajo realizado en esta materia por Daniel Fernández García en “Generación de gotas aisladas
en chorros líquidos” [19], se pueden inferir las siguientes conclusiones acerca del estudio que se
encuentra recogido en estedocumento.
• El estudio se ha beneficiado de una mejora progresiva de los códigos debida al análisis crítico
de los sucesivos resultados, detectando comportamientos no físicos y refinando la elección
de parámetros. Además, el enfoque aportado en la sección 5, no sólo reafirma esta idea, si no
que dota de mayor profundidad al análisis.
• La alta capacidad de computación y fiabilidad del software Basilisk se han visto demostradas.
Todos los resultados numéricos comparados con estudios teóricos o experimentales han
mostrado insignificantes discrepancias.
• La premisa inicial de que existe un régimen bajo el cual el propio chorro induce la generación
gotas de manera periódica, se ha visto reforzada tras los resultados obtenidos. Para valores del
número de Weber entre 3 y 4, se tienen alentadoras previsiones alrededor de este fenómeno.
• El fenómeno de autoexcitado necesita de un chorro que se encuentre totalmente desarrollado.
En vista de lo obtenido para pequeños tiempos de exposición, se ha determinado que estos
resultados son inconcluyentes tras observar los efectos ante largos periodos de simulación.
• La gráfica donde se representa la longitud de ruptura media frente al número de Ohnesorge,
establece un diferenciado cambio régimen ante un exceso de viscosidad. Esta tendencia
muestra una posible transición de autoestimulación a rotura ruidosa con este parámetro, que
debe ser comprobada con estudios posteriores.
Dada la gran complejidad física del problema, es imprescindible mencionar que este estudio
cualitativo debe tomarse como una primera aproximación al problema de chorros autoexcitados;
ya que sería necesario realizar futuros análisis con vistas a obtener mayor exactitud en el modelo
predictivo, además de resultados experimentales que evidenciasen la veracidad de estos resultados
numéricos.
En particular, además de la detección de frecuencias asociadas a modos autosostenidos, explorada
en este trabajo, puede ser de interés la detección de patrones espaciales que muestren la presencia
55
56
del modo inestable que se propaga aguas abajo y la de un modo estable que viaja desde la zona de
rotura hasta el orificio de salida, cerrando el bucle de autoexcitación.
Es por ello, que las conclusiones aquí halladas, esperan servir como referencia y punto de partida
a próximos proyectos que decidiesen abordar este fenómeno. Además de suponer una guía a cual-
quiera que quisiera informarse e iniciarse en el uso de softwares de CFD a través de Basilisk.
57
Apéndice A
Códigos Basilisk
A.1 Código Basilisk Inicial
En esta sección del anexo, se expone el código preliminar utilizado en el barrido inicial para
determinar la zona donde se tiene el régimen de autoestimulación.
/*
* =====================================================================
*
* Filename: selfExcitment.c
*
* Description: Simulation of jets
*
* CC=’mpicc -D_MPI=_ -DDISPLAY=_ -DLEVEL=_ -DADAPT=_ -DLEVELMAX=_ -
DLEVELMIN=_ -DWEBER=_ -DTEND=_ -DDTMAX=_ -DDTOUT=_ -DJACOBI=_ -DLY
=_ -DFILTERED=_ -DSNAPSHOTS=_’ make selfExcitement.tst
* MPI = number of processes
* DISPLAY = 1, play controls, start running immediately.
* -1: play controls, initially paused.
* 0: no play controls, start running
immediately.
* LEVEL = initial refinement level (default 6)
* ADAPT = 1-> adaptive refinement (default 0)
* LEVELMAX = max refinement level (default 7)
* LEVELMIN = min refinement level (default 3)
* WEBER = Weber number (required)
* TEND = computational total time (required)
* DTMAX = maximum size of the time step (default 1.e.3)
* DTOUT = time interval to output data (default 1.0)
* LY = size of the box (default 16.0)
* SNAPSHOTS = 1 -> save snapshots images and simulation (
default 0)
* JACOBI = 1 set Jacobi option for the solver (default
0)
58
A.1 Código Basilisk Inicial 59
* FILTERED = 1 -> takes the average of the tracer at the
corners of each cell (default 0)
* Version: 1.0
* Created: 28/05/19 15:00
* Revision: none
* Compiler: qcc
*
* Author: Pedro A. Vazquez
* Company:
*
*
* =====================================================================
*/
#include "axi.h"
#include "navier-stokes/centered.h"
#define mu(f) (clamp(f,0,1)*(mu1-mu2) + mu2)
#define rho(f) (clamp(f,0,1)*(rho1-rho2) + rho2)
//#define mu(f) (mu1*mu2/(mu1 + clamp(f,0,1)*(mu2-mu1)))
//#define rho(f) (rho1*rho2/(rho1 + clamp(f,0,1)*(rho2-rho1)))
#include "two-phase.h"
#include "tension.h"
#include "view.h"
#include "tag.h"
//#include "axistream.h"
//#include "adapt_wavelet_leave_interface.h"
//#include "drop_stat.h"
#define PI (3.14159265359)
/* Parameters */
#define R 1.0
#ifndef LY
#define LY (16.0)
#endif
#define LX (LY)
#ifndef WEBER
#define WEBER 3.0
#endif
#ifndef BETA
#define BETA (sqrt(WEBER))
#endif
#ifndef TEND
#define TEND 10.
#endif
A.1 Código Basilisk Inicial 60
#ifndef DTOUT
#define DTOUT 1.0
#endif
#ifndef DTMAX
#define DTMAX (1.e-3)
#endif
#ifndef ADAPT
#define ADAPT 0
#endif
#ifndef LEVEL
#define LEVEL 6
#endif
#ifndef LEVELMAX
#define LEVELMAX 7
#endif
#ifndef LEVELMIN
#define LEVELMIN 3
#endif
#ifndef SNAPSHOTS
#define SNAPSHOTS 0
#endif
#define RHOW (1000.)
#define RHOA (1.18)
#define MUW (0.0009)
#define MUA (0.0000185)
#define ROD (0.001)
#define GAMMA (0.073)
#define OH (MUW/sqrt(RHOW*GAMMA*ROD))
/**
Global variables */
double xFront;
double adapt_f = 1.e-2;
double adapt_u = 3.e-2;
/****************************************
C FUNCTIONS
*************************************** */
/**
Boundary conditions for velocity at inlet */
float Ud ( double xp, double yp) {
A.1 Código Basilisk Inicial 61
if ( yp<=4.*R )
return BETA;
else
return 0.;
}
/****************************************
FIELDS AND BOUNDARY CONDIIONS
*************************************** */
// Boundary conditions at inlet
u.n[left] = dirichlet( Ud(x, y) );
u.t[left] = dirichlet(0);
p[left] = neumann(0.);
f[left] = dirichlet( (y<=R));
// Boundary conditions at outlet
u.n[right] = neumann(0);
u.t[right] = neumann(0.);
p[right] = dirichlet(0);
pf[right] = dirichlet(0);
f[right] = neumann(0);
// Boundary conditions at top
u.n[top] = dirichlet(0.); // slip
u.t[top] = neumann(0);
f[top] = neumann(0.);
/****************************************
MAIN PROGRAM
*************************************** */
int main (int argc, char * argv[])
{
DT = DTMAX;
TOLERANCE = 1.e-4;
// Initialize the grid
size (LY);
init_grid (1 << LEVEL);
origin (0, 0);
// fluid parameters
rho1 = 1., mu1 = OH;
rho2 = RHOA/RHOW, mu2 = OH*MUA/MUW;
f.sigma = 1.;
stokes = false;
// print configuration
if ( pid()==0 ) {
A.1 Código Basilisk Inicial 62
fprintf(stdout,"# npe (number of processes)= %d\n", npe());
fprintf(stdout,"# WEBER : %g\n", WEBER);
fprintf(stdout,"# BETA : %g\n", BETA);
fprintf(stdout,"# OH : %g\n", OH);
fprintf(stdout,"# LY : %g\n", LY);
fprintf(stdout,"# DTMAX : %g\n", DTMAX);
fprintf(stdout,"# LEVEL : %d\n", LEVEL);
if (ADAPT!=0) {
fprintf(stdout,"# LEVELMAX : %d\n", LEVELMAX);
fprintf(stdout,"# LEVELMIN : %d\n", LEVELMIN);
}
}
run();
//MPI_Barrier (MPI_COMM_WORLD);
}
/****************************************
EVENTS
*************************************** */
/* Init event */
event init (t = 0) {
//CFL = 0.45;
// initial state
foreach() {
u.x[] = ( BETA );
u.y[] = 0.;
//f[] = ( (y<R) ? 1. : 0.);
f[] = 0.;
}
boundary ((scalar *){f,u});
}
#if 0
event acceleration (i++) {
face vector av = a;
foreach_face(x) {
av.x[] += BO;
}
}
#endif
// remove bubbles
event remove_bubbles_droplets ( i+=1 ) {
A.1 Código Basilisk Inicial 63
remove_droplets (f, minsize = 10, bubbles = true);
//remove_droplets (f, minsize = 10, bubbles = false);
}
#if 0
// Save interface
event interface(i+=10)
{
char name[200];
sprintf(name, "int-%02d",pid());
FILE * fp = fopen(name,"w");
//face vector s[];
output_facets (f, fp);
fclose(fp);
}
#endif
// Compute position of front of main jet
event front (i+=10)
{
scalar m[];
double THR = 1.e-2;
foreach()
m[] = (f[]>THR);
tag (m);
scalar c[];
foreach()
c[] = m[]==1 ? f[] : 0;
boundary({c});
scalar pos[];
position (c,pos,{1.,0.});