Solución de la ecuación de Rayleigh-Plesset por medio del método del elemento finito

Abstract

In this article, the authors have presented numerical solutions of the Rayleigh-Plesset equation which describes the evolution of cavitating bubbles. In order to do that, they consider Finite Element Method Galerkin; this simulation is performed for an inviscid and incompressible fluid in an uniform temperature field with constant surface tension, and the cavitation model into the which the pressure inside bubbles is equal to the fluid vapor pressure. Thus, in this problem is considered the Dirichlet boundary problem, and the authors obtained criteria for the boundary conditions at the cavitation phenomenon through to the which give rise to the bubble growing.

Headnote

Resumen

En este trabajo se plantean soluciones numéricas a la ecuación de Rayleigh- Plesset que describe la evolución de las burbujas en la cavitación. Para ello, se considera el MEFG (Método del Elemento Finito de Galerkin); tal simulación se realiza en un fluido invíscido e incompresible en un campo de temperatura uniforme, una tensión superficial esencialmente constante, y el modelo de ca- vitación en el flujo siendo la presión interna de las burbujas igual a la presión de vapor del fluido. De esta manera, para el problema se considera el problema de Dirichlet y se obtienen los criterios de frontera que auspician el fenómeno de cavitación a través del crecimiento de las burbujas o cavidades.

Palabras clave: cavitación, ecuación de Rayleigh-Plesset, método del elemen- to finito de Galerkin, presión de vapor, fluido invíscido, fluido incompresible.

Aspectos relevantes

* Se plantea una solución numérica a la ecuación de Rayleigh-Plesset por el método numérico del elemento finito de Galerkin (MEFG). * La implementa- ción del MEFG se realizó en tres etapas: el pre-procesamiento donde se modela el problema; la solución de ecuaciones donde se encuentra los valores para la incógnitas del problema y el post-procesamiento donde se visualizan y anali- zan los resultados encontrados en la segunda etapa de la implementación del MEFG.

Solution of the Rayleigh-Plesset Equation Through the Finite Element Method

Abstract

In this work we present numerical solutions of the Rayleigh-Plesset equation which describes the evolution of cavitating bubbles. In order to do that, we consider FEMG (Finite Element Method Galerkin); this simulation is performed for an inviscid and incompressible fluid in an uniform temperature field with constant surface tension, and the cavitation model into the which the pressure inside bubbles is equal to the fluid vapor pressure. Thus, in this problem is considered the Dirichlet boundary problem, and we obtained criteria for the boundary conditions at the cavitation phenomenon through to the which give rise to the bubble growing.

Key words: cavitation, Rayleigh-Plesset equation, Galerkin's finite element method, vapor pressure, inviscid fluid, incompressible fluid.

(ProQuest: ... denotes formulae omitted.)

1 Introducción

El proceso de formación de burbujas o cavidades en un fluido como con- secuencia del decrecimiento en la presión a una temperatura constante es conocido como cavitación[1], La dinámica de tales burbujas es de in- terés principal en el área de la hidráulica, dado la presencia de éstas puede ocasionar una disminución en la eficiencia y potencia de bombas e hidroturbinas; adicional a ruidos, vibración y erosión en materiales de dichas turbomáquinas [2],[3],[4],[5]. Es así, como el estudio de éste fe- nómeno, comenzado por Euler hacia 1754, y consolidado por Rayleigh hacia 1917 y Plesset en la década de los cincuenta, ha sido una rama activa de la mecánica de fluidos, lo cual ha generado rápidos desarrollos en esta área, más aun por el progreso simultáneo de las computadoras [6],[7],[8],[9] y de métodos experimentales [10],[11],[12], Sin embargo, da- da la complejidad del fenómeno, todavía no es posible tener del todo una teoría bien consolidada que describa en gran detalle cada aspecto de la cavitación, principalmente por la cantidad de fenómenos involucrados como cambios de fase, cambios abruptos en las variables termodinámi- cas, condiciones que favorecen o desfavorecen el desarrollo de cavidades, entre otros. En este trabajo, se solucionó numéricamente la ecuación de Rayleigh-Plesset, la cual describe la dinámica de las burbujas involu- cradas en el fenómeno de cavitación, por el método del elemento finito de Galerkin (MEFG). Tal solución se realizó, por medio de una apro- ximación lineal para el radio de la burbuja y usando una discretization del dominio dado en 59, 299 y 599 elementos así como considerando diferentes condiciones de Dirichlet para el problema estudiado.

2 Aspectos teóricos

En un fluido, la transición de fase de líquido a gas puede ocurrir al presentarse dos fenómenos. El primero de ellos es el de "ebullición ", el cual se presenta debido a que se aumenta la temperatura en el sistema a una presión constante como se muestra en la Figura 1.

Los valores de temperatura y presión a la cual ocurre la ebullición son conocidos como temperatura y presión de saturación. Ahora bien, si se sigue la flecha de izquierda a derecha a una presión constante (presión de saturación) el líquido se evaporará completamente. En el mismo sentido, en el fluido se puede presentar la aparición de la fase de vapor pero ya no como una consecuencia de la ebullición sino por el fenómeno conocido como cavitación, en el cual ocurre que ahora en el fluido el valor de la presión en algunas partes o regiones del fluido cae por debajo del valor de la presión de vapor de dicho fluido y se evidencian en él pequeños núcleos o espacios vacíos de gas, llamados burbujas o cavidades. Como puede ahora observarse en la Figura 2, si se decrece en el valor de la presión a una temperatura constante el líquido pasa a ser vapor.

Una vez se presentan las cavidades en el fluido, se da un proceso dinámico, el cual consta de crecimiento, decrecimiento y colapso. Esta dinámica, viene descrita por la ecuación de Rayleih-Plesset(ERP) [13] dada por

... (1)

Siendo R(r,t) el radio de la burbuja, R la derivada temporal de la burbuja , pb la presión al interior de la burbuja, Poo la presión del fluido que rodea a la burbuja, p la viscosidad, a la tensión superficial y p la densidad del fluido respectivamente. Como el estudio se realiza en un fluido invíscido, se tiene que p = 0 y (1) se reduce a:

... (2)

Donde (2) es la ecuación diferencial que se solucionará en la modela- ción numérica.

3 Modelación numérica

En general, la implementación del método del elemento finito cuenta con tres etapas: el pre-procesamiento, la solución de ecuaciones y el post- procesamiento. En la primera de ellas, que se describe a continuación, trata básicamente con el modelo matemático del problema, la geometría, la generación de la malla y la aplicación de las condiciones de Dirichlet y/o Neumann dadas. En esta misma sección, se describe la etapa de la solución de las ecuaciones para el problema y en la siguiente sección es dada la etapa del post-procesamiento, la cual básicamente consiste en visualizar y analizar los resultados encontrados.

Para la solución numérica del problema por medio del MEFG se utilizan los datos dados por la Tabla 1. El proceso se realiza usando elementos lineales y el enmallado se genera a lo largo de la dirección radial de crecimiento de la burbuja en el flujo como se observa en la Figura 3. La implementación del MEFG se hace al escogerse como región del enmallado el mostrado en la Figura 3 al discretizarse la región con 60 nodos y 59 elementos, entre 0,0001 m y 0,006 m, tomando como longitud de cada elemento (/q =0.0001m); para la región discretizada con 299 elementos y 300 nodos entre 0,0001 m y 0,03 m y para la región discretizada con 599 elementos y 600 nodos 0,0001 m y 0,06 m con la misma longitud del elemento de la región con 59 elementos.

Para dar solución a la ecuación dada por (2), se realiza la sustitución u = clR/clt, lo cual conduce a

... (3)

Así, se organiza (3) para expresar el residuo R(A') de la ecuación diferencial, el cual está dado por

... (4)

En el estudio llevado a cabo, se realiza una aproximación lineal para R, de tal forma que la solución en cada elemento puede ser escrita como una combinación lineal de funciones denominadas funciones de prueba [14], las cuales son usualmente polinomios. De esta manera, la aproxi- mación para la solución viene dada por:

... (5)

Siendo N el número de nodos. Las funciones de prueba se encuentran por medio de las funciones de interpolación de Lagrange [15], y están dadas por

... (6)

Donde N0 y Nr son las funciones de prueba para cada nodo en cada elemento y hi = Rr - Ro es la longitud de cada elemento. De esta manera, se utiliza el método del elemento finito de Galerkin (MEFG), el cual multiplica el residuo por las funciones de prueba [16]. Una vez se realiza esto, se efectúa la integral sobre el elemento y se iguala a cero; es decir, la forma ponderada

... (7)

Al usar (4) en (7), esta última se convierte en

... (8)

... (9)

Las ecuaciones dadas por (8) y (9) aquí dadas para las funciones Nq(R) y Nr(R), tienen que ser integradas entre la longitud de cada elemento al usar las funciones de prueba dadas por (6). Lo anterior conlleva a las respectivas ecuaciones de los elementos, las cuales dan como incógnitas los valores de u0 y uR en los nodos, y vienen dadas por:

... (10)

El siguiente paso en el proceso de solución una vez se ha llegado a 40), consiste en hacer la combinación o el así llamado "ensamble" de las ecuaciones de los elementos para construir un conjunto de ecuaciones algebraicas simultáneas. Dicho sistema, viene dado de forma general por

... (11)

Donde [K] representa la matriz de rigidez global y sus elementos están dados por la parte encerrada entre corchetes cuadrados de (10). Así mismo, {u} es el vector que representa las incógnitas del problema o los valores de las incógnitas en los nodos, y el vector {/} representa el vector fuente dado por el lado derecho de (10).

En la construcción de la matriz de rigidez global, se hace uso de la tabla de conectividad del elemento [17] dada por la Tabla 2, para encontrar las matrices de rigidez local y así la matriz global de rigidez estar dada por [18]

... (12)

En (12), N representa el número de nodos y e el número de elementos.

Una vez se ha realizado dicho ensamble, se tiene como anteriormente se dijo, un sistema algebraico cuya representación general viene dada por (11). Así, dicho sistema se soluciona por medio del proceso iterativo de Gauss-Seidel [19] obteniéndose los valores para u, los cuales son utili- zados para realizar la integración numérica de u = dR/dt por medio de la regla trapezoidal [20], para finalmente dar solución a la ERP.

4 Evolución de las burbujas

El crecimiento de una burbuja típica en el fluido fue determinado a tra- vés de la solución de la ecuación de Rayleigh-Plesset mediante el MEFG. De modo que, los valores obtenidos para la evolución del tamaño de la burbuja son presentados en la Figura 4, para lo cual la región (burbuja) fue discretizada en 59 elementos con condiciones de Dirichlet: presión en el flujo de 1000 Pa, presión al interior de la burbuja de 2340 Pa, es decir, es utilizado en el flujo el modelo de cavitación pb = pv, donde pv es la presión de vapor del fluido, y se supone que dicha presión no varía una vez se inicia el crecimiento de la cavidad en el fluido cuya temperatura es de 293 K. El radio máximo alcanzado por la burbuja es del orden de 20 veces su radio inicial, es decir, 0,002 m. Una vez es alcanzado dicho valor, se inicia el decrecimiento del cual se tiene información solamente hacia los 150 ms debido a las oscilaciones que se presentan en el radio. Como se puede observar, dichas oscilaciones se podrían asociar desde un punto de vista dinámico a colisiones de la burbuja en el flujo con otras burbujas o pequeñas partículas presentes en el fluido. De otra parte, las oscilaciones también pueden estar asociadas a cambios bruscos en las variables termodinámicas de la burbuja, como por ejemplo que la supo- sición inicial del modelo de cavitación no sea del todo válida, puesto que como la cavidad está creciendo continuamente la presión varíe debido a que el volumen de la cavidad no permanece invariante en el tiempo. Así mismo, es posible que dichos cambios termodinámicos abruptos se asocien al cambio de fase que se vuelve a presentar nuevamente en el flujo, en este estadio de gas a líquido; y en donde la variación de can- tidades como la presión, entalpia, energía, entropía a nivel local y en la interfase de la burbuja presenta saltos, por lo que se puedan llegar a tener, por ejemplo, procesos ele carácter irreversibles en el fenómeno de la cavitación, los cuales pudiesen llegar a ser descritos en principio por la teoría de las fluctuaciones. Sin embargo, es importante mencionar que dicho cambio de fase es de primer orden y entonces su estudio se diese por la ecuación de Clausius-Clapeyron, lo cual evidencia la complejidad para poder describir el fenómeno estudiado debido a las razones de pro- ducción de entropía [21] ya la influencia de la caída de presión en la producción de entropía a nivel local [22],

Así mismo,en la Figura 4 se ha mostrado el tiempo en el cual se da el colapso de la cavidad o burbuja en el flujo y se encuentra que dichos tiempos de vida están alrededor de 150 ms para la discretization realizada con una presión constante de 1000 Pa en el fluido y del orden de los 30 ms cuando la presión está representada por campo oscilante sinusoidal, 1000 * sen[w * t] donde la frecuencia de oscilación (w) fue de 1 Hz ; dicho campo permite evidenciar el efecto de un campo de presión dependiente del tiempo en el crecimiento de la burbuja, como en el caso de la cavitación acústica o vibratoria [23]. Así, la dinámica de las cavidades será diferente, lo cual sería un indicio de la influencia de la presión sobre la frontera de la burbuja.

En la Figura 5, se representa la evolución del radio de la burbuja o cavidad para la región discretizada con 299 elementos. El colapso en ambos casos se presenta en tiempos cercanos a los 150 ms cuando las condiciones para la presión en el fluido son de 1000 Pa y 1000* sen[w*t\ para una frecuencia ele oscilación de 1 Hz. De igual manera, el salto en el radio se da en un tiempo no mayor a 10 ms, a diferencia del caso discutido anteriormente y mostrado en la Figura 4 y cuyo orden estuvo alrededor de 27 veces el radio inicial, esto es, 0,0027 m. Esto puede ocurrir principalmente, por una mejora en la aproximación dado de que se esperaría una mayor precisión en los resultados una vez se ha discretizado la región con más elementos. En un tiempo posterior al inicio del decrecimiento en la burbuja, el radio de ésta presenta oscilaciones aún más fuertes y que se evidencia al usar una mayor cantidad de elementos. De igual manera, dicho comportamiento oscilatorio se debe a los cambios bruscos en las variables termodinámicas como se concluyo para la Figura 4. La característica principal que se encuentra es el tamaño alcanzado por la burbuja en el flujo y el tiempo de vida de ésta.

Una forma similar para la evolución del radio es presentada en la Figura 6, en la que la región ha sido discretizada en 599 elementos. El colapso para la cavidad se presenta entre los 250 y 260 ms cuando la presión en el fluido fue de 1000 Pa y 1000 * se/? [ur * í] para una frecuencia de oscilación de 1 Hz. El salto en el radio es aún más evidente, el cual no es superior a los 8 ms y cuya longitud alcanzada fue de 29 veces el radio inicial, esto es 0,0029 m. La mejora en la aproximación hace que sea más evidente las fluctuaciones en las variables termodinámicas, y a diferencia de lo concluido para la Figura 4, principalmente dichas oscilaciones son debidas a procesos de carácter irreversible como flujos de entropía y cambios en la presión a nivel local más no a algunos de carácter dinámico, los cuales pueden presentarse, pero su contribución podría no ser tan significativa en comparación con los cambios experimentados a nivel local y saltos en la interfase de la burbuja de las variables termodinámicas. Se evidencia que los cambios en la dinámica de la burbuja, se presentan esencialmente al igual que para la Figura 4 y 5 en los tiempos de vida y en la máxima longitud alcanzada por el radio de las burbujas.

5 Conclusiones

Se determina la evolución del radio de la burbuja por medio del MEFG para una discretization de 59, 299 y 599 elementos para la región (bur- buja). Se encuentra que para las condiciones de Dirichlet en la presión del fluido, los tiempos de vida para la cavidad en el flujo son diferentes como se muestra en las Figuras 4, 5 y 6. Sin embargo, para los casos de la Figura 5 y 6, los tiempos de vida son similares lo cual puede estar relacionado a una mayor precisión en la solución al discretizar la re- gión con más elementos y que la frecuencia de oscilación para la presión dependiente del tiempo (1000 * sen[w * ¿]) fue pequeña, (1 Hz).

Para la Figura 4, se encuentra un radio máximo alcanzado por la burbuja en el flujo de 0,002 m, es decir, del orden de 20 veces el radio inicial. Para la Figura 5, el radio máximo alcanzado es de 0,0027 m, siendo en este caso del orden de 27 veces el radio inicial y 1,35 veces el radio que alcanzó la burbuja para el caso de la Figura 4. Así mismo, para la Figura 6 se encuentra un radio máximo de 0,0029 m, el cual es del orden de 1,07 y 1,45 veces el radio alcanzado por la burbuja para el caso de la Figura 5 y 4, respectivamente. Como era de esperarse, si la burbuja alcanza mayores radios, ésta tendrá tiempos de vida mayores como se evidencia en las Figuras 4, 5 y 6.

En todos los casos, se presentan oscilaciones para el radio de la bur- buja, lo cual puede estar asociado a las fluctuaciones y saltos en los valores las diferentes variables termodinámicas (energía, entropía, ental- pia, presión, etcétera), puesto que acá se da la transición de fase que se presenta en el fluido-de vapor a líquido-, lo cual puede hacer que se lleguen a presentar procesos de carácter irreversible a nivel local, por las razones de producción de entropía[21] y la variación a nivel local de la presión[22]. Así mismo, al discretizarse la región con más elementos, se encontró mayores oscilaciones para el radio de la burbuja para las Figuras 5 y 6 y por ende una pérdida más notoria en la información de la evolución del radio en el flujo como era de esperarse. La aproximación que se realizó para la solución en el elemento al usar MEFG fue lineal; sin embargo, los resultados predicen de buena manera el comportamien- to que se esperaría ocurriese en el fenómeno de la cavitación de acuerdo con el crecimiento, decrecimiento, oscilaciones que se presentan duran- te el tiempo de vida de las burbujas y colapso, en concordancia con lo mostrado por otros autores para la solución del problema haciendo uso de otros métodos numéricos [24], [25],[26].

Agradecimientos

Este trabajo se realizó en el marco del proyecto de investigación "Hi- drodinámica de Flujos con Vorticidad ", apoyado por el Centro de In- vestigaciones y Desarrollo Científico (CIDC) de la Universidad Distrital Francisco José de Caldas.

References

Referencias

[1] E. Brennen C, Cavitation and Bubble Dynamics. New York: Oxford University Press, 1995. 148

[2] C. Tropea and A. Yarin, Springer Handbook of Experimental Fluid Mechanics. Würzburg: Springer, 2007. 148

[3] N. Berchiche, J.-P. Franc, and J.-M. Michel, "A cavitation erosion model for ductile materials," J. Fluids Eng., vol. 124, no. 3, pp. 601-606, 2002. 148

[4] J. Hundemer and M. Abdel-Maksoud, "Prediction of tip vortex cavitation in- ception on marine propellers at an early design stage," in 7th International Symposium on Cavitation, Michigan-USA, 2009, pp. 870-876. 148

[5] A. Philliph and W. Lauterborn, "Cavitation erosion by single laser-produced bubbles," J. Fluid Mech., vol. 361, no. 1, pp. 75-116, 1998. 148

[6] T. Tao Xing, Z. Li, and S. Frankel, "Numerical Simulation of Vortex Cavitation in a Three-Dimensional Submerged Transitional Jet," J. Fluids Eng., vol. 124, no. 4, pp. 714-725, 2005. 149

[7] D. Fuster, G. Agbaglah, C. Josserand, S. Popinet, and S. Zaleski, "Numerical simulation of droplets, bubbles and waves: state of the art," Fluid Dyn. Res., vol. 41, p. 065001 (24pp), 2009. 149

[8] Y. A. O. Xiong-liang and Z. A-man, "A numerical investigation of bubble dy- namics based on the potential-flow theory," Journal of Marine Science and Ap- plication, vol. 5, no. 4, pp. 14-21, 2006. 149

[9] M. Ohta, D. Kikuchi, and M. S. Yutaka Yoshida, "Robust numerical analysis of the dynamic bubble formation process in a viscous liquid," Int. J. Multiph. Flow, vol. 37, no. 9, pp. 1059-1071, 2011. 149

[10] X. Liu and J. Katz, "Cavitation phenomena occurring due to interaction of shear layer vortices with the trailing corner of a two-dimensional open cavity," Phys. Fluids, vol. 20, no. 4, p. 041702, 2008. 149

[11] G. F. Oweis, J. Choi, and S. L. Ceccio, "Dynamics and noise emission of laser induced cavitation bubbles in a vortical flow field," J. Acoust. Soc. Am., vol. 115, no. 3, pp. 1049-1058, 2004. 149

[12] M. T. S. Tabar and Z. Poursharifi, "An Experimental Study of Tip Vor- tex Cavitation Inception in an Axial Flow Pump," World Academy of Scien- ce, Engineering and Technology, vol. 73, no. 1, pp. 527-530, 2011. 149

[13] D. Joseph, T. Fuñada, and J. Wang, Potential Flows of Viscous and Viscoelastic Fluids. New York: Cambridge University Press, 2007. 150

[14] Y. X-S., A First Course in Finite Element Analysis. United Kingdom: Luviner Press, 2007. 153

[15] T. Chung, Computational Fluid Dynamics. United Kingdom: Cambridge Uni- versity Press, 2002. 153

[16] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method. New York: Dover, 2009. 153

[17] T. Chandrupatla and A. Belegundu, Introduction to Finite Elements in Engi- neering, 2nd ed. U.S.A.: Prentice Hal, 1998. 155

[18] D. Logan, A First Course in the Finite Element Method. U.S.A.: PWS Publis- hing Company, 1992. 155

[19] D. Young, Iterative Solution of Large Linear Systems. New Yor: Dover, 2003. 156

[20] P. Davis and P. Rabinowitz, Methods of Numerical Integration, 2nd ed. New York: Dover, 2007. 156

[21] Z. Bilicki, M. Giot, and R. Kwidzinski, "Fundamentals of two-phase flow by the method of irreversible thermodynamics," Int. J. Multiph. Flow, vol. 28, no. 12, pp. 1983-2005, 2002. 157, 160

[22] R. Revellin, S. Sthépane, Sammer, and J. Bonjour., "Local entropy generation for saturated two-phase flow," Energy, vol. 34, no. 9, pp. 113-1121, 2009. 157, 160

[23] K. S. Suslick and D. J. Flannigan, "Inside a Collapsing Bubble: Sonolumines- cence and the Conditions During Cavitation," Annu. Rev. Phys. Chem., vol. 59, pp. 659-683, 2008. 157

[24] H. Alehossein and Z. Qin, "Numerical analysis of Rayleigh-Plesset equation for cavitating water jets," Int. J. Numer. Meth. Engng., vol. 72, no. 7, pp. 780-807, 2007. 160

[25] C. Crowe, Multiphase Flow Handbook. USA: CRC Press, 2006. 160

[26] J.-P. Franc and J.-M. Michel, Fundamentals of cavitation. USA: Kluwer Aca- demic Publishers, 2004. 160