Análisis de la influencia del número de venecianas en intercambiadores de calor compactos

Analysis of the influence of louver number in compact heat exchangers

Alberto Menéndez-Pérez^I, Rubén Borrajo-Pérez^I, Daniel Sacasas-Suarez^II

^I Universidad Tecnológica de La Habana. Centro de Estudio de Tecnologías Energéticas Renovables, CETER. La Habana, Cuba
^II Centro de Inmunoensayo, CIE. La Habana. Cuba

RESUMEN

El trabajo presenta un estudio numérico de intercambiadores de calor compactos tubo-aleta venecianas. Se determina la influencia del número de venecianas sobre las características termo hidráulicas del intercambiador. Se utilizó un modelo numérico en 3D dentro del régimen laminar y para una variación del número de Reynolds, basado en el diámetro hidráulico, entre 120 y 1200. Variando la cantidad de venecianas aguas arriba y aguas abajo de la veneciana central, se determina que, para menores números de Reynolds las aletas con 2 venecianas, se obtiene el  mayor coeficiente de transferencia de calor. Por otra parte, al aumentar la velocidad del aire a la entrada del modelo, se obtienen mejores resultados de transferencia de calor para las geometrías con mayor número de venecianas.

Palabras claves: aletas venecianas, intensificación de transferencia de calor, simulación numérica.

ABSTRACT

The paper presents a numerical study of compact louvered fin -and- tube heat exchangers. The influence of the louver number in the thermo-hydraulic characteristics of the heat exchanger is determined. A 3D numerical model within laminar flow conditions, for a Reynolds number variation spanning from 120 to 1200, based on the hydraulic diameter. Changing the louver number up and downstream from the central louver, it was found that for smaller Reynolds number, the 2 louver configuration had the highest heat transfer value. On the other hand, as the inlet velocity at the model increases, higher values of heat transfer coefficient appear for geometries having more louvers.

Key words: louvered fins, heat transfer enhancement, numerical simulation.

INTRODUCCIÓN

Entre las superficies extendidas, las más usadas son las aletas planas, venecianas, alternadas y las onduladas así como las posibles combinaciones de estas. Los generadores de vórtices se han sumado a estas geometrías como elementos intensificadores de la transferencia de calor. El comportamiento de los flujos es un aspecto decisivo en la transferencia de calor, por lo que existen dos métodos que permiten alterar los patrones de este, según se plantea en T’Jo en et al [4]. Estos métodos se definen como: el método de  alteración del flujo principal y la introducción de flujos secundarios. Puede ocurrir el caso en que ambos mecanismos estén presentes. La diferencia radica en que en el primer método; las características principales son manipuladas a través de cambios geométricos, mientras que en el segundo se modifica la estructura del flujo local.

Los modelos de aletas venecianas son una de las superficies extendidas intensificadas más avanzadas y se forman esencialmente cortando la placa de metal de la aleta a intervalos y girando las tiras de metal formando un ángulo que re-direcciona el flujo [5].

Existen combinaciones de intercambiadores de calor con aletas venecianas y tubos elípticos, planos o circulares. Las geometrías con tubos planos posee mejor eficiencia de las aletas que las geometrías con tubos circulares [3]. Las aletas venecianas provocan generalmente altos valores de transferencia de calor por unidad de área, pero en cambio generan grandes caídas de presión. La razón por la que las venecianas constituyen una técnica compleja de analizar se debe a que poseen un diverso número de parámetros geométricos. Uno de los parámetros que distinguen a las venecianas es el ángulo de inclinación de las venecianas (louvers), por esta razón, muchos autores se han dedicado a investigar y establecer los rangos óptimos de estos ángulos.
[1-3]

En intercambiadores que presentan tubos circulares o elípticos, se generan vórtices longitudinales cuando el flujo pasa alrededor del tubo y estos son llamados vórtices de herradura. La principal característica de estos es que generan altos valores de transferencia de calor en la zona delantera de los tubos y alrededor de los mismos [6, 7]. Sin embargo, la transferencia de calor detrás del tubo es pobre debido a la recirculación del flujo en la zona muerta o trasera del tubo.

Para números de Reynolds bajos, el espesor de la capa limite es tal que las venecianas quedan embebidas dentro de esta, forzando así al fluido a cambiar el canal por donde circula. Al aumentar el número de Reynolds, disminuye el flujo que cambia de canal y aumenta el que mantiene la dirección principal. Lo anterior motivado por la presencia de una capa limite más fina.

Karthik et al [8] concluyen, del estudio de un sistema de calentamiento de aire con agua usando aletas venecianas con tubos elípticos que, aumentar la velocidad de entrada del aire produce un mayor incremento del coeficiente de traspaso de calor global que el obtenido por aumentar el flujo másico de agua por los tubos.

Sanders [9] plantea que la mayoría de los estudios relacionados a este tema se realizan en dos dimensiones ignorándose así los efectos creados en la superficie de los tubos. Esto se debe a que recientes artículos han concluido que los análisis de flujos tridimensionales próximos a la pared de los tubos no influyen en el funcionamiento del equipo intercambiador.

Algunos autores concluyen que usando las aletas venecianas con diferentes grados de inclinación cada una, se pueden obtener buenos resultados para pequeños valores de Reynolds. El estudio consiste en comparaciones, variando el paso entre aletas y el ángulo de inclinación de las venecianas para un particular arreglo donde el número de Reynolds toma valores entre 100 y 1000 aproximadamente [1, 7, 8]. Un dato muy interesante es que para valores pequeños del paso entre aletas y grandes valores del ángulo de inclinación se logran grandes y fuertes desviaciones de flujo de aire.

Según Ameel et al [11], el comportamiento de los intercambiadores de calor como caja negra es analizado comúnmente mediante el método de la Logarithmic Mean Temperature Difference (LMTD) o mediante el método de las NTU (Number of Transfer Units).

Existen investigaciones donde se demuestra que las venecianas re-direccionan el flujo de aire en la dirección paralela a sus propios planos [12,13], donde se utiliza una visualización del flujo en modelos escalados para llegar a esa conclusión.

Achaichia y Cowell [14] concluyen en una investigación que los flujos no deseados entre aletas son provocados cuando el paso entre venecianas es combinado con grandes espaciamientos entre aletas y pequeños ángulos de las venecianas. El anterior criterio es también soportado por Xiaoping et al [15] quienes afirman que se obtienen mejores resultados usando geometrías con pequeños pasos entre venecianas (1,6 mm) y pequeños ángulos de las mismas, 26º.

Kang et al [16] arriban a la conclusión de que el empleo de arreglos de venecianas y múltiples tubos circulares produce entre un 10 % y un 15 % más de desempeño de traspaso de calor que si se emplean tubos planos para los mismos flujos volumétricos y velocidades.

Aunque la revisión bibliográfica muestra la existencia de múltiples estudios sobre este tipo de aletas, existe poca información de la influencia de la cantidad de venecianas sobre las características termo-hidráulicas de estos intercambiadores de calor. Precisamente este trabajo determina la influencia del número de venecianas sobre el calor transferido y la caída de presión en un rango de Reynolds entre 120 y 1200.

MATERIALES Y MÉTODOS

En la figura 1 se muestra una vista frontal de una aleta veneciana donde se aprecian, entre otras, el ángulo de inclinación de las venecianas y la altura de las mismas. La figura 1 corresponde con una aleta que posee un número de venecianas igual a 2 a ambos lados de la veneciana central, comúnmente llamada veneciana redireccionadora.

Las ecuaciones que rigen el estudio de un fenómeno donde interesa determinar el movimiento de un fluido y el intercambio de calor de este con las superficies que este contacta son la ecuación de continuidad, las ecuaciones de cantidad de movimiento en cada uno de los ejes y la ecuación de la energía. Estas ecuaciones en el dominio computacional, para flujo incompresible, con propiedades constantes y en estado estacionario sin disipación viscosa y en régimen laminar, pueden ser expresadas en el mismo orden en que se mencionan, como sigue, ecuaciones 1, 2 y 3:

La solución de las ecuaciones 1, 2 y 3 se logra para un dominio computacional que se muestra en la figura 2. El dominio, como se observa, ha sido extendido en la región de entrada y la de salida del modelo. Fue extendido 7 veces la separación entre aletas en la dirección de entrada y 7 veces el diámetro menor del tubo en la dirección de salida. La necesidad de tener un perfil de velocidad uniforme y unidimensional en la entrada del modelo, así como evitar la existencia de flujo reverso en la sección de salida es el motivo de estas extensiones.

La interpolación utilizada para la energía fue de segundo orden del tipo upwind pues esta es más efectiva cuando no existe una malla totalmente normal al flujo. En el acoplamiento presión-velocidad fue empleado el método simple y un esquema de primer orden upwind en las ecuaciones de momento y el modelo estándar en la presión. La aleta en la región central del dominio se toma como un sólido, mientras que los canales sobre y bajo la misma son considerados como regiones fluidas. Dos condiciones de contorno que se necesitan son, la temperatura en el tubo y los parámetros de entrada del fluido al canal. La temperatura de los tubos se considerará igual a la del fluido refrigerante que circula por dentro de estos, por lo que se desprecia el valor del coeficiente de traspaso de calor por la parte interior del tubo dado su elevado valor cuando se compara al que debe existir por la parte externa del tubo. La velocidad de entrada es un dato que puede variar pues el estudio de diferentes números de Reynolds precisa del cambio de la velocidad en la entrada. Las condiciones de contorno en general por regiones se resumen a continuación:

A la entrada del modelo, ecuación 4:

En la superficie del tubo tendremos, ecuación 6:

En las regiones laterales del dominio y consideradas como simétricas se cumple que (ecuación 7):

A la salida del modelo, ecuación 8:

El modelo trabaja bajo régimen de flujo laminar con bajas velocidades del aire en el canal, en previsión de la existencia de molestos ruidos cinéticos. Se modela además en estado estacionario. El flujo a analizar es tridimensional, con los campos de velocidad y de temperatura desacoplados, lo que garantiza la independencia entre ambas variables. La velocidad del flujo a la entrada del modelo se varía entre 0,5 y 5 m/s. La temperatura del aire a la entrada del modelo es 300 K y en la pared del tubo se considera igual a 286 K. El flujo es tratado como incompresible.

La región de la geometría y el canal son malladas con elementos de tipo Tetraedros híbridos (TGrid) en el software Gambit 2.4.6. Esta es mallada de forma que la densidad de elementos fuera capaz de producir una solución independiente de la malla, ya que la diferencia entre los resultados de caída de presión y coeficiente de transferencia de calor no excede el 1,6 % entre dos tamaños de malla sucesivos. Esta no posee ningún elemento con volúmenes de acrecentada esbeltez (skewness), ni ningún elemento con volumen negativo o invertido, pues ambos conspirarían contra la convergencia de las iteraciones. Para el mallado se utilizaron funciones de talla en las caras más complejas, estas se muestran en la figura 3.

Validación del modelo numérico

Con el objetivo de considerar válidos los resultados se debe certificar el método de obtención utilizado para su determinación. El método empleado aquí en la certificación consiste en construir un modelo computacional con dimensiones y características correspondientes con uno que ya había sido investigado y sus resultados están disponibles en la literatura. Se seleccionó el modelo estudiado por Han et al [17], que a la vez fue validado contra resultados experimentales obtenido por Wang et al [18, 19]. Han et al y afirman que las desviaciones promedio del factor de fricción (f) y del factor de Colburn (j) entre su modelo numérico y su validación en base a los experimentos de Wang et al son de un 9,3 % y un 4,5 % respectivamente. Los experimentos de los autores se basaron en una geometría similar de aleta veneciana con tubo elíptico pero usando las venecianas con longitudes variables alrededor del borde del tubo. Otra diferencia es que el aire entrante al túnel de viento posee una temperatura de 308 K y la temperatura de la pared de los tubos es de 353 K. Dicho estudio tiene todas las dimensiones distintas al modelo estudiado en este trabajo. Se usan entonces los mismos tipos de flujos, interpolaciones, condiciones de contornos y criterios de convergencia para la energía y el momentum. También se implementa un acoplamiento entre la velocidad y la presión mediante el algoritmo SIMPLEC.

En la figura 4 se observa que para los valores de Δp el comportamiento es bastante aproximado, aunque el modelo numérico subestima este parámetro ligeramente. El comportamiento del coeficiente de transferencia de calor es también aproximado por el modelo numérico. La diferencia entre los comportamientos de estas curvas pueden ser ocasionadas por la inexactitud en la reproducción de la geometría mostrada en el trabajo de Han et al, la cual no ofrece los detalles necesarios para la exacta reproducción del modelo ni la forma constructiva correcta. Este trabajo sobrestima el valor de h para velocidades superiores a 1,5 m/s, mientras que subestima h para velocidades inferiores a 1,5 m/s si se compara con los de Han et al [17].

El promedio de las desviaciones entre los resultados para el coeficiente pelicular de traspaso de calor global y la caída de presión del trabajo de Han et al [17], son del orden de 5,22 % y 21,61 % respectivamente. Luego, a pesar de las diferencias antes mencionadas, puede notarse como la mayoría de los valores quedan abarcados dentro de la región de incertidumbre, que delimita las barras de error en la figura 4. En correspondencia con lo anterior podemos afirmar que es válido el método utilizado en el trabajo para lograr la simulación. Puede establecerse entonces que el método empleado, las simplificaciones, así como las condiciones de contorno establecidas en este trabajo son válidos para alcanzar los objetivos propuestos.

Reducción de datos

El número de Reynolds se determinó usando la velocidad en la sección mínima del canal u_min mientras la longitud característica es el diámetro hidráulico del canal D_h., ecuación 9:

El calor transferido en la superficie de intercambio puede ser calculado con el cambio en temperatura del aire entre las secciones de entrada y salida del modelo (T_out - T_in), el flujo de masa m_a, y el calor específico del aire C_pa, según la ecuación 10:

Cuando el fluido experimenta un cambio de fase por el interior de tubos, es práctica común considerar un elevado valor para el coeficiente pelicular de transferencia de calor interno. Al asumir esto el coeficiente global será obtenido solo considerando el coeficiente de transferencia por la parte externa y la conducción dentro de la pared del tubo. La temperatura de la pared interna de los tubos se considera constante y con el mismo valor que la del refrigerante que circula. El calor transferido puede también ser calculado a través de la conocida ecuación que involucra el logarítmico (LMTD), conociendo el área de transferencia A_f, el coeficiente global de transferencia , y la mencionada diferencia de temperaturas, ecuación 11:

El factor de corrección del LMTD, F, fue considerado nulo a causa de que uno de los fluidos mantiene su temperatura constante. EL coeficiente global de transferencia de calor es calculado considerando la igualdad de los calores expresados por las ecuaciones 10 y 11. Q = Q_h

La eficacia de la aleta η₀ está involucrada y es a su vez una función del coeficiente de transferencia de calor. La eficacia se calcula en función de la eficiencia de la aleta η y el área total de transferencia A₀, ecuación 12:

La eficiencia de la aleta para una aleta rectangular se determina usando el método aproximado desarrollado por Schmidt para aletas circulares. La eficiencia de la aleta se expresa según, ecuación 13:

El valor de m es calculado con la conductividad térmica de la aleta (k_f) y su espesor (f_t), ecuación 14

Donde X_m es la mitad del espaciamiento transversal y X_L es calculado según, ecuación 17:

La eficacia de la aleta y el coeficiente global de transferencia de calor tienen una formulación implícita por lo tanto se necesita de un proceso iterativo. La igualdad de calores de las ecuaciones 10 y 11 es la que debe ser satisfecha. Para una geometría fija existe solo un par de valores de estas magnitudes, eficiencia de aleta y coeficiente global de transferencia, que cumplen esta condición. Luego el factor de Colburn puede obtenerse a partir del número de Prandtl y el factor de fricción según, ecuación 18:

Donde A_min es el area mínima de la sección de paso y A_t es el area de transferencia de calor.

RESULTADOS y DISCUSIÓN

Es costumbre, como fue mencionado antes, en la literatura actual presentar el comportamiento de los intercambiadores de calor utilizando números adimensionales y en consecuencia se muestran, en las figuras 5 y 6, los valores del Factor de Colburn j y el Factor de fricción f como función del Número de Reynolds. Las leyendas de las figuras 5 y 6 se refieren, a la cantidad de venecianas a cada lado de la veneciana central, que produce un cambio de dirección en el ángulo de las mismas. Luego un numero 3 en la escala significaría un total de 6 venecianas adicionales a la central.

Puede notarse como para pequeños números de Reynolds las aletas con cantidad de venecianas (C_L) igual a 2 presentan los mejores resultados. Mientras, cuando crece el número de Reynolds no existen diferencias significativas en el comportamiento de los diferentes modelos. Sin embargo, en la figura 7 cuando se muestra el coeficiente de transferencia de calor en función de la velocidad de entrada, se mantiene la tendencia de, a menor cantidad de venecianas, mejores resultados a velocidades menores. Por otro lado, al aumentar la cantidad de venecianas, el comportamiento se invierte y comienzan a mostrar mejores resultados los modelos con esta característica.

Una probable explicación para este fenómeno será que cuando aumenta el número de venecianas, también lo hace la cantidad de superficies donde se desarrollan capas límites nuevas (bordes frontales de la veneciana). Cuando el número de venecianas es menor, estas consiguen acercar el fluido frio que circula por el centro del canal, a la superficie de la aleta. Lo anterior es posible pues, a menores números de venecianas implica venecianas mayores en longitud. 

En la figura 8 se pone de manifiesto que el mecanismo de transferencia de calor más importante, cuando las velocidades son mínimas, es el desarrollo de la capa limite en el frente de la aleta y en las venecianas, pues se puede observar la forma en que las venecianas quiebran la capa límite del fluido, lo cual influye en el incremento del calor intercambiado. Además se observa claramente la región de bajos coeficientes de transferencia de calor asociados a la zona muerta de los tubos. Las diferencias de calor intercambiado en la región de las venecianas se deben a la incidencia del flujo sobre las superficies

La figura 9 muestra las líneas de corriente generadas en la cara de entrada al canal para una extensión de la condición de simetría del modelo. Se observa claramente las diferencias de los vórtices de herraduras en el primer y segundo tubo. En el segundo tubo, estos vórtices tienden a contornear menos el tubo, provocando una región de recirculación mayor que la que se observa tras los tubos de la primera fila. Se observa además cómo se manifiesta el efecto de aceleración que experimenta el flujo al pasar por la segunda fila, disminuyendo la recirculación en la zona trasera del tubo.

Así mismo, la figura 10 muestra los perfiles de temperatura en la parte superior de la aleta. La dirección principal del flujo es de derecha a izquierda. Puede observarse que en la medida que el fluido posee una mayor velocidad a la entrada del modelo, la temperatura media de la aleta es más elevada y además existe menos uniformidad en la salida del fluido. Al ser esta superficie un intercambiador de calor donde el aire está a mayor temperatura que la superficie de la aleta, existe una mayor transferencia de calor, pues la diferencia de temperatura entre la superficie de la aleta y el fluido es mayor.  Se observa en la salida de la aleta que la temperatura nunca es constante en toda la superficie, destacando que a mayor velocidad, la temperatura posee mayores gradientes.

Limitaciones del trabajo

Los resultados del trabajo son válidos solamente para aletas de similar geometría y en el entorno de las dimensiones de la aquí estudiada, siempre y cuando el régimen de flujo sea laminar.

CONCLUSIONES

Se creó un modelo computacional capaz de reproducir el comportamiento termo hidráulico, mostrado en la literatura consultada, de una superficie de intercambio de calor de aletas venecianas y tubos circulares. Esta superficie fue utilizada como vehículo para la certificación del método.

Entre los resultados fundamentales se determinó, que para el modelo de dos filas de tubos, el flujo tiende cada vez más a separarse de la cara lateral a medida que aumenta la velocidad, mientras que en el primer tubo el comportamiento es el contrario, debido a la presencia del segundo.

El coeficiente de transferencia de calor se encontró aumentando cuando el número de venecianas es menor, resultado que es acentuado para los menores números de Reynolds.  Cuando se estudiaron las velocidades mayores, fueron los modelos con el mayor número de venecianas los que presentaron un mejor comportamiento termohidraúlico.

REFERENCIAS

1. Wang CC, Kuan-Yu C, Jane-Sunn L. An experimental study of the air-side performance of fin-and-tube heat exchangers having plain, louver, and semi-dimple vortex generator configuration. International Journal of Heat and Mass Transfer. 2015;80:281-7.

2. Wang CC, Kuan-Yu C, Yur-Tsai L. Investigation of the semi-dimple vortex generator applicable to fin-and-tube heat exchangers. International Journal of Heat and Mass Transfer. 2015;88:192-7.

3. Boris D, Anica T, Lenic K. Numerical investigation of heat transfer enhancement in a fin and tube heat exchanger using vortex generators. International Journal of Heat and Mass Transfer. 2014;78:662-9.

4. T’Joen C, Huisseune H, De Paepe M. Interaction between mean flow and thermo-hydraulic behavior in inclined louvered fins. International Journal of Heat and Mass Transfer. 2010;54:826-37.

5. Powar VS, Mirza MM. Performance of Louver Fin Pattern as Extended Surface Used To Enhance Heat Transfer - A Review. Int Journal of EngineeringResearch and Applications. 2013;3:1409-13.

6. Zhong Y, Jacobi AM. Experimental study of louver-fin flat-tube heat exchanger performance under frosting conditions. In: International Conference on Enhanced; Hoboken, New Jersey, USA; 2005. Citado octubre de 2015] Disponible en: http://dc.engconfintl.org/heatexchangerfall2005/29.

7. Salviano LO, Dezan DJ, Jurandir IY. Optimization of winglet-type vortex generator positions and angles in plate-fin compact heat exchanger: Response Surface Methodology and Direct Optimization. International Journal of Heat and Mass Transfer. 2014;82:373-87.

8. Pooranachandrana K, Ali Khan SI, Ramalingamd V. Experimental and numerical investigation of a louvered fin and elliptical tube compact heat exchanger. Thermal Science. 2015;19(2):679-92.

9. Sanders PA. Effects of Louver Lengthand Vortex Generators to Augment Tube Wall Heat Transfer in Louvered Fin Heat Exchangers. Virginia, Estados Unidos: Virginia Polytechnic Institute and State University; 2005. [Citado: diciembre 2015] Disponible en: http://hdl.handle.net/10919/35189.

10. Huisseune H, T’Joen C, De Paepe M. Influence of the louver and delta winglet geometry on the thermal hydraulic performance of a compound heat exchanger. International Journal of Heat and Mass Transfer. 2012;57:58-72.

11. Ameel B, Degroote J, De Paepe M. Interaction effects between parameters in a vortex generator and louvered fin compact heat exchanger. International Journal of Heat and Mass Transfer. 2014;77:247-56.

12. Schmidt TE. Heat transfer calculations for extended surfaces. Refrig Eng. 1949;12:351–7.

13. Ameel B, T’Joen C, De Paepe M. On fin efficiency in interrupted fin and tube heat exchangers. International Journal of Heat and Mass Transfer. 2013;60:557-66.

14. Achaichia A, Cowell TA. Heat Transfer and Pressure Drop Characteristics of Flat Tube and Louvered Plate Fin Surfaces. Experimental Thermal and  Fluid  Science. 1988;1:147-57.

15. Xiaoping T. CFD Simulation and Experimental Study on Air side Performance for MCHX. In: International Refrigeration and Air Conditioning Conference; Indiana, EEUU: Purdue University; 2010. [Citado Disponible en: http://docs.lib.purdue.edu/iracc/1023.

16. Kang H, Jacobi AM, Minkyoo L. Air-Side H eat Transfer Performance of Louver Fin and Multi-Tube H eat Exchanger for Fuel-Cell Cooling Application, in Air-Side H eat Transfer Performance of Louver Fin and Multi-Tube Heat Exchanger for Fuel-Cell. In: International Refrigeration and Air Conditioning Conference; Indiana, USA: Purdue University; 2012. [Citado octubre 2015] Disponible en: http://docs.lib.purdue.edu/iracc/1337.

17. Han H, Ya-Ling H, Ming W. A numerical study on compact enhanced fin-and-tube heat exchangers with oval and circular tube configurations. International Journal of Heat and Mass Transfer. 2013;65:686-95.

18. Wang CC, Lee CJ, Chang CT. Some aspects of plate fin-and-tubeheat exchangers: with and without louvers. J Enhanced Heat Transfer 6. 1999;6:357-68.

19. Wang CC, Tsai YM, Lu DC. Comprehensive study of convex-louver and wavy fin-and-tube heat exchangers. J Enhanced Heat Transfer 12, . 1998;12:423-30.

Recibido: 15 de abril de 2016. ]]> Aceptado: 21  de julio de 2016.

Alberto Menéndez-Pérez, Universidad Tecnológica de La Habana. Centro de Estudio de Tecnologías Energéticas Renovables, CETER. La Habana, Cuba
Correo electrónico: amenendez@mecanica.cujae.edu.cu