Metodología de diseño y dinámica computacional de fluido del impulsor y voluta de un compresor centrífugo

1. Introducción

En este artículo se muestra la combinación de las metodologías de Aungier [2], [3] y Lüdtke [4] para el diseño de impulsores centrífugos, con el fin de disminuir el uso de suposiciones e interpolaciones gráficas. Un trabajo de combinación de metodologías disponibles fue realizado por Meroni, et al. [5] demostrando la viabilidad del proceso al mantenerse una compatibilidad de las metodologías; sin embargo, su trabajo se enfocó en el diseño de un compresor centrífugo y no en la creación de una nueva metodología. Algunas metodologías modernas para el diseño de los impulsores centrífugos usan algoritmos de diseño genético por ejemplo Casey et al.[6], por su parte, Chu et al [7] presenta el uso de un algoritmo genético que predice el comportamiento de un compresor centrífugo. Análogamente, Hildebrand et al [8] examinan cual es el factor que más afecta el comportamiento de un compresor centrífugo, compilando los resultados en una base de datos con la cual se dan los lineamientos para el diseño de un impulsor. A pesar de los resultados satisfactorios obtenidos por los autores mencionados, el proceso de diseño del impulsor de un compresor centrífugo sigue llevando al diseñador a no tener contacto con toda la información disponible, limitando así el proceso de optimización. Por lo tanto, en el presente trabajo se crea una metodología de diseño para impulsores de compresores centrífugos la cual combina dos metodologías existentes y añade el principio fundamental de las turbomáquinas, con el fin delimitar el uso de ecuaciones empíricas, evitar las interpolaciones gráficas y dar al usuario la posibilidad de manipular las parámetros y ecuaciones con el fin de adaptar los diseños a situaciones particulares obteniendo beneficios con lo cual, se pueden obtener beneficios en dinámica de fluidos, termodinámica del proceso de compresión, y manufactura del impulsor.

Adicionalmente, se diseñó una voluta apropiada para el impulsor, siguiendo las indicaciones de Pfleiderer [9] quien predice las condiciones geométricas de la voluta al mantener un momento angular mayor a 1. El impulsor y la voluta diseñados se analizaron con simulación mediante modelo RANS (Reynolds Averaged Navier- Stokes) de manera que se puede analizar la viabilidad del diseño en fases preliminares.

Por último, se realizó la validación de las simulaciones al replicar los resultados Mojaddam et al [1] haciendo uso de la metodología de simulación expuesta en este artículo para determinar la desviación de los datos por medio de la comparación con el porcentaje de error.

2. Diseño del impulsor

El diseño sigue el diagrama de flujo presentado en la figura 1 en la cual se puede observar de forma general el procedimiento seguido por los autores para diseñar y simular el impulsor.

Las metodologías de diseño fueron elegidas porque presentan límites de diseño similares. Aungier [3] describe que su metodología puede ser seleccionada si el coeficiente de flujo ∅ se encuentra en el rango de 0.003 y 0.2 y una relación de presiones totales de 3.5. Por su parte, Lüdtke [4] limita el coeficiente de flujo desde 0.01 hasta 0.15 y la relación de presiones totales de extiende hasta 6. Interceptando los parámetros iniciales de ambos autores se encuentra el impulsor diseñado debe tener un coeficiente de flujo mayor a 0.01 y menor a 0.15, mientras que la relación de presiones totales debe ser menor a 3.5. Bajo estos parámetros se seleccionaron las condiciones iniciales, que incluyen, presión de entrada, temperatura de entrada, flujo másico, y relación de presiones totales, además de la velocidad angular.

En el Handbook de compresores centrífugos [10] se muestra que la velocidad angular del compresor adquiere valores entre 15000-20000 [RPM]; sin embargo, se afirma que estos valores pueden ser sobrepasados, como lo muestra Gutiérrez [11] en su recopilación para determinar las pérdidas, donde todos los compresores analizados superan el rango previamente descrito.

2.1. Propiedades termodinámicas

Determinar en una etapa temprana de diseño propiedades como eficiencia y cabeza isentrópicas es necesario, para lo cual se computaron las propiedades termodinámicas de salida de acuerdo con las necesidades, es decir, relación de presiones, temperatura de salida, flujo másico y el estado termodinámico de entrada. Este proceso determinó si se requería más de 1 etapa de compresión; ya que no es recomendable un trabajo mayor a 40[kJ/kg] por etapa [4].

Rolls-Royce [12], muestra que la presión estática puede ser aumentada como máximo en un 50% del total en el impulsor. Con lo cual, una de las propiedades termodinámicas de salida es determinada, usando la ecuación de Redlich-Kwong, modificada por Aungier [2], ecuaciones 1 a 5 se definió completamente el estado termodinámico. Es importante resaltar, que las ecuaciones de estado son función de 3 propiedades intensivas independientes, de las cuales, si no se da una restricción en la temperatura, se debe estimar teniendo en cuenta que se limita la energía que se puede brindar al fluido a 40[kJ/kg] por etapa.

Ecuación1

Ecuación 2

Ecuación 3

Ecuación 4

Ecuación 5

Dónde P denota presión, R la constante universal de los gases, T la temperatura, 𝑣 el volumen específico, 𝑇_𝑅 la temperatura reducida, 𝑇𝑐 la temperatura crítica del gas, 𝑝𝑐 la presión crítica del gas y 𝑝𝑣 representando la presión de vapor del gas a una temperatura reducida de 0.7

La cabeza y eficiencia politrópicas se determinaron por medio de las ecuaciones 6 y 9. Con el fin de corroborar que la restricción de cabeza isentrópica se mantiene y la eficiencia que arroja el proceso está dentro de los rangos aceptados por Lüdtke [4] y Aungier [2]; [3].

Ecuación 6

Ecuación 7

Ecuación 8

Ecuación 9

Donde,𝑇1 es la temperatura de entrada, 𝑃1la presión de entrada, 𝑇2 la temperatura de salida, 𝑃2 la presión de salida y el promedio de las relaciones de capacidades específicas a presión y volumen constante de entrada y salida 𝑘, la cabeza politrópica 𝑦𝑝 y la diferencia de entalpías de entrada y salida ℎ1 𝑦 ℎ2 respectivamente.

2.2. Dimensionamiento.

Se inició por determinar la velocidad lineal en la salida 𝑢. del impulsor, haciendo uso de la ecuación 9, la cual fue deducida de la investigación realizada por Lüdtke [4].

Ecuación 10

El radio externo 𝑟. del impulsor se halló usando la ecuación 11 que muestra la relación entre velocidad lineal y velocidad angular 𝑁 puesto que la velocidad angular es uno de los parámetros de entrada.

Ecuación 11

Se calculó entonces el coeficiente de flujo con la ecuación 12, de manera que se encuentre dentro del rango estipulado. Lüdtke [4] Muestra que los bajos coeficientes de flujo llevan a impulsores ineficientes debido a que al reducir la sección transversal por la cual circula el fluido se aumentan las pérdidas de energía.

Ecuación 12

Siendo 𝑚˙ el flujo másico. Este valor adimensional puede tomar cualquier valor que se encuentre dentro del rango descrito, sin embargo, es mostrado por Lüdtke [4] y Aungier [2]; [3] que las mejores condiciones de operación son alcanzadas cuando el valor del coeficiente de flujo se acerca a 0.08.

La metodología sigue al calcular el número de aspas 𝑧, ecuación 13, puesto que esto permitió hallar una de las relaciones trigonométricas del triángulo de velocidades.

Ecuación 13

Donde 𝑃𝑟𝑡 corresponde a la razón de presiones totales. Esta ecuación junto con las ecuaciones 14 a 17, presentadas por Aungier [2], permiten usar relaciones trigonométricas para calcular las dimensiones generales del impulsor, dentro de las cuales se encuentran ángulo de entrada𝛽1 ángulo de salida 𝛽2, radio de entrada 𝑟1.

Ecuación 14

Ecuación 15

Ecuación 16

Ecuación 17

Es importante destacar que la ecuación 14 es un límite teórico al radio de entrada, es decir, esta ecuación se usa con el fin de determinar que el radio de entrada no es igual al radio de salida.

Trazando los respectivos triángulos de velocidad, los elementos anteriormente mencionados fueron derivados, se hace la aclaración de que el radio de entrada que se está hallando corresponde a un radio medio, es decir, el área de entrada alude a la forma de toroide, donde se tiene un radio interior y exterior, el radio de entrada anteriormente calculado se sitúa justo en la mitad del radio interior y exterior. Repetir el cómputo de los triángulos de velocidad usando los radios interior y exterior completa el dimensionamiento del aspa.

Para la altura del impulsor (𝐴𝑙𝑖) se utilizó la ecuación 18 [3], la longitud de arco 𝑎𝑐 requerida para el canal se computó con la ecuación 19, la altura del canal de salida 𝑎𝑙𝑐 con la ecuación 20 y el espesor del aspa 𝑏𝑡 se determina con la ecuación 21 [13].

Ecuación 18

Ecuación1 9

Ecuación 20

Ecuación 21

2.3. Cálculo de pérdidas

Determinar las pérdidas de energía del fluido permite al diseñador modificar la geometría en fases tempranas, ahorrando tiempo en simulación, manufactura y experimentación. Existen diferentes ecuaciones que predicen los mecanismos de pérdida dentro del impulsor, sin embargo Gutiérrez [11], estudió y comparó diferentes ecuaciones con el fin de encontrar el set que más se acerca a la realidad, encontrando que la mejor combinación de ecuaciones da como resultado un error máximo de 0.6% y un error mínimo de 0.1%.

Debido al uso de modelos CAD las pérdidas por tolerancias fueron obviadas, lo cual no genera ningún inconveniente con el cálculo de las otras pérdidas de energía debido a que son ecuaciones independientes entre sí. Con lo cual, los mecanismos de recirculación, incidencia del fluido, carga de las aspas y rugosidad de superficie fueron tenidos en cuenta. Las ecuaciones predicen pérdidas de energía que son descontadas de la cabeza politrópica anteriormente calculada.

2.3.1. Pérdidas por recirculación

La recirculación, denotadas para este trabajo como ℎ𝑝𝑟, esta asociada con un flujo que sale del impulsor y reingresa, se calculó por medio de las ecuaciones 22-24 [11].

Ecuación 22

Ecuación 23

Ecuación 24

Donde 𝛼. representa el ángulo de salida real del fluido en radianes, 𝐷𝑓 es el factor de difusión, 𝑣𝑟𝑡1𝑒 es la componente tangencial de la velocidad relativa en el radio exterior de la entrada, 𝑣𝑟𝑡2 la componente tangencial de la velocidad relativa en la salida; ℎ𝑎𝑒𝑟𝑜 corresponde a la entalpía aerodinámica, 𝑟.𝑒 es el radio externo de la entrada y 𝑐𝑝1, 𝑐𝑝2 son las capacidades de calor específico a la entrada y salida respectivamente [11].

2.3.2. Pérdidas por incidencia del fluido

Se representa en este trabajo como 𝒉𝒑𝒊𝒇 compete al ajuste del flujo al ángulo de entrada del impulsor, además tiene en cuenta el cambio de velocidad abrupto del fluido cuando tiene el primer contacto con el aspa [11]. Se usa la ecuación 25 para calcular este mecanismo de pérdida.

Ecuación 25

Donde 𝑓𝑖𝑛𝑐 es el factor de incidencia, el cual puede tomar valores entre 0.5 y 0.7 [11]. Se asumirá el tope superior, simulando las condiciones de peor operación.

2.3.3. Pérdidas por carga de las aspas

Denotada como ℎ𝑝𝑐𝑎, es la pérdida generada por un gradiente de velocidad negativo que se crea entre el fluido y el impulsor, que tiene como repercusión un incremento de la capa límite de velocidad [11]. Se calculó con la ecuación 26.

Ecuación 26

2.3.4. Pérdidas por rugosidad de superficie

Toma la abreviatura de ℎ𝑝𝑟𝑠, y su influencia se predice con las ecuaciones 27 a 30, y compete a la pérdida de energía derivada de los esfuerzos cortantes en las paredes del impulsor [11].

Ecuación 27

Ecuación 28

Ecuación 29

Ecuación 30

Donde, 𝑙𝑏 es la longitud del canal por el cual se desplazará el fluido, estimada por medio de la medición sobre el modelo CAD generado, y 𝐷ℎ𝑦𝑑 es conocido como el diámetro hidráulico, con 𝐴 siendo el área del canal y 𝑝el perímetro del área que se encuentra en contacto directo con el fluido, cuyo valor se promediará entre la entrada y salida debido al cambio de sección transversal y 𝐶𝑓 es el coeficiente de fricción superficial, el cual debe ser estimado con la ecuación de Colebrook [14], tomando un valor de rugosidad superficial de 0.0000015 [m] debido al proceso de rectificado al que se someten las piezas posterior al maquinado en frio [15].

3. Diseño de la voluta

La metodología expuesta por Pfleiderer [9], para el diseño de voluta, requiere generar secciones radiales parciales de la voluta las cuales se muestran de forma general en la figura 2.

A cada una de estas secciones se le computó un momento angular que debe permanecer constante para todas las secciones y debe tener un valor superior a 1 [9], usando la ecuación 31 donde M denota el momento anular, 𝑅𝑣 el radio desde el centro de rotación hasta el centroide de la voluta y 𝑉𝑦 la velocidad del fluido en dicha sección.

Ecuación 31

Debido a que el fluido de trabajo es compresible al ser aire ideal, es decir, con una composición química de 71% nitrógeno, 29% oxígeno y sin humedad [16], y se esperan cambios drásticos en la velocidad del fluido al incrementar la sección transversal, se hace necesario modelar el fluido como compresible. Las propiedades del flujo se relacionan con el momento angular por medio de un balance de masa sobre cada sección, el cual toma la forma de la ecuación 32.

Ecuación 32

Donde 𝑉, 𝜌, 𝐴 denotan la velocidad, densidad y área respectivamente, siendo x e i son los subíndices para las propiedades de la sección anterior y propiedades del fluido que entra del impulsor, estos elementos se pueden visualizar adecuadamente en la figura 1 donde se hace una vista en detalle de lo que compete a una sección. Así mismo, 𝑃𝑜 𝑦 𝑇𝑜 son la presión y temperatura de estancamiento respectivamente; 𝑘 es la relación de capacidades específicas; 𝑐𝑝 es la capacidad específica a presión constante a 75[°C] debido a que su variación respecto al cambio de las propiedades se considera despreciable. R es la constante del gas, en nuestro caso es aire.

La ecuación 32 debe ser iteró usando el método de la secante presentado en la ecuación 33 con el fin de obtener el valor de la velocidad y con esta determinar las propiedades termodinámicas y geométricas de la sección transversal.

Ecuación 33

Debido a su rápida convergencia y evitar el uso de aproximaciones numéricas para las derivadas en caso de usar el método de Newton-Raphson se consideró que el método de la secante para la solución numérica genera menos errores relativos y es más sencillo de aplicar [17].

4.Modelado y simulación del impulsor

4.1Modelado del impulsor

En la tabla 1 se presentan las propiedades generales con la cuales debe contar el impulsor, solo se permiten modificaciones a la velocidad angular para este caso de estudio, debido a que de esta forma se disminuyen considerablemente las iteraciones que se deben realizar para obtener una solución. Tampoco se incluyen restricciones de tamaño, temperatura o resistencia de materiales.

Usando las ecuaciones 1 a 9 se definen completamente las propiedades termodinámicas, además se determina en una fase temprana si la eficiencia máxima es aceptable. En la tabla 2 se presentan los resultados obtenidos, los cuales pertenecen a los rangos establecidos por los autores e esta investigación después de revisar la información presentada por Aungier [2]; [3] y Lüdtke [4].

El siguiente paso en la metodología de diseño fue determinar los parámetros geométricos del impulsor, es decir, las dimensiones específicas, para lo cual se usaron las ecuaciones (10-21) además de las relaciones trigonométricas necesarias para resolver los triángulos de velocidad, los resultados son mostrados en la tabla 3.

Los elementos de la tabla 3 describen el impulsor en sus rasgos más generales, es importante destacar que 𝑟.𝑖 compete al radio de entrada interno y 𝑟.𝑒 al radio de entrada externo, puesto que la entrada del impulsor está delimitada por 2 círculos concéntricos. Para determinar estos valores fue necesaria la ecuación 34, la cual relaciona el área de entrada, y una distancia equivalente 𝑥 desde el radio de entrada 𝑟. y los radios de entrada interno y externo.

Ecuación 34

Computar los ángulos de entrada para los radios de entrada interno y externo se realiza por medio de los triángulos de velocidad. Algunos investigadores como Mojaddam, Et al. [1], Aungier [2], Lüdtke [4], Meroni Et al. [5], confían en definir 3 ángulos de entrada para obtener por completo el perfil del aspa. Sin embargo, en trabajos más modernos presentados por Casey Et al. [6] y Hildebrand, Et al. [8] muestran que más ángulos de entrada deben ser definidos.

Para generar el modelo con las especificaciones calculadas se usa ANSYS BladeGen, complemento con el cual se puede obtener un modelo CAD de las aspas, a las cuales se les agrega la cubierta superior “shroud” e inferior “hub”. Haciendo uso de SolidWorks 2019. La figura 3 muestra un render del modelo final obtenido.

4.2. Simulación

El proceso de simulación inició con el mallado de la geometría, para obtener mejores resultados disminuyendo el costo computacional se tomó ventaja de la simetría rotacional que presenta el impulsor, con lo cual se puede simular un canal de flujo que contenga un aspa como se muestra en la figura 4. El mallado de esta geometría se da por medio de TurboGrid, complemento de ANSYS, y se impone un límite de 512000 elementos o nodos. Un análisis de sensibilidad para determinar el refinamiento de malla es realizado en conformidad con los lineamientos de la ASME [18], y revela que se requieren alrededor de 506000 elementos para obtener resultados fiables sin superar el límite impuesto. Los puntos críticos donde se debe refinar la malla son, el borde de ataque y las paredes del aspa, debido a que son los puntos que tienen mayor incidencia en el comportamiento del fluido.

Las condiciones de frontera para la simulación fueron una presión total de entrada, con un valor de 100[kPa]; los vectores de velocidad se mantienen perpendiculares a la entrada y una temperatura de entrada constante de 25°C; flujo másico en la salida del impulsor, el cual tiene un valor variable desde 0.2 hasta 0.6 [kg/s] para poder generar las curvas de desempeño. Se seleccionó un modelo de turbulencia tipo SST debido a que combina la precisión del modelo 𝑘 − 𝜔 en cercanías a las paredes del modelo y el desarrollo de la turbulencia en del modelo 𝑘 − 𝜀 [19],[20]; Para simular el comportamiento de la compresibilidad del gas se usa un modelo de energía total, activando los efectos de alta velocidad.

En la figura 5 se muestran las curvas características obtenidas por medio de simulación, en el eje principal se muestran las razones de presión total de la salida y la presión total de la entrada. En el eje secundario están asociadas las eficiencias isentrópicas y la potencia consumida. Es importante destacar de la figura 5 que la relación de presiones tiene el comportamiento esperado al presentar una tendencia al descenso, sin embargo, la eficiencia aumenta después del punto de diseño considerado teóricamente como el punto de mayor eficiencia, el error es causado porque la densidad de malla para flujos mayores al punto de máxima eficiencia no es la adecuada ya que estos son puntos de alta turbulencia.

Para validar los resultados obtenidos por medio de la simulación, se replicaron los resultados obtenidos por Mojaddam Et al.[1], usando la metodología expuesta en esta investigación. Los resultados son presentados en la tabla 4 en los cuales se muestran las comparaciones de eficiencia isentrópica obtenidas en [1] y la respectiva validación realizada por los autores.

5. Modelo y simulación de la voluta.

5.1. Modelo de la voluta

Se usó Matlab para iterar la ecuación 32 con el método de la secante como se ha expuesto anteriormente. Se determinó que el número de secciones necesarias es 100 puesto que es la forma de mantener la continuidad de la curvatura en el modelo CAD realizado en SolidWorks. La ecuación 32 determina la componente axial de la velocidad del flujo, con la cual se definió el estado termodinámico y la geometría de la voluta. En la figura 6 se muestra la variación del radio para la respectiva sección.

5.2. Simulación de la voluta

Con la información presentada en la figura 6 se usó SolidWorks para generar una geometría con la cual se simuló el comportamiento del fluido. Se usó una línea tangente a la circunferencia interna de la voluta, con el fin de tener una guía para la salida del fluido. Se agregó una sección adicional para la simulación, en Spaceclaim, para evitar errores numéricos por la violación de las condiciones de frontera establecidas. Se usó una malla de tipo hexaédrica con capas de inflación sobre las fronteras geométricas para obtener mejores resultados en la capa límite y en el desarrollo de la turbulencia a partir de esta. Un análisis de sensibilidad muestra que 225000 elementos y 85700 nodos son suficientes para realizar simulaciones adecuadas. En la figura 7 se presenta el modelo que se usó para mallado.

Debido a que la unión del impulsor y la voluta supera el límite de elementos impuesto por los autores, se decidió simular ambos elementos por aparte, manteniendo velocidades, temperatura y presión estática de la salida del impulsor como entrada de la voluta, este acercamiento genera errores considerando el alto número de Reynolds [21]. Se usó el modelo de turbulencia tipo SST, debido a que combina los modelos de turbulencia 𝑘 − 𝜔 y 𝑘 − 𝜀 cuyas propiedades en el desarrollo de la simulación han sido mencionadas anteriormente y hacen del modelo SST la mejor opción, además diversos autores lo usan para la simulación de volutas [15]; [20]; [21]. Se usó el flujo másico como condición de salida, para mantener la continuidad. Se emplearon monitores de presión estática y presión total para determinar convergencia de las variables.

Los resultados obtenidos son presentados en la tabla 5 donde se comparan las presiones totales antes y después de pasar por la voluta con el fin de determinar las pérdidas en el flujo.

Se observó que las pérdidas de presión en la voluta son ínfimas, sin embargo, otros análisis muestran que el perfil de velocidades obtenido no es de la mejor calidad, como se observa en la figura 8. Óptimamente se busca un perfil anular sobre la sección mostrada, donde la velocidad máxima este en el centro. Los resultados pueden ser revisados teniendo en consideración las turbulencias generadas en el impulsor.

6. Conclusiones

Se diseñó el impulsor de un compresor centrífugo implementando una combinación de metodologías que permiten limitar el uso de relaciones empíricas, así como eliminar el uso de interpolaciones gráficas que puedan llegar a generar fuentes de error.

Se diseñó la voluta de un compresor centrifugo haciendo uso de la conservación de momento angular, metodología expuesta en [3]; [10]; [15] para lo cual se solucionaron ecuaciones de flujo compresible por medio de métodos numéricos.

Se simuló la geometría diseñada con la metodología presentada en el documento y se calculó el error relativo entre los cálculos realizados para el flujo másico de 0.4[𝑘𝑔/𝑆] y las simulaciones. Se obtuvo un valor de 3,26% de error para las eficiencias; 3.31 % de error para la temperatura; 13% de error para la presión total de salida y 7,14% de error para la presión estática de salida.

Se graficaron las presiones totales, potencias consumidas y eficiencias en función de los flujos másicos, y se encontró el mismo comportamiento reportado en la literatura para las presiones totales y las potencias consumidas. Simulaciones adicionales son requeridas para confirmar la eficiencia.

Se calculó, por medio de la simulación, que la energía perdida al pasar por la voluta es 0,75% de la presión total que tiene el flujo.

Se validó la metodología de simulación por medio de la comparación de los datos obtenidos por Mojaddam Et al [1] con datos generados por el autor, encontrándose que los valores de error para las eficiencias isentrópicas corresponden a 3.89% mientras que los errores promedio para la relación de presiones totales fueron de 11.69%. Se advierte que los errores promediados se obtuvieron excluyendo los resultados obtenidos para un flujo másico de 0.08[𝑘𝑔/𝑠], debido a que este representa un dato estadísticamente no valido.

Se observó que 2 vórtices son los encargados de generar un flujo no uniforme en la salida de la voluta, se advierte que estos resultados se pueden modificar si se llega a obtener el acople entre la voluta y el impulsor en una sola simulación; puesto que las interacciones entre las turbulencias derivadas del proceso de aceleración en el impulsor son obviadas para la presentada en este documento.

Referencias

Información adicional

Cómo citar: J.F. Rincón, G.F. García, “Metodología de diseño y dinámica computacional de fluido del impulsor y voluta de un compresor centrífugo,” Rev. UIS Ing., vol. 19, no. 1, pp. 49-58, 2020. doi: 10.18273/revuin.v19n1-2020004

Revista UIS Ingenierías
ISSN: 1657-4583
Vol. 19
Num. 1
Año. 2020

Metodología de diseño y dinámica computacional de fluido del impulsor y voluta de un compresor centrífugo

Juan FelipeGabriel Fernando Rincón-FrancoGarcía-Sánchez
Universidad Pontificia BolivarianaUniversidad Pontificia Bolivariana,ColombiaColombia