Derivación de la ecuación de Rayleigh-Plesset.

...

[pic 5]

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

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 postprocesamiento. 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 (hi =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.

[pic 6]

Figura 3: discretización del dominio y enmallado para la burbuja por medio de elementos unidimensionales

[pic 7]

Para dar solución a la ecuación dada por (2), se realiza la sustituciónu = dR/dt, lo cual conduce a[pic 8]

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

[pic 9]

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, las cuales son usualmente polinomios. De esta manera, la aproximación para la solución viene dada por:

[pic 10]

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

[pic 11]

Donde N0 y NR son las funciones de prueba para cada nodo en cada elemento y hi = RR − R0 es la longitud de cada elemento. De esta manera, se utiliza el método del elemento finito de Galerkin, el cual multiplica el residuo por las funciones de prueba. Una vez se realiza esto, se efectúa la integral sobre el elemento y se iguala a cero; es decir, la forma ponderada

[pic 12]Al usar (4) en (7), esta última se convierte en[pic 13]

Las ecuaciones dadas por (8) y (9) aquí dadas para las funciones N0(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:

[pic 14]

El siguiente paso en el proceso de solución una vez se ha llegado a (10), 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

[pic 15]

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 {f} 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 dada por la Tabla 2, para encontrar las matrices de rigidez local y así la matriz global de rigidez estar dada por

[pic 16]

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

[pic 17]Tabla 2: tabla de conectividad de los elementos. Los nodos 1 a N representan losnodos globales mientras que los nodos 1 y 2 son los nodos locales para cada elemento

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 utilizados 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.

Evolución de las burbujas

El crecimiento de una burbuja típica en el fluido fue determinado a travé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 P a, presión al interior de la burbuja de 2340 P a, 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 lasvariables termodinámicas de la burbuja, como por ejemplo que la suposición inicial