ES2933618A1

ES2933618A1 - Métodos de optimización de diseños asistidos por aprendizaje automático y micromezcladores mecánicos diseñados con dichos métodos

Info

Publication number: ES2933618A1
Application number: ES202130456A
Authority: ES
Inventors: Ortiz Francisco Javier Granados; Casanova Joaquín Ortega
Original assignee: Universidad de Malaga
Current assignee: Universidad de Malaga
Priority date: 2021-05-20
Filing date: 2021-05-20
Publication date: 2023-02-10

Abstract

Métodos de optimización de diseño asistidos por aprendizaje automático y micromezcladores mecánicos diseñados con dichos métodos. La presente invención se refiere a un método de optimización de diseño asistido por aprendizaje automático que comprende generar un conjunto de datos iniciales mediante técnicas computacionales de bajo coste o mediante datos de la literatura, crear un predictor para guiar el proceso de optimización de manera eficiente en la creación de un modelo sustituto, y refinar dicho modelo sustituto apoyado en el predictor. Adicionalmente, la invención refiere micromezcladores mecánicos diseñados usando dicho método, para, mediante desprendimiento de vórtices, conseguir la máxima eficiencia de la mezcla a la par que la mínima caída de presión.

Description

DESCRIPCIÓN

Métodos de optimización de diseño asistidos por aprendizaje automático y micromezcladores mecánicos diseñados con dichos métodos

CAMPO DE LA INVENCIÓN

Esta divulgación hace referencia a métodos de optimización de diseño asistidos por aprendizaje automático de utilidad para procesos de diseño en ingeniería mecánica a cualquier escala. Esta divulgación hace referencia además al diseño de micromezcladores mecánicos utilizando tales métodos.

ANTECEDENTES DE LA INVENCIÓN

En aplicaciones cotidianas de ingeniería mecánica, a menudo es complejo lograr un diseño óptimo para múltiples objetivos, debido a los costes relacionados con la fabricación y prueba de diferentes prototipos. Por esta razón, el uso de herramientas computacionales es una práctica recomendable.

Lograr un mezclado adecuado es crucial en muchas aplicaciones médicas, químicas, biológicas e industriales, como en dispositivos médicos para prevenir la formación de trombosis, la mejora de reacciones en aplicaciones químicas y farmacéuticas, cámaras de combustión de turbinas de gas o el aumento de transferencia de calor usando generadores de vórtice y nanofluidos. Además, la micromecánica y microfluídica es un campo de investigación que ha ganado mucha atención en los últimos años, debido a sus múltiples aplicaciones en las áreas mencionadas. Entre algunas de estas aplicaciones se pueden destacar el uso de dispositivos LOC (Lab On a Chip), que son pequeños dispositivos que integran una o más funciones de laboratorio en un solo chip. Dentro del campo del mezclado de fluidos, este tipo de dispositivo tiene ciertas ventajas, como el bajo consumo de fluidos, tiempos de respuesta rápidos debido a sus pequeñas dimensiones y la posibilidad de producción en masa. Como ejemplos de dispositivos LOC se pueden mencionar aquellos utilizados en análisis clínicos para la detección temprana del cáncer o para la detección de enfermedades como el VIH, la malaria o muchas otras patologías.

En el frente opuesto a los microdispositivos pasivos, donde no se requiere suministrar energía adicional (excepto aquella para lograr superar la diferencia de presión para el movimiento del fluido), se encuentran los microdispositivos activos. En los microdispositivos activos es necesario un consumo de energía adicional para mejorar la eficiencia de dichos dispositivos. Los parámetros de interés para el diseño óptimo del mezclador pasivo son las dimensiones del pilar, en términos de la relación de aspecto AR y la relación de bloqueo BR, y el número de Reynolds Re. La posición axial Lu del pilar con respecto a la entrada no es relevante, ya que se asume que el flujo está completamente desarrollado a lo largo del canal 2D de ancho H. También se asume que la longitud aguas abajo es lo suficientemente larga para que no se note ningún impacto en la simulación debido a la condición de contorno a la salida.

El mecanismo fluidomecánico responsable de una mezcla adecuada en este dispositivo es el desprendimiento de vórtices generado aguas abajo del pilar. Este conocido comportamiento del flujo se ha estudiado durante décadas. Algunos autores se han centrado, por ejemplo, en el uso de diferentes formas de pilares, en el análisis de la frecuencia de desprendimiento, en la emisión de ruido, en el efecto de mezcla de especies/calor, etc. A pesar de que el desprendimiento de vórtices se ha estudiado ampliamente en la literatura, todavía hay muchas investigaciones recientes y en curso sobre la física de flujos en relación a este fenómeno. Por ejemplo, en algunos trabajos se encuentran relaciones interesantes para demostrar cómo el número de Reynolds crítico, el número de Reynolds crítico individual y el número de Roshko están influenciados por la fracción sólida. Otros investigadores han orientado recientemente sus esfuerzos a estudiar la supresión del desprendimiento de vórtices cuando este es problemático. Sin embargo, a pesar de la gran evidencia de investigación de este fenómeno en la literatura, todavía hay aplicaciones de este problema de flujo bajo estudio para su mejora.

Un punto de partida para la presente invención es poder predecir si el desprendimiento de vórtices puede tener lugar o no. Es importante saber de antemano si un diseño conducirá a este fenómeno, por la importancia en el proceso de mezcla y las limitaciones en el AR, BR y Re. Dado que la potencial aplicabilidad del dispositivo es para micromecánica y microfluídica en aplicaciones médicas, farmacéuticas o biológicas, el Re máximo, debido a las dimensiones muy reducidas del dispositivo, estará restringido a Re < 200. Por otro lado, AR y BR están restringidos por la geometría del pilar y el propio canal. Por tanto, bajo estas limitaciones, puede resultar complejo lograr un buen diseño. Para predecir la aparición de desprendimiento de vórtices mediante esta invención, se utilizan algoritmos de clasificación de aprendizaje automático. Debe destacarse que la presente investigación tiene como objetivo lograr el desprendimiento de vórtices para lograr la mejora de la mezcla. Sin embargo, los modelos predictivos de aprendizaje automático proporcionados en este documento se pueden utilizar para evitar diseños que conduzcan a la aparición de vórtices debido a posibles efectos no deseados, por ejemplo en la eficiencia de los elementos aguas abajo, generar resonancia o emisión de ruido. Dentro de este contexto de diseño, el uso de este tipo de modelos predictivos puede ser un tema de interés en el modelado de Gemelos Digitales, que consiste en la integración de diseños virtuales y físicos para compartir datos entre ellos de manera continua. Esto convierte a los Gemelos Digitales en una de las tecnologías más prometedoras para revolucionar la fabricación inteligente en la industria 4.0.

Concretamente, Regresión Logística y Bosque Aleatorio son las técnicas de clasificación utilizadas en la presente invención. La Regresión Logística (LR, del inglés Logistic Regression) se utiliza como modelo para asignar una probabilidad, que es un número que varía entre 0 y 1. La LR tiene la siguiente forma:

donde f es el vector del espacio de variables de diseño y la función de regresión z ( f) se modela como:

donde at son los coeficientes a estimar y z¿(f) son los términos de la parte literal de los monomios. Estos pueden diseñarse para una regresión lineal o no lineal. z¿(f) corresponde al espacio de las variables de diseño & en el caso de modelos de regresión lineal de primer orden. Los polinomios de orden superior o las transformaciones de variables también se pueden probar en la regresión logística lineal. Una vez calculada la probabilidad, se puede definir un umbral por encima del cual los casos se pueden clasificar como positivos (1) y viceversa para definir los negativos (0). Este umbral suele ser el valor intermedio 0.5, pero puede cambiar según la aplicación. La LR es un modelo predictivo muy popular en las ciencias médicas y biológicas, debido a la facilidad para interpretar los resultados, al contrario de lo que sucede con, por ejemplo, los métodos de aprendizaje profundo. Sin embargo, no se han encontrado aplicaciones en la literatura de la LR en aerodinámica industrial similares a nuestro estudio.

Otro modelo de predicción interesante, pero más complejo, es el Bosque Aleatorio (RF, del inglés Random Forest). RF es un método de árbol de decisión avanzado que consiste en una combinación de árboles en una forma de aprendizaje en conjunto, y mediante la creación de estos árboles predictores basados en un muestreo aleatorio (es decir, un bootstrap con o sin reemplazo, lo que significa volver a seleccionar o no las muestras ya seleccionadas) sobre el set de datos original. El aprendizaje en conjunto significa utilizar diferentes predictores entrenados para la misma tarea para combinarlos. El método del RF se ha establecido en la literatura del aprendizaje automático como un método muy robusto y preciso, convirtiéndose en uno de los algoritmos de predicción más populares en la era actual de Big Data, como se puede ver en la gran cantidad de trabajos en dicha literatura.

Hasta donde conocen los inventores, no existen trabajos en la literatura sobre el uso de estas técnicas de clasificación para predecir la generación de desprendimiento de vórtices, lo cual es de vital importancia en esta invención para lograr la mezcla deseada. Los predictores también permiten decidir si simular o no un nuevo caso (punto de muestra) para las funciones a evaluar en el proceso de optimización.

El aprendizaje automático es un campo amplio, que incluye más aplicaciones que la clasificación. Muchos autores han aplicado el aprendizaje automático en la optimización basada en modelos sustitutos o en la cuantificación de la incertidumbre (UQ, del inglés Uncertainty Quantification). Entre estas técnicas utilizadas en optimización o UQ se pueden encontrar los Procesos Gaussianos o Interpolación de Kriging, que se utilizan en esta invención para generar modelos sustitutos, también conocidos como superficies de respuesta o metamodelos. Estos sustitutos son el punto de partida para el proceso de optimización. Kriging tiene como objetivo obtener una aproximación de un modelo exacto de la forma:

donde y es el modelo exacto, y es el modelo sustituto y e el error entre el modelo sustituto y el modelo exacto, todos definidos en el espacio de parámetros de diseño f . Kriging se implementa en dos pasos: primero se genera una función de regresión / ( f ) basada en el conjunto de datos, y que está destinada a capturar la mayor variación en los datos (es decir, la “tendencia” general). En segundo lugar, se construye un proceso gaussiano Z (f) , obteniendo la expresión:

donde / ( f ) es el vector de dimensiones k x 1 que contiene las funciones de la base de la regresión [/" (O f 2(O ... f$ ( f ) ] y ^y¿ se refiere a los coeficientes de la regresión. Dependiendo del orden de la regresión, el método se llama Kriging Universal (f(^) es un polinomio multivariante de orden n > 1), Kriging Ordinario (f(£) = Yo, con y o como coeficiente desconocido a estimar) o Kriging Simple ((£) = constante).

Una vez que los modelos sustitutos están disponibles como funciones sin apenas coste computacional en su evaluación, deben explorarse para lograr el diseño óptimo. Para este propósito, se utiliza el popular algoritmo NSGA-II, denominado algoritmo genético de ordenación no dominada II, que es una versión mejorada del algoritmo NSGA original. NSGA-II es un poderoso algoritmo genético (GA, del inglés Genetic Algoríthm) con una formulación elitista. El elitismo significa que cada solución se "clasifica" de acuerdo con el nivel de no dominación, en aras de apuntar a soluciones no dominadas de manera más directa. Los GA son populares debido a su capacidad para imitar las ideas evolutivas de la selección natural y la genética de acuerdo con la teoría de la "supervivencia del más apto". Esto hace que el método sea muy útil como algoritmo de búsqueda adaptativo para problemas de optimización, incluso para optimización no lineal en problemas de ingeniería complejos. Estos algoritmos son robustos y tienen un buen rendimiento en escenarios de optimización complejos. La interpolación de funciones de base radial (RBF, del inglés Radial Basis Functions) también es un método popular en la literatura para definir metamodelos adecuados para la optimización.

DESCRIPCIÓN DE LA INVENCIÓN

El proceso de optimización descrito aquí, denominado optimización de diseño asistido por aprendizaje automático (MLADO, del inglés Machine Learning-Aided Design Optimisation), se desarrolla de acuerdo con el diagrama presentado en la Fig. 4. El marco MLADO consiste en usar simulaciones CFD, un modelo predictivo, un modelo sustituto y un algoritmo de optimización. Sobre los datos iniciales considerados para construir los modelos sustitutos, se tiene que considerar si estos son suficientes datos y si los modelos sustitutos tienen suficiente calidad. En principio, se debe decidir dónde colocar los nuevos puntos de datos (es decir, nuevas simulaciones CFD) en función del error cuadrático medio (MSE), R2 o cualquier otra medida estadística proporcionada por el método sustituto considerado. Al incluir un predictor para clasificación, el proceso para lograr modelos sustitutos útiles puede ser guiado por los datos (data-dríven) y ahorrar importantes recursos computacionales. Si el resultado del predictor es que la nueva configuración potencial es válida, entonces se considera con solidez como un nuevo punto para la mejora del modelo sustituto y por tanto se simula.

La primera etapa del proceso es adquirir datos suficientes para construir un sustituto inicial. Aquí se puede utilizar cualquier Diseño de Experimento (DoE, del inglés Design o f Experíment), dependiendo de la naturaleza del problema, el presupuesto computacional, etc. Sin embargo, se recomienda explorar inicialmente al menos tres puntos por dirección (i.e. variable en el espacio de diseño), ya que la respuesta puede ser no lineal. Tras esto, se construyen los modelos sustitutos. Como es probable que los modelos sustitutos iniciales sean aproximaciones un poco groseras, algunas regiones deben refinarse. Para tal tarea, el error cuadrático medio (MSE, del inglés Mean Square Error) indica si es necesario refinar o no. Si es necesario un refinamiento, el predictor de desprendimiento de vórtice (modelo de clasificación) filtra si la nueva configuración es útil o no. Si es así, el punto se puede simular y agregar al conjunto de datos sustitutos compuesto por Ns muestras.

Debe destacarse que Np (número de muestras del predictor de desprendimiento de vórtice o predictor de configuración) y Ns (número de muestras para el modelo sustituto) no tienen por qué coincidir. Esto es así porque la suavidad (smoothness) de la respuesta de la predicción de la configuración deseada puede ser mayor que la suavidad en e.g. los sustitutos de Kriging de la eficiencia del proceso. Esto dependerá de las variables de interés y del problema que se esté considerando. Además, el predictor de configuración deseado no requiere del entrenamiento con simulaciones completas. Tan solo es necesario diferenciar entre casos válidos y no válidos, y esto se puede hacer, por ejemplo, a partir de simulaciones de baja fidelidad o datos externos.

El predictor de configuración deseado es un algoritmo previamente entrenado. Este predictor se puede entrenar con datos de cualquier naturaleza (simulaciones de problemas no estacionarios atendiendo a su comportamiento como simulación estacionaria, mallas más gruesas, simulaciones de baja fidelidad, soluciones analíticas, bases de datos, otra literatura, etc.). Una vez en el bucle del MLADO, actualizar el modelo predictivo para clasificación es una opción recomendable, pero hay que tener precaución. El aprendizaje del algoritmo de clasificación con datos desequilibrados es un problema a destacar. Si el predictor de configuración deseado se actualiza automáticamente cuando se simulan completamente nuevas muestras no consideradas en el entrenamiento del clasificador, entonces puede aparecer un desequilibrio en el entrenamiento. Dado que los nuevos datos de entrenamiento se simularían en función de la predicción de configuraciones válidas, se aumentaría solo el número de casos de estas configuraciones válidas, aumentando así el desequilibrio en el entrenamiento. Solo los falsos positivos del predictor mejorarían el equilibrio en el conjunto de datos, ya que la simulación incluiría así datos con configuraciones no válidas. Por este motivo, no se recomienda actualizar automáticamente el clasificador y es más recomendable utilizar un predictor bien entrenado desde el principio si es posible, el cual debería actualizarse solo según criterio experto.

Después de la ejecución del bucle, se puede dejar que el marco de optimización funcione automáticamente hasta que se logre cierto nivel de refinamiento. El criterio de convergencia para esto está abierto al investigador y se puede controlar manualmente. Cuando se logran unos modelos sustitutos suficientemente buenos, se usa el algoritmo de optimización para evaluarlos y encontrar un frente de Pareto para el problema multiobjetivo. Este marco MLADO podría extenderse y automatizarse a cualquier proceso de diseño de ingeniería sujeto a clasificación, así como a otros procesos de diseño similares en ingeniería mecánica a cualquier escala, si se incluyen estrategias computacionales de parametrización de formas. Una característica interesante para agregar al marco MLADO puede ser incluir una opción de transformación de malla (mesh morphing). De hecho, con esta herramienta se podría definir un modelo de orden reducido (reduced order modelling) del cálculo e inspeccionar los resultados rápidamente y en tiempo real.

El marco ha demostrado ser una forma eficiente para el diseño óptimo de dispositivos micromecánicos de mezcla que consisten en un pilar rectangular confinado en un microcanal diseñado para generar desprendimiento de vórtices. Un clasificador de Bosque Aleatorio (RF) está entrenado para predecir qué configuración geométrica puede provocar desprendimiento de vórtices. Posteriormente, se investiga el problema de optimización multiobjetivo, que consiste en minimizar la potencia de bombeo requerida y maximizar la eficiencia de mezcla bajo algunas restricciones de diseño. Si se desean datos de entrenamiento adicionales para los sustitutos, el clasificador RF se puede usar para predecir si vale la pena o no simular la nueva configuración, evitando ejecutar casos computacionales intensivos irrelevantes y acelerando el proceso basado en datos. Se simulan los diseños óptimos resultantes del uso del algoritmo genético NSGA-II en los sustitutos y se muestra su rendimiento. Las configuraciones geométricas óptimas, incluso para condiciones de mezcla muy desfavorables y un número de Reynolds medio-bajo de 200, dan una eficiencia de mezcla máxima de alrededor del 50% a un coste de potencia de bombeo muy bajo en un canal corto, mejorando a los dispositivos existentes en la literatura, especialmente en mezcla por unidad de longitud, ya que se logra una buena mezcla incluso para un microcanal corto de longitud L = 5 y la interacción con un único objeto.

BREVE DESCRIPCIÓN DE LAS FIGURAS

Figura 1: Croquis de la geometría.

Figura 2: Validación y comparación con otros autores: Turki et al. (2003) [95], Patil y Tiwari (2008) [96] y Sharma y Eswaram (2004) [97].

Figura 3: Detalle de la malla óptima alrededor del pilar con AR = 2.

Figura 4: Marco de optimización MLADO.

Figura 5: Curva ROC del poder de clasificación de la Regresión Logística en (a) datos de test y (b) datos de entrenamiento.

Figura 6: Error OOB del algoritmo Bosque Aleatorio.

Figura 7: Representación de casos simulados (valores reales) con y sin desprendimiento de vórtices y el límite de separación para a) AR = 1, b) AR = 0.5, c) AR = 0.25 y d) AR = 0.125.

Figura 8: Modelos sustitutos para r¡ mostrados para diferentes regímenes de Re. Los cuadrados rojos son configuraciones sin desprendimiento de vórtice y los puntos azules son configuraciones con desprendimiento de vórtice.

Figura 9: Modelos sustitutos para r¡ mostrados para diferentes valores de AR. Los cuadrados rojos son configuraciones sin desprendimiento de vórtice y los puntos azules son configuraciones con desprendimiento de vórtice.

Figura 10: Error cuadrático medio (MSE) de los modelos sustitutos r¡ mostrados para diferentes regímenes de Re. Los cuadrados rojos son configuraciones sin desprendimiento de vórtice y los puntos azules son configuraciones con desprendimiento de vórtice.

Figura 11: Error cuadrático medio (MSE) de los modelos sustitutos n mostrados para diferentes valores de AR. Los cuadrados rojos son configuraciones sin desprendimiento de vórtice y los puntos azules son configuraciones con desprendimiento de vórtice.

Figura 12: Modelos sustitutos de r¡ mediante el uso de diferentes correlaciones para AR = 1.

Figura 13: Modelos sustitutos de r¡ mediante el uso de diferentes correlaciones para Re = 120.

Figura 14: Modelos sustitutos para (n ) mostrados para diferentes valores de AR. Los cuadrados rojos son configuraciones sin desprendimiento de vórtice y los puntos azules son configuraciones con desprendimiento de vórtice.

Figura 15: Modelos sustitutos para (n ) mostrados para diferentes regímenes de Re. Los cuadrados rojos son configuraciones sin desprendimiento de vórtice y los puntos azules son configuraciones con desprendimiento de vórtice.

Figura 16: Modelos sustitutos para (n ) mostrados para diferentes valores de AR. Los cuadrados rojos son configuraciones sin desprendimiento de vórtice y los puntos azules son configuraciones con desprendimiento de vórtice.

Figura 17: Modelos sustitutos para (n ) mostrados para diferentes regímenes de Re.

Figura 18: (a) Frente de Pareto obtenido por el algoritmo NSGA-II para las funciones objetivo y (n ) con un tamaño de población de 500 muestras. (b) Valores de los parámetros óptimos. A, B y C son puntos óptimos seleccionados para estudio como referencia.

Figura 19: Frente de Pareto obtenido por el algoritmo NSGA-II para las funciones objetivo r¡ y (n ) construyendo los sustitutos solo con datos VS = 1. Figura 20: Frentes de Pareto de optimización multiobjetivo con diferente número de muestras Ns para el modelo sustituto.

Figura 21: Representación de los puntos óptimos seleccionados en los sustitutos r¡ y (n ) para un valor constante de Re = 200.

Figura 22: Config. A: a) Fracción másica. b) Líneas de corriente. Config. B: c) c) Fracción másica. d) Líneas de corriente. Config. C: e) Fracción másica. f) Líneas de corriente. En todas las figuras, las líneas de corriente están coloreadas de acuerdo con la velocidad total máxima (rojo) y la velocidad total mínima (azul).

DESCRIPCIÓN DETALLADA DE LA INVENCIÓN Y DE UNA DE SUS REALIZACIONES

Aproximación Numérica y Geometría Computacional

La geometría del microdispositivo seleccionado para la implementación de la invención se muestra en la Fig.1: un pilar rectangular bidimensional (2D) (que mide h y l metros de largo en la dirección x e y, respectivamente) se encuentra en la línea central de un microcanal recto de H metros de ancho. Esto significa una relación de anchura de pilar a canal, es decir, una relación de bloqueo BR para h , y una relación de longitud a ancho de pilar, es decir, la relación de aspecto AR para l; el pilar está centrado en el origen de las coordenadas (x, y), estando la entrada del canal a una distancia Lu metros aguas arriba de la cara frontal del pilar, mientras que el canal tiene L metros de longitud.

La entrada del canal consiste en un flujo laminar completamente desarrollado formado por dos fluidos idénticos con diferentes valores de magnitud escalar Y a mezclar entre sí (Y puede ser, por ejemplo, temperatura o fracción másica), coeficiente de difusividad escalar Dc (para la temperatura Dc sería difusividad térmica a, mientras que para la fracción de masa este parámetro sería la difusividad molecular D), y con densidad y viscosidad conocidas, p y ^, respectivamente (v = ^ /p es la viscosidad cinemática). El flujo ingresa a la entrada del canal con un perfil de velocidad parabólico con una velocidad media U. Por la mitad superior de la sección de entrada entra el fluido con valor nulo de la magnitud escalar Y, es decir, Y = 0 en x = - (Lu l / 2), 0 < y < H / 2, mientras que el otro fluido entra por la mitad inferior con un valor unidad de la magnitud escalar, es decir, Y = 1 en x = - (Lu l / 2), - H / 2 < y < 0. Con este tipo de condiciones de contorno para Y se asume que la variable de interés se reescalaría siempre para que sus valores se muevan entre 0 y 1. Aguas abajo, la calidad de la mezcla de los fluidos dependerá del número de Reynolds, el coeficiente de difusividad escalar y la relación de aspecto y bloqueo de la geometría. Por lo tanto, bajo los supuestos de movimiento transitorio y fluido incompresible, los campos de velocidad v = (u ,v) y presión p del flujo 2D en la geometría descrita anteriormente están gobernados por las ecuaciones de continuidad y momento que, en una notación adimensional, pueden ser escritas, respectivamente, como

V • v = 0,

(5)

mientras que la mezcla se rige por la ecuación escalar de convección-difusión, escrita como

La longitud, velocidad, presión y tiempo característicos utilizados en el estudio son H,U,pU2 y H/U, respectivamente. Re en la ecuación (6) es el número de Reynolds definido, como es habitual en los problemas de flujo del canal, como Re = pUH / ^, mientras que Pe en la ecuación (7) es un número similar al Peclet definido como Pe = UH/DC , que es una relación entre dos tiempos característicos: el difusivo, H2/D C; y el convectivo, H/U. Los valores altos de Pe indican que el tiempo necesario para la mezcla por difusión es muy largo y las condiciones de mezcla muy desfavorables. El número de Peclet también se puede escribir como Pe = Re rC, donde rC, definido como rC = v /D C , es una relación entre dos mecanismos de transporte molecular (rC sería el número de Prandtl Pr o el número de Schmidt Se en caso de que el escalar Y fuera la temperatura o la fracción de masa, respectivamente). Se puede ver que, por un lado, Pe depende tanto del fluido como de la geometría, mientras que, por otro lado, rc solo depende del tipo de fluido, siendo este una relación entre dos mecanismos de difusión a nivel térmico/molecular: los de difusión de momento y de temperatura/masa. Para un fluido dado y conocido con propiedades físicas constantes, si las dimensiones de la geometría cambian, Pe y Re pueden cambiar, a pesar del uso del mismo fluido, mientras que rc permanecerá constante. De hecho, un valor alto y muy desfavorable de rc = 104 se ha mantenido constante durante todo el estudio para que haya mayor dependencia de la geometría y menos de la difusión molecular/térmica entre fluidos.

En cuanto a la geometría, las medidas adimensionales son: 1 unidad de ancho (H = 1) transversal, 5 unidades longitudinal (L = 5) y el pilar se ubica a 1 unidad de la entrada del canal (Lu = 1). Además, la relación de bloqueo del canal BR y la relación de aspecto del cilindro rectangular AR son, respectivamente, BR = h / H y AR = l / h. Por lo tanto, para valores dados de Re,AR y BR, el desprendimiento de vórtices del pilar puede aparecer, o no, aguas abajo. En consecuencia, el flujo sería inestable o constante. Por consiguiente, para cada par (AR,BR), hay un número de Reynolds crítico Recr por encima del cual el flujo es inestable y oscilatorio: Recr = Recr (AR, BR). Para un Re < Recr dado, el flujo sería constante, y también las fuerzas de sustentación y arrastre del fluido sobre el pilar. Estas fuerzas, en notación adimensional, se pueden escribir como

donde Cl y Cd son los coeficientes de sustentación y arrastre, respectivamente, y donde h , como es habitual en los estudios de caudales alrededor de pilares, se ha utilizado como longitud característica para tener coeficientes comparables. El comportamiento oscilatorio del Cl debido al desprendimiento del vórtice del pilar cuando Re > Recr, nos permite identificar tanto el valor crítico Recr como la frecuencia de oscilación / del Cl. Esta frecuencia se representa en notación adimensional mediante el número de Strouhal: St = fh /U . Además de esto, si una magnitud adimensional arbitraria g oscila con el tiempo, su valor promediado en el tiempo para un período de oscilación St se puede calcular como

1 i-ío S t"1

<9> = s h { 9{t' )d t '-

(9)

siendo t0 cualquier referencia temporal para el análisis del flujo periódico.

Las pérdidas de carga son un aspecto importante en los flujos circulando a través de un canal. Estas pérdidas definen la potencia de bombeo necesaria para que el fluido sea capaz de fluir a lo largo de dicho canal. Concretamente, esta es la potencia necesaria para mantener un determinado caudal q en el conducto, para el que deben superarse las pérdidas de presión y viscosidad en todo el canal (lo cual se puede evaluar como la diferencia de presión de entrada-salida multiplicada por el caudal másico). De ahí, la potencia de bombeo adimensional, denotada por n, se puede evaluar como

n = A pq,

( 10)

donde Ap es la caída de presión entre la entrada y la salida del canal. Debido a las magnitudes características elegidas, q = 1, y n se ha hecho adimensional usando pU3h. Teniendo en cuenta todas las especificaciones anteriores, finalmente la potencia de bombeo se define por n = Ap.

El parámetro más importante en el diseño del dispositivo es la evaluación de la calidad de la mezcla. Esta característica debe evaluarse en la sección de salida, definiéndose como la eficiencia de mezcla p, en%:

( 11)

donde o es la desviación estándar de la magnitud escalar correspondiente (temperatura o fracción másica) en la salida del canal y omax es la desviación estándar máxima en la entrada: 0.5 en el presente caso. Dado que 0 < Y < 1 y 0 < o < 0.5; o = 0 significa distribución uniforme de Y en la salida, logrando una mezcla completa (^ = 100%), mientras que o = 0.5 significaría que no hay mezcla alguna (^ = 0%).

Las simulaciones numéricas realizadas para simular el funcionamiento del dispositivo como micromezclador se han llevado a cabo con el software comercial ANSYS-Fluent. Con ese fin, el conjunto de Ecs. (5) - (7) se han resuelto numéricamente mediante la formulación basada en presión, con métodos de discretización espacial de segundo orden. Por otro lado, se utilizó el algoritmo SIMPLE (Método semi-implícito para ecuaciones ligadas a presión) para abordar el acoplamiento velocidad-presión. Además, las simulaciones numéricas se ejecutaron con las siguientes condiciones de contorno. La entrada era una condición de contorno velocity-inlet con un perfil de velocidad parabólico impuesto y una fracción másica/temperatura conocida. Para todas las paredes sólidas (pilar y canal), se impuso una condición de contorno de tipo w allcon condición antideslizante (non-slip) con flujo nulo de calor/masa a través de las paredes. A la salida se impuso una presión atmosférica constante mediante la condición de contorno pressure outlet. Esto permite un flujo inverso y no desarrollado, por lo que se considera más adecuado que la condición de contorno de tipo outflow disponible en el software para el problema en estudio. Además, con respecto a la integración temporal de las Ecs. (5) - (7), se ha utilizado un método implícito de segundo orden. El intervalo de tiempo At se eligió para tener el número máximo de Courant por debajo de 1.

Para mayor fiabilidad computacional, las investigaciones numéricas deben realizar tanto un estudio de validación como de convergencia de malla. Estos proporcionan una descripción general de la malla más óptima y permiten cuantificar la incertidumbre de discretización (Grid Convergence Index, GCI). Ambos estudios se resumen en la tabla 1 y la figura 2, respectivamente. Se han estudiado diferentes magnitudes características como St, (Cd), valor pico a pico del coeficiente de elevación (Clpp) y su raíz cuadrada media (Clrms), para tres mallas diferentes con tamaño de celda uniforme (dsi = {0.0125,0.025,0.05} h, i = 1,2,3). Esto se ha probado para cuatro valores de la relación de bloqueo (BR = {1 /8 ,1 /5 ,1 /4 ,1 /3 }), un número de Reynolds fijo (Re = 100) y una relación de aspecto fija

= 1 (pilar cuadrado). Mediante el cálculo de las mencionadas magnitudes de salida se puede estimar su GCI y su Extrapolación de Richardson, cuyo valor corresponde a ds ^ 0 (ver Fig. 2). El GCIi i, i nos permite cuantificar la incertidumbre de discretización de cada magnitud para las dos mallas más finas, es decir, para i = 1,2. Teniendo en cuenta que el índice de refinamiento de la malla es 2, y utilizando un factor de seguridad Fs de 1.25, la incertidumbre de la malla más fina osciló entre el 0.5% y el 2.3%, mientras que para la red media osciló entre el 1.8% y el 6.5%, como se puede ver en la Tabla 1, donde se resumen todos los valores de GCI. Como se ilustra en la tabla, la malla media con ds2 da valores de discretización de incertidumbre razonables y podría elegirse como la óptima. Sin embargo, para reducir aún más la incertidumbre de la discretización, finalmente se ha utilizado un ds menor que ds2: ds = 0.02h. Esto significa un aumento del 20% en el número de nodos de la cuadrícula a lo largo del pilar de ancho h. En la Fig. 3 se muestra una vista detallada de la malla uniforme óptima alrededor del pilar. En cuanto a la bondad de la malla óptima desde el punto de vista de la mezcla convectivodifusivo, cabe mencionar que se debe prestar especial atención a resolver adecuadamente las escalas de longitud más pequeñas en el campo de magnitud escalar F, especialmente en líquidos con Fc normalmente muy grande. En este caso, el tamaño más pequeño de la escala de longitud de campo de fracción másica se puede estimar utilizando la escala de longitud de Batchelor, que en forma dimensional se denota como 1B y se puede escribir como

% Avc/

Ae = 7 f r

( 12)

con Avel, la escala de longitud de las estructuras de velocidad más pequeñas. Suponiendo que para el flujo viscoso de este estudio 2vel ~ H, con rc = 104 ,ÁB ~ 0.01 H = 0.01 h BR~ ", lo que significa que la situación más desfavorable sería la de mayor BR. Dado que en nuestra investigación el valor máximo de BR sería 0.5, esto significa que 1B ~ 0.02h, que es el orden de magnitud de la malla final y óptima elegida.

Los cálculos en el presente trabajo también se han comparado con otros trabajos en la literatura para su validación. En la Fig. 2 se muestra una comparación entre los resultados obtenidos aquí por el modelo numérico y los dados por diferentes autores de la literatura. La buena concordancia entre la presente simulación y los resultados anteriores, además del reducido GCI previamente analizado, evidencia que las consideraciones numéricas en la metodología seguida son las adecuadas para el problema bajo estudio.

i j

Tabla 1: Resultados del estudio de convergencia de malla con los valores de GCI para las mallas indicadas.

Marco de optimización de diseño asistido por Machine Learning (MLADO)

En aras de lograr un dispositivo de mezcla eficiente y óptimo, en esta invención se realiza un estudio de optimización. El proceso de optimización, denominado Optimización del Diseño Asistido por Aprendizaje Automático (MLADO) en este documento, se desarrolla de acuerdo con el diagrama presentado en la Fig. 4. El marco MLADO consiste en utilizar simulaciones CFD, un modelo predictivo, un modelo sustituto y un algoritmo de optimización. Sobre la base de los datos iniciales considerados para construir los sustitutos, se tiene que considerar si estos son suficientes puntos de datos y si los modelos sustitutos tienen suficiente precisión. En principio, se debe decidir dónde colocar los nuevos puntos de datos (i.e. nuevas simulaciones CFD) en función del error cuadrático medio (MSE), R2 o cualquier otra medida estadística proporcionada por el método del modelo sustituto considerado. Sin embargo, dado que la configuración óptima estará orientada al desprendimiento de vórtices para lograr una mezcla adecuada, los diseños sin desprendimiento de vórtices son irrelevantes. Por lo tanto, al incluir un predictor de desprendimiento de vórtice, el proceso para lograr modelos sustitutos útiles puede basarse en datos (datadríven) y ahorrar importantes recursos computacionales. Si el resultado del predictor es que la nueva configuración potencial produce desprendimiento de vórtices, entonces se considera como un nuevo punto para la mejora de los modelos sustitutos y se procede a su simulación.

La primera etapa del proceso consiste en adquirir datos suficientes para construir un sustituto inicial. Aquí se puede utilizar cualquier Diseño de Experimento (DoE), dependiendo de la naturaleza del problema, el presupuesto computacional, etc. Sin embargo, se recomienda explorar inicialmente al menos tres puntos por dirección (variable del espacio de diseño), ya que la respuesta puede ser no lineal. Tras esto, se construyen los modelos sustitutos. Como es probable que los sustitutos iniciales sean poco precisos, algunas regiones pueden tener que refinarse. Para tal tarea, el error cuadrático medio (MSE) indica si es necesario refinar o no. Si el refinamiento es necesario, entonces el predictor de desprendimiento de vórtices (modelo de clasificación) filtra si la configuración es un caso de desprendimiento de vórtices (ES = 1). De ser así, el punto se simula y agrega al conjunto de datos para crear los modelos sustitutos, compuesto por Ns muestras.

Debe destacarse que Np (número de muestras del predictor de desprendimiento de vórtices) y Ns (número de muestras para el modelo sustituto) no es necesario que coincidan. Esto es así porque la suavidad (smoothness) de la predicción del desprendimiento de vórtices puede ser mayor que la suavidad en e.g. los modelos sustitutos Kriging de la eficiencia de mezcla. Esto dependerá de las variables de interés y del problema que se esté considerando, ya que, por ejemplo, la respuesta de la caída de presión no tiene variaciones importantes y es bastante suave. Además, el predictor de desprendimiento de vórtices no requiere el entrenamiento de simulaciones completas. Tan solo es necesario diferenciar entre los casos VS = 0 y VS = 1, y esto se puede hacer, por ejemplo, a partir de simulaciones de baja fidelidad o datos externos.

El predictor de desprendimiento de vórtices es un algoritmo previamente entrenado. Este puede ser entrenado con datos de cualquier naturaleza (simulaciones de problemas no estacionarios atendiendo a su comportamiento como simulación estacionaria, mallas más gruesas, simulaciones de baja fidelidad, soluciones analíticas, bases de datos, otra literatura, etc.). Dado que el desprendimiento de vórtices es un problema en el que las configuraciones con y sin comportamiento oscilatorio pueden agruparse fácilmente, se ha observado que pocas simulaciones son suficientes para entrenar bien al predictor. Una vez en el bucle de trabajo del MLADO, actualizar el modelo predictivo para clasificación es una opción recomendable, pero hay que tener cuidado con esto. El aprendizaje potencial sobre datos desequilibrados es un problema a señalar. Si el predictor de desprendimiento de vórtices se actualiza automáticamente cuando se simulan completamente nuevas muestras no consideradas en el entrenamiento del clasificador, entonces puede aparecer un desequilibrio en el entrenamiento. Dado que los nuevos datos de entrenamiento se simularían en base a la predicción de VS = 1, se incrementaría solo el número de casos de VS = 1, aumentando así el desequilibrio en el entrenamiento. Solo los falsos positivos del predictor mejorarían el equilibrio en el conjunto de datos, ya que la simulación incluiría datos con VS = 0. Por esta razón, no se recomienda actualizar automáticamente el clasificador y es preferible usar un predictor bien entrenado desde el principio, el cual debe actualizarse solo de acuerdo bajo criterio experto.

Después de la ejecución del bucle, se puede dejar que el marco de optimización funcione automáticamente hasta que se logre cierto nivel de refinamiento. El criterio de convergencia para esto está abierto al investigador, y en la presente invención se ha controlado manualmente. Tras lograr un modelo sustituto suficientemente bueno, se usa el algoritmo de optimización para evaluarlo y encontrar un frente de Pareto para el problema multiobjetivo. Este marco MLADO podría extenderse a cualquier proceso de diseño de ingeniería sujeto a clasificación. Es interesante señalar que, en la presente invención, para cada variación geométrica, se debe generar una nueva malla. Por tanto, una característica interesante a agregar al marco MLADO puede ser incluir una opción de transformación de malla (mesh morphing). En este escenario, se podría definir un modelo de orden reducido del cálculo y la solución completa rápidamente inspeccionada en tiempo real. Tener el flujo representado en la misma topología de la malla también podría ayudar a definir métodos automáticos de detección de ciertas características del flujo.

Modelo de Clasificación para la Predicción de Desprendimiento de Vórtices

Es interesante poder predecir la presencia de desprendimiento de vórtices en ciertas aplicaciones debido a los efectos deseables/indeseables en el rendimiento de un sistema. En la presente invención, se desea el desprendimiento de vórtices para lograr una mezcla eficaz de la correspondiente magnitud escalar Y.

Con el fin de clasificar qué configuración puede (ES = 1) o no (ES = 0) conducir al desprendimiento de vórtices, se han desarrollado y probado dos modelos predictivos en un conjunto de datos de 80 simulaciones del pilar confinado en el canal. El conjunto de datos de 80 simulaciones es el más completo utilizado en este trabajo, pero es computacionalmente un poco costoso de crear. Sin embargo, como para comparar con la optimización tradicional también se desarrollarán 80 simulaciones, las 80 simulaciones completas ya están disponibles y se decide examinar y probar la construcción de modelos predictivos para este conjunto de datos completo ya disponible como referencia. Estos consisten en casos con sus parámetros de configuración en los rangos 120 < Re < 200, 0.2 < BR < 0.5 y 0.125 < AR < 1, mientras que, como se dijo anteriormente, rc = 104. La razón de tener un valor tan alto es tener un caso de mezcla muy desfavorable, donde la difusión de la magnitud escalar entre los fluidos es muy pobre y, a menos que se usen mecanismos adicionales, la eficiencia de mezcla también sería muy pobre. Esto ayuda a identificar cómo el desprendimiento de vórtices y las oscilaciones de estela aguas abajo mejora la mezcla de fluidos. Entre las 80 configuraciones simuladas, 27 eran desprendimientos sin vórtice (33,75%), por lo que el conjunto de datos no presentaba un desequilibrio preocupante. La muestra completa se ha dividido en dos partes (de manera estratificada, asegurando un % similar de casos VS = 0 y VS = 1 en ambas muestras): datos de entrenamiento (80% de la muestra, 64 casos) y datos de prueba (20% de la muestra, 16 casos). Debe destacarse que la determinación de si cada caso corresponde a VS = 1 o VS = 0 se basa en el coeficiente de oscilación de sustentación. Aunque esto podría detectarse automáticamente, esta asignación se ha realizado en una etapa de post-procesamiento por simplicidad, ya que el análisis no se realiza en tiempo real.

El primer modelo probado fue una Regresión Logística (LR). Este modelo es sencillo y suele tener buena capacidad predictiva. El modelo predictivo se ha entrenado con los datos de entrenamiento, obteniendo la LR que se muestra en la Tabla 2. Se puede observar que todas las variables incluidas en el modelo son significativas (es decir, son estadísticamente relevantes en el modelo). Si alguna variable tiene un valor p-valor superior a 0.05, se debe considerar si se descarta dicha variable del modelo.

Tabla 2: Modelo de regresión logística para la predicción de desprendimiento de vórtice.

El modelo LR debe validarse con los datos de prueba (denominados comúnmente como test), para garantizar la independencia en la construcción y evaluación. Esto se hace mediante el análisis de la llamada Curva de Operación de Recepción (ROC, del inglés Receiver Operating Characterístic). El área delimitada por esta gráfica se conoce como Área Bajo Curva (AUC, del inglés Area Under Curve) en la literatura y se usa como medida de precisión en los modelos de clasificación. El valor AUC = 1 significa que todos los casos están clasificados correctamente. Un valor de AUC = 0.5 significa que el modelo funciona como una decisión aleatoria y, por lo tanto, el modelo no es fiable para clasificar. En la Fig. 5 (a) se muestra la curva ROC relativa a los datos de test. El modelo LR logra un AUC = 1, pudiendo predecir con precisión todas las configuraciones de desprendimiento de vórtices en los datos de test. A pesar de que no es recomendable, el poder predictivo también se probó con los mismos datos utilizados para entrenar el modelo (datos de entrenamiento), con un AUC = 0.9859 (como se ve en la Fig. 5 (b). Una prueba adicional para demostrar la independencia y que no hay sobreajuste es dividir la muestra y realizar una validación cruzada de 5 partes (5-fold cross-validation). Esta prueba consiste en dividir los datos en 5 partes y tras ello el modelo se entrena con 4/5 partes de los datos y se valida con el 1/5 restante. Este proceso se repite alternando las 5 particiones, obteniendo así 5 modelos diferentes evaluados en 5 datos de prueba diferentes. La división de las muestras se ha realizado con estratificación (es decir, tratando de preservar un porcentaje de casos equilibrado entre VS = 0 y VS = 1). Los resultados se muestran en la Tabla 3, donde se puede ver una buena consistencia en un alto poder predictivo, ya que el modelo falla en predecir muy pocos casos.

Table 3: Validación cruzada con 5 partes del modelo de regresión logística.

Aunque el modelo de regresión logística desarrollado es muy preciso y la interpretabilidad del método es una ventaja, existen métodos de clasificación más sofisticados. Ejemplos de clasificadores populares y de alta precisión son Bosque Aleatorio (RF), Redes Neuronales en Aprendizaje Profundo (Deep Learning) o métodos de boosting como el Adapting Boosting (AdaBoost), Gradient Boosting Machine (GBM) o XGBoost. Entre estos, se ha seleccionado un algoritmo RF para crear un modelo predictivo para el desprendimiento de vórtice. Este método es muy popular en la literatura científica y es bien conocido que proporciona alta precisión y robustez, como se mostrará a continuación. Por tanto, no es necesario intentar otros algoritmos de mayor complejidad en este trabajo. Sin embargo, pueden ser considerados para otras aplicaciones.

Para una descripción teórica, se define como 0 el vector de hiperparámetros k que representan el algoritmo de RF. Cada árbol Jk se genera utilizando cada vector aleatorio 0 k y el conjunto de datos de entrenamiento Dt , lo que conduce a A (X) = J(X,Qk). Un clasificador de RF consta de una colección de clasificadores estructurados en árbol { ] (X ,0 k), k = 1 ,..., ntrees], donde los (0 k} son vectores aleatorios independientes distribuidos de manera idéntica y cada árbol determina un valor de predicción para entrada Xt , donde ntrees es el número de árboles que construyen el predictor RF. En los problemas de clasificación, RF es, por tanto, un conjunto de varios modelos y la clase CRF predicha a partir de los ntrees árboles de decisión se selecciona como el voto mayoritario de una clase c, es decir:

ntrees

Dado que RF es un conjunto de varios árboles de decisión, el algoritmo requiere definir el número de árboles ntrees. No existe una limitación teórica para esto, pero con un número de árboles elevado, los recursos computacionales se incrementarían innecesariamente. Se puede obtener una pista del número óptimo de árboles a partir del error Out-Of-Bag (OOB). Este parámetro es una medida de error que resulta de la evaluación de los árboles de muestras de datos de entrenamiento que no se incluyeron en el bootstrapping de los árboles, por lo que esto se hace para cada iteración y árbol consiguiente. Por encima de un cierto número de árboles, el OOB no se ve afectado. Otro parámetro ajustable en el Bosque Aleatorio es el número de variables que se muestrean aleatoriamente en cada división, m try. Esto es particularmente relevante para un conjunto de datos con una gran cantidad de variables. Debido a que en la presente invención solo hay 3 variables, el valor de m try se ha fijado en 3. El error OOB con respecto al número de árboles se muestra en la Fig. 6 para ntrees = 500. El error se estabiliza desde 400 árboles aproximadamente. Por lo tanto, se mantiene un total de 500 árboles, ya que se ha observado que la diferencia en el tiempo de cálculo era casi nula. Se ha permitido el muestreo con reemplazo para el algoritmo bootstrap del Bosque Aleatorio, ya que los datos de entrenamiento tenían un tamaño muy limitado para entrenar todos los árboles de decisión subyacentes. El predictor Bosque Aleatorio tiene un poder predictivo muy alto, sin falsos positivos ni falsos negativos para la validación de los datos de test. Estos resultados se muestran en las matrices de confusión de la Tabla 4. La matriz de confusión muestra la salida del predictor (valor predicho) versus el valor real (referencia). Esto corresponde a un AUC = 1, ya que todos los casos se predicen correctamente. Debe tenerse en cuenta que el conjunto de datos no se divide en datos de test y de entrenamiento, ya que el algoritmo de RF ya está dividiendo los datos por dentro y probando varios modelos de árbol de decisión. Por lo tanto, no tiene sentido dividir más aún los datos para hacer pruebas. No obstante, aunque no resulte muy relevante, como se hizo para el modelo LR, se desarrolla una validación cruzada de 5 partes para el predictor de RF, cuyos resultados se presentan en la Tabla 5. Durante la validación cruzada se puede observar que la precisión no es AUC = 1 para algunas particiones de muestra, pero aún así sigue siendo muy alto. El predictor RF ha funcionado de manera similar al predictor LR, por lo que se puede concluir que ambos se recomiendan en esta aplicación en particular, aunque en otros problemas de predicción más compleja sus rendimientos pueden diferir. En nuestra opinión, se prefiere usar RF porque se pueden encontrar varias aplicaciones en la literatura donde supera a LR debido a que es más sofisticado.

Tabla 4: Matriz de confusión para el predictor Bosque Aleatorio.

Tabla 5: Validación cruzada de 5 particiones del predictor Bosque Aleatorio.

Tanto los modelos LR como RF antes mencionados se basan en un total de 80 simulaciones, pero para seleccionar este número de simulaciones como referencia se desarrolló un estudio previo. Dado que en la presente invención hemos elegido el modelo de RF como definitivo, el análisis del número de muestras necesarias (similar a un análisis de convergencia) se muestra solo para este modelo. Para el estudio, se requiere un conjunto de datos inicial "suficientemente completo” para entrenar el predictor, que hemos establecido en un DoE de Np = 36 puntos casi equiespaciados como la combinatoria de fíe = (120,160,200], BR = {0.2,0.3,0.4,0.5} y AR = (0.125,0.5,1). El RF se entrena con este conjunto de datos y se utiliza para predecir siguientes niveles de refinamiento. En este estudio, se han desarrollado cuatro niveles de refinamiento a partir de las Np = 36 muestras: a) Agregación de fíe = 140 a la combinatoria (ahora Np = 36 para clasificar 48 casos), b) Agregación de ^ fí = 0.25 a la combinatoria de el refinamiento anterior (ahora Np = 48 para clasificar 64 casos), c) Agregación de fíe = 180 a la combinatoria del refinamiento anterior (ahora Np = 64 para clasificar 80 casos), y finalmente d) Agregación de ^ fí = 0.8 a la combinatoria del refinamiento anterior (ahora Np = 80 para clasificar 100 casos). La evolución de la precisión de la predicción se muestra en la Tabla 6, donde se puede observar que cuando el predictor se entrena con Np = 80 simulaciones, los datos de referencia de 100 simulaciones se predicen sin errores de clasificación.

Mediante la observación de la evolución de las predicciones, se puede decidir detenerse en un cierto número de simulaciones Np, elegido como número definitivo de simulaciones para entrenar al predictor. A partir de ese valor, Np se puede aumentar (es decir, el predictor de RF se puede actualizar dentro de MLADO) dependiendo del error cuadrático medio (MSE) o cualquier otra métrica de los sustitutos de Kriging, si así se desea. Debe recordarse que Np (número de muestras del predictor de RF) y Ns (número de muestras para el modelo sustituto) no necesitan coincidir, como se demostrará más adelante. Por ejemplo, el predictor se puede entrenar con simulaciones básicas (simulaciones estacionarias, mallas más gruesas, simulaciones de menor fidelidad, datos externos, etc.) ya que solo es necesario determinar el V5. Por otro lado, los modelos sustitutos de Kriging (o cualquier método sustituto/de interpolación utilizado) requieren la simulación transitoria completa y precisa, ya que el cálculo preciso de y n es imprescindible. Afortunadamente, este método MLADO ayuda en la decisión de simular o no estas costosas simulaciones.

De los resultados dados en la Tabla 6 se pueden extraer conclusiones útiles. El procedimiento univariado anidado (anidado quiere decir que se mantienen las muestras anteriores y se van añadiendo sobre estas) seguido en el análisis de la convergencia es solo un ejemplo de cómo proceder y se pueden considerar otras opciones, como el muestreo anidado aleatorio o las cuadrículas dispersas (sparse grids). Sin embargo, el problema bajo estudio no ha requerido de un análisis de convergencia más formal, ya que incluso con pocas muestras, el algoritmo RF ha proporcionado una capacidad predictiva sobresaliente. Incluso se ha comprobado que el RF entrenado con Np = 48 predijo las 100 muestras con solo 5 clasificaciones erróneas (dos falsos positivos y tres falsos negativos). Este desempeño sobresaliente se debe a que la activación del desprendimiento de vórtices no espera cambios repentinos, ya que es un problema de respuesta suave. Sin embargo, en otras aplicaciones, la situación puede no ser tan favorable. El presente análisis sugiere que para otras aplicaciones en la predicción de desprendimiento de vórtices, se puede predecir el rendimiento con un número relativamente bajo de muestras de datos de entrenamiento. Es importante señalar que las clasificaciones erróneas observadas en la Tabla 6 corresponden en realidad a valores en el límite de separación de VS = 0 y VS = 1. Para esta aplicación, estos problemas en la predicción no son motivo de preocupación. Primero, porque la línea que separa y agrupa VS = 0 y VS = 1 es difusa (¿debería considerarse una oscilación muy débil como VS = 0 o VS = 1?). Para ser conservadores, en esta demostración incluso un caso de prueba de oscilación débil se marca con VS = 1. Y segundo, porque los puntos con alta eficiencia de mezcla en teoría están lejos de esa región difusa (las zonas lejanas corresponden a desprendimiento de vórtice más fuerte), por lo que la clasificación errónea en las proximidades del límite de separación entre VS = 0 y VS = 1 no son una gran preocupación. Las clasificaciones erróneas se muestran en la Tabla 7, y en la Fig. 7 se puede observar que estos puntos están cerca del límite de separación.

Datas referencia

Valores predichos

a) b)

c) <1) Tabla 6: Evolución de la precisión del algoritmo RF en la predicción del siguiente nivel de refinamiento como muestras de test. a) Matriz de confusión en la clasificación de 48 muestras (RF entrenado con Np = 36, denominado RF1). b) Matriz de confusión en la clasificación de 64 muestras (RF entrenado con Np = 48, denominado RF2). c) Matriz de confusión en la clasificación de 80 muestras (RF entrenado con Np = 64, denominado RF3). d) Matriz de confusión en la clasificación de 100 muestras (RF entrenado con Np = 80, denominado RF4).

r CFD)

Tabla 7: Muestras mal clasificadas de la Tabla 6 y predicción de RF2 sobre 100 casos.

Aplicación de un Marco de Optimización de Diseño Asistido por Aprendizaje Automático al Diseño de un Micromezclador Mecánico

El primer paso consiste en construir los sustitutos de Kriging a partir de los datos. Como datos de entrenamiento hemos considerado inicialmente ambos casos con y sin desprendimiento de vórtices, debido a que esto ayuda a entender visualmente lo que sucede en los modelos surrogados. El cambio en los modelos sustitutos es insignificante en las áreas relevantes (donde la mezcla es importante). La visualización de ambos casos en la misma figura utilizando diferentes marcadores permite comprender visualmente regiones diferenciadas entre configuraciones con y sin desprendimiento de vórtices. Sin embargo, el error cuadrático medio y la predicción del RF determinan la necesidad real de refinar la distribución de puntos de datos.

En esta implementación (y para el método MLADO en general), no es necesario que coincidan el número de muestras utilizadas en el modelo sustituto y el predictor de desprendimiento de vórtices. Esto es estratégico y uno de los beneficios del método: el algoritmo de clasificación se puede entrenar sobre otros datos (simulaciones estacionarias, mallas gruesas, correlaciones empíricas, simulaciones de baja fidelidad, simulaciones laminares, datos externos, etc.). A partir de esos datos, se puede clasificar y decidir qué simular (simulación completa) a continuación, para crear los modelos sustitutos para la optimización. Por ejemplo, en la presente invención, el conjunto de simulaciones antes mencionado para entrenar el predictor de RF consistió en un DoE de puntos equidistantes ejecutados en FLUENT en estado estacionario. Esto fue suficiente para observar cuándo aparece la oscilación en las fuerzas de sustentación y luego marcar cada caso como VS = 0 o VS = 1. El predictor de RF entrenado con 80 simulaciones ha sido elegido en nuestra optimización del micromezclador, ya que tales datos están ya disponibles del estudio de optimización tradicional, por lo que no tiene sentido utilizar predictores de menor calidad. Sin embargo, un predictor con menos de 80 simulaciones sigue siendo preciso y útil para predecir el comportamiento oscilatorio. El marco MLADO puede ser muy interesante cuando el procedimiento de optimización comienza con un pequeño conjunto de datos o cuando hay una buena cantidad de datos disponible para construir sustitutos, pero el refinamiento adicional es costoso. La primera opción, que comienza con pocas muestras, es donde se puede aprovechar más el MLADO, ya que el problema se transforma en una solución eficiente basada en datos para guiar el proceso (data-driven) de refinamiento. Para demostrar la capacidad del método en el diseño del micromezclador, en esta sección, primero se presenta un proceso de optimización sin MLADO pero usando DoE casi equidistante para modelos sustitutos de Kriging. Esto se utilizará como referencia. En segundo lugar, el MLADO se aplicará desde una etapa de diseño inicial de Ns = 36 simulaciones. Estos resultados se compararán con los de referencia para cuantificar los recursos computacionales ahorrados.

Enfoque Tradicional para la Optimización de un Micromezclador

Es frecuente ver en la literatura el uso de modelos sustitutos de Kriging para la optimización en CFD, donde el DoE consiste en un número importante de costosas simulaciones en términos de recursos computacionales. Entonces, la decisión de simular una nueva y costosa muestra adicional depende únicamente del MSE.

Se ha desarrollado un estudio de optimización a partir de un conjunto de individuos DoE casi equidistante, que consiste en la combinatoria de Re = {120,140 160,180,200}, BR = {0.2, 0.3,0.4,0.5} y AR = {0.125,0,25,0,5,0,8,1}. Esto corresponde a una muestra total de 100 simulaciones. Sin embargo, debido a que 100 simulaciones es demasiado grande para nuestros recursos computacionales, algunas simulaciones van a ser descartadas. Por experiencia, se sabe que = 0.8 no es un candidato útil, por lo que solo se seleccionan 80 simulaciones en el DoE, descartando las combinaciones con = 0.8.

En la implementación actual, la toolbox de DACE [108] se utiliza para construir los sustitutos de Kriging. Estos se basan en las muestras de entrada. Para la interpolación es aconsejable normalizar los datos de entrada, pero esta toolbox realiza la normalización automáticamente. Para el método del Proceso Gaussiano (Kriging), se debe especificar la función de correlación para dos puntos en las ubicaciones x¿ y x@, C(0¿; x¿,x@). Esta función es importante para construir un buen estimador de Kriging. La función depende solo de la distancia entre dos puntos: cuanto menor es la distancia, mayor es la correlación. Esto ayuda a lidiar con agrupaciones (clusters) de datos. Por otro lado, si se aumenta la distancia entre los dos puntos, la correlación cae hasta cero. En la toolbox DACE hay varias opciones disponibles como función de correlación (aunque también se pueden construir funciones personalizadas). Las correlaciones probadas en este trabajo se muestran en la Tabla 8.

Tabla 8: Funciones de correlación testadas para los modelos sustitutos, definidas en [108].

En las Figs. 8 y 9 se pueden ver los modelos sustitutos, que se muestran para valores constantes de Re y AR, respectivamente. De manera análoga, en las Figs. 10 y 11, también se dan los errores cuadráticos medios (MSE). En estas figuras, los cuadrados rojos son los casos de datos simulados sin desprendimiento de vórtice, mientras que los puntos azules son los casos con desprendimiento de vórtice. Estos modelos sustitutos corresponden a un modelo de Kriging universal de regresión de segundo orden con correlación exponencial generalizada (GEXP). En las Figs. 12 y 13 también se muestran los sustitutos de r¡ por medio de otras correlaciones, solo para AR = 1 y solo para Re = 120, respectivamente, ya que estos son los sustitutos más complicados. Como puede verse en la Fig. 12, la forma ondulada de la función puede corregirse mediante una correlación lineal. Sin embargo, las pendientes repentinas en ciertos puntos producen un sustituto con áreas afiladas poco realistas. El uso de splines con Kriging Ordinario proporciona buenos sustitutos si se observa para AR = 1 en la misma figura, pero la realidad es que no tiene un mejor desempeño que el GEXP, como se observa con el número de Reynolds fijado en Re = 120 en la Fig. 13. También se observó que el MSE es considerablemente mayor que los valores mostrados en las Figs. 10 y 11. Las correlaciones esféricas funcionaron mejor con Kriging Universal con regresión de segundo orden, pero produjeron modelos sustitutos afilados similares a los de la correlación lineal. La correlación gaussiana tuvo un desempeño similar al GEXP, ya que en realidad es una exponencial generalizada de segundo orden, pero la ondulación es más notoria que para GEXP, además de exhibir un MSE más grande. Por lo tanto, se encontró que la correlación más apropiada era el GEXP, aunque se nota cierta ondulación en la vecindad de t] = 0. No obstante, dado que los valores óptimos de eficiencia de mezcla deben estar lejos de r¡ = 0, esto no es relevante.

De manera similar, en las Figs. 14-17 se muestran los sustitutos y MSE para la potencia de bombeo adimensional (n). Dado que estos puntos de datos de entrenamiento tienen una respuesta suave, el modelado se simplificó a los sustitutos de Kriging Ordinario con una correlación exponencial simple por cuestiones de simplicidad.

Mediante los modelos sustitutos generados es posible encontrar una configuración óptima para la geometría del pilar en el canal. Para este estudio de optimización, existen dos funciones objetivo: la eficiencia de mezcla r¡ y la potencia de bombeo (n). El objetivo es maximizar r¡ y minimizar (n). Cuando se tratan estos objetivos por separado, restringidos por los rangos dados de AR, BR y Re en la presente invención, las mejores configuraciones se dan en la Tabla 9.

Estas configuraciones se determinan a partir de las configuraciones CFD simuladas. Como se ve en la tabla anterior, cuando r¡ es alto, (n ) también es alto y viceversa. Por lo tanto, los valores óptimos de compensación deben calcularse mediante un enfoque multiobjetivo. Para ello se utiliza el algoritmo NSGA-II comentado con anterioridad. Los valores candidatos óptimos del frente de Pareto se dan en la Fig. 18. Se puede ver que los valores del Re permanecen agrupados alrededor de Re = 200, por lo que en la práctica ese parámetro podría fijarse en 200 y, por lo tanto, el espacio de diseño se reduciría a BR y AR. De la figura, cualquiera de los puntos mostrados en verde con marcador circular es una posible solución óptima al problema. Entre estos, los puntos A, B y C son finalmente seleccionados como candidatos. Puede observarse que si se elige el punto A, se daría prioridad a lograr una buena eficiencia de mezcla a un mayor coste de potencia de bombeo. Por el contrario, si el Punto C es la opción definitiva, se estaría más interesado en el bombeo a baja potencia en cierto detrimento de la eficiencia de la mezcla. El punto B sería una solución intermedia. El grupo de puntos cerca de ^r= 0, como se mostrará más adelante, corresponde a casos sin desprendimiento de vórtices. Dado que las configuraciones de interés son, en principio, las relacionadas con el desprendimiento de vórtices, los modelos sustitutos podrían haberse desarrollado utilizando solo los datos de entrenamiento de las simulaciones con desprendimiento (ES = 1). Para demostrar cómo afectaría esto a los resultados finales, se ha probado el proceso de optimización con VS = 1 solamente. El frente de Pareto resultante se muestra en la Fig. 19, que es prácticamente idéntico al frente de Pareto en la Fig. 18 cuando > 10, por lo que el grupo de puntos óptimos cerca de = 0 parece ser debido a configuraciones sin desprendimiento de vórtices.

Optimización del Diseño Asistido por Aprendizaje Automático de un Micromezclador

El interés real del uso del método MLADO es ahorrar recursos computacionales. Para este objetivo, el marco se ha probado en un conjunto inicial de pocas muestras de datos que constan de tres puntos en las coordenadas del espacio de diseño AR y Re, y cuatro en la coordenada BR. Este es un total de Ns = 36 simulaciones de la combinatoria de Re = {120,160,200}, BR = {0.2,0.3,0.4,0.5} y AR = {0.125,0.5,1}. Tras esto, los modelos sustitutos se refinan de acuerdo al MSE y al predictor de desprendimiento de vórtices RF entrenado con Np = 80 (RF4), pero solo aquellos casos con ES = 1 se agregarán al conjunto de datos de población original de Ns = 36 individuos. Debe recordarse que cualquiera de los predictores de desprendimiento de vórtices RF1, RF2 y RF3 podría ser utilizado, tal como se mencionó anteriormente (estos clasifican erróneamente solo algunas configuraciones cerca del límite de separación crítico), y también que estos no necesitan cálculos completos para el entrenamiento. Por lo tanto, preferimos usar el RF4, ya que es el mejor entrenado y no tiene sentido usar un predictor menos preciso, aunque los frentes de Pareto que se obtienen son esencialmente los mismos.

Por lo tanto, siguiendo los pasos sintetizados en la Fig. 4, y las configuraciones para construir modelos sustitutos explicadas en secciones previas, el primer paso es construir el sustituto a partir de los datos iniciales Ns = 36 y decidir si refinar. Debido a que el DoE era de puntos (casi) equidistantes, las áreas de mayor MSE son aquellas donde el espacio de diseño está menos discretizado. Estas áreas son, por ejemplo, aquellas en las áreas de Re = 140, AR = 0.25, Re = 160 y AR = 0.8, ya que tanto el Re como el AR tenían solo 3 puntos en cada dirección y la respuesta es altamente no lineal para r¡. Estas direcciones se pueden refinar de manera univariante o multivariante. Dado que solo se considerará VS = 1 para refinar el espacio de diseño, un enfoque univariante puede ser suficiente para beneficiarse del MLADO. Sin embargo, especialmente para DoEs más sofisticados, un muestreo anidado puede ser más eficiente.

De acuerdo con el diagrama de la Fig. 4, existe un bucle de ejecución bajo ciertas condiciones, que se dejan abiertas a la aplicación concreta. La primera ejecución del bucle se ha configurado para incrementar los Ns = 36 iniciales con simulaciones combinatorias de Re = 140, pero solo se incluirán aquellas que conduzcan a VS = 1. Por lo tanto, la nueva muestra de simulaciones es Ns = 42. Aunque el bucle en la Fig. 4 sugiere que el experto es quien decide cuándo dejar de realizar el refinamiento del espacio de diseño bajo algunos criterios, uno puede, e.g. implementar los algoritmos de optimización en cada etapa de refinamiento y observar algún tipo de convergencia. Este es el proceso seguido en esta implementación.

Como segunda ejecución del bucle, el MSE sugiere refinar, por ejemplo, las posiciones AR = 0.25, debido a la gran separación entre AR = 0.125 y AR = 0.5. El predictor RF se utiliza para predecir cuáles son los casos combinatorios VS = 1 a agregar a las muestras anteriores del modelo sustituto, lo que da como resultado un número total de Ns = 55 muestras. Una tercera ejecución del bucle conduce a un refinamiento en la dirección Re = 180, nuevamente considerando solo el VS = 1 predicho. Esto conduce a un número total de Ns = 68 individuos. El problema se puede refinar aún más, pero como se observa en la Fig. 20, los resultados del procedimiento de optimización muestran cierta convergencia. Los problemas de optimización multiobjetivo son complejos de representar (los frentes de Pareto generalmente se construyen para ilustrar los óptimos), pero en la optimización de un solo objetivo, el candidato óptimo es un solo individuo, por lo que el procedimiento se puede adaptar, e.g. para comprobar la convergencia de un óptimo. Los métodos basados en gradientes también podrían usarse en MLADO para guiar mejor la colocación de nuevas muestras en regiones potenciales de candidatos óptimos.

En la Fig. 20 se puede observar que incluso con pocas muestras de Ns, el frente de Pareto está muy cerca de la referencia para r¡ < 35% y (n ) < 2.25. Sin embargo, por encima de esta región, se observa que es necesario algún refinamiento para el espacio de diseño de los modelos sustitutos Ns = 55 en comparación con la referencia de Ns = 80 puntos. Por lo tanto, el refinamiento en Re = 180 era recomendable (los sustitutos mejoran, ya que la discretización entre Re = 160 y Re = 200 era demasiado vasta y los óptimos se ubican en Reynolds altos). Se podría considerar un mayor refinamiento entre Re = 180 y Re = 200, pero el MSE no es grande. Por tanto, un total de Ns = 68 parece suficiente para construir modelos sustitutos fiables.

Cabe señalar que el proceso de optimización puede ser de hecho más eficiente. Cuando se considera el conjunto de datos inicial de Ns = 36, si se implementa el algoritmo de optimización, no se observan candidatos óptimos a números de Reynolds bajos tras analizar el frente de Pareto, por lo que el refinamiento en Re = 140 podría descartarse. Además, especialmente para espacios de diseño más grandes donde por intuición es complicado identificar las mejores regiones, el predictor RF se puede usar para explorar si el número de casos VS = 1 aumenta o disminuye cuando se varía el Re (o cualquier parámetro). En este caso, el número de casos con VS = 1 va disminuyendo a medida que nos acercamos a Re = 140. Dicho esto, de manera similar al proceso antes mencionado, los casos combinatorios de AR = 0.25 y Re = 180, junto con sus casos predichos como VS = 1 mediante el algoritmo RF, se pueden utilizar para refinar el conjunto de datos de los modelos sustitutos. Estos dos refinamientos conducen a un número total de solo Ns = 59 muestras, produciendo resultados idénticos a los del análisis anterior con Ns = 68 y similares a los resultados de referencia con Ns = 80 (ver Fig. 20).

Finalmente, se puede concluir que ciertos recursos computacionales se pueden ahorrar mediante el enfoque MLADO. Para cuantificar los costes totales en comparación con la optimización tradicional, se considerará el coste total de cada simulación. Las simulaciones que resultaron en flujos estacionarios (ES = 0) necesitaron aproximadamente 1800 iteraciones a una velocidad de 2 segundos por iteración para converger, por lo que el tiempo total transcurrido fue de aproximadamente 1 hora por simulación. Por otro lado, las simulaciones que exhibieron un flujo oscilatorio (ES = 1) fueron cálculos transitorios, que convergieron después de varios pasos de tiempo con aproximadamente 40E3 iteraciones a una velocidad de 2 segundos. Es decir, un tiempo total transcurrido de 22,2 horas por simulación. Los costes totales se comparan en la Tabla 10, donde Ns = 100 sería el coste total del modelo sustituto completo con Re = {120,140,160,180,200}, BR = {0.2,03,0.4,0.5} y

AR = {0.125,0.25,0.5,0.8,1}.

₁ N _{1 s} Simulaciones VS = 0 Simulaciones VS = 1 Total horas simulación

Tabla 10: Costes computacionales (medidos como tiempo de simulación) en función del número de simulaciones.

A estos costes hay que añadir los relacionados con el entrenamiento del algoritmo Random Forest. Sin embargo, en otros estudios con datos externos disponibles (correlaciones empíricas, conjuntos de datos relevantes de datos experimentales, etc.) este coste no está incluido (los recursos computacionales en el entrenamiento de un predictor de RF son insignificantes). Como se vio anteriormente, incluso un algoritmo entrenado con pocas muestras es preciso. Por ejemplo, el algoritmo entrenado con el conjunto de datos inicial de Np = 36 simulaciones completas (RF1) clasificó erróneamente solo 7 de las 100, y sería entrenado sin costos adicionales. El algoritmo entrenado con Np = 48 (RF2) clasificó erróneamente solo 5 muestras de 100, y solo necesitó 12 muestras adicionales obtenidas mediante simulaciones estacionarias parcialmente convergidas, ya que los casos VS = 1 muestran oscilaciones en simulaciones configuradas como estacionarias aún no siéndolo. Esto permite detectar el comportamiento sin necesidad de una convergencia completa en régimen transitorio. Dado que la detección de desprendimiento de vórtices en estas simulaciones requirió aproximadamente 2/3 del tiempo necesario para obtener una simulación estacionaria completamente convergida, el método es aún más eficiente que el enfoque tradicional. No es necesario utilizar un predictor entrenado con más de Np = 36, ya que también se ha visto anteriormente que las clasificaciones erróneas tienen lugar en el límite de separación, que no es una región de interés en términos de alta eficiencia de mezclado. Por consiguiente, en la Tabla 10 se puede observar que el MLADO con un conjunto final de Ns = 59 puede reducir el tiempo computacional hasta un 27% con respecto a una optimización estándar con Ns = 100, y un 18% con respecto a optimización estándar con Ns = 80. Sin embargo, destacamos que en la presente invención, las simulaciones a evitar no son de hecho las costosas (ES = 0 son simulaciones estacionarias). Para otras aplicaciones, es decir, diseños en los que se busca suprimir el desprendimiento de vórtice, el proceso de optimización estaría destinado a diseños con VS = 0. En esta situación, se recomienda especialmente el uso del MLADO, ya que las simulaciones a evitar serían aquellas con VS = 1, logrando ahorros computacionales más relevantes.

Como se ve en la Fig. 20, los puntos óptimos seleccionados A, B y C son los mismos con y sin el MLADO. Antes de proceder a simular los casos seleccionados, es útil comprobar mediante el modelo predictivo si se puede producir el desprendimiento de vórtices. Si no hay desprendimiento de vórtices, entonces el rendimiento de las configuraciones de baja eficiencia puede ser discutible. Las configuraciones de los tres puntos A, B y C se introdujeron en el predictor Bosque Aleatorio y este predijo desprendimiento de vórtices para todos. La descripción de los tres puntos se da en la Tabla 11 (que incluye también el rendimiento simulado con CFD), y su representación sobre los modelos sustitutos se muestra en los isocontornos en la Fig. 21, donde en aras de una buena visualización, se ha supuesto que para los puntos A, B y C el valor óptimo de Re es exactamente Re = 200. En la Tabla 11 se puede observar que los sustitutos de Kriging predijeron un desempeño muy cercano al escenario simulado. Este es un buen indicador de la precisión de los modelos sustitutos para explorar configuraciones no simuladas. Las visualizaciones relevantes de CFD de los puntos óptimos seleccionados se muestran en la Fig. 22.

Opt imo candidato Afí

predicho (II) predicho r/ CFD (íl) CFD Punto A 0.130 0.500 199.956 49.032 % 3.288 48.96 % 3.284 Punto B 0.131 0.370 199.875 30.305 % 1.875 30.2 % 1.85 Punto C 0.251 0.203 199.879 12.751 % 0.896 13.4 % 0.891

Tabla 11: Puntos óptimos del problema de optimización multiobjetivo. Los valores de ^ y (n) predichos corresponden a los valores proporcionados por la evaluación de los sustitutos. CFD ^ y (n) son los valores de la simulación CFD de los puntos óptimos candidatos.

Para concluir, en la Tabla 12 se muestra una comparación del desempeño con respecto a trabajos previos en la literatura (que se relacionan con diferentes geometrías). Debido a que los trabajos previos utilizados con fines comparativos son sobre la mezcla de fluidos con diferentes concentraciones de fracción de másica, rc sería el número de Schmidt: rc = Se. Por lo tanto, y para una comparación justa, tanto los números de Reynolds como los de Schmidt en los estudios que se van a comparar deben ser los mismos, pero cada autor ha utilizado unos valores diferentes. Sin embargo, la comparación aún puede ser valiosa para comprender la capacidad de diseño óptimo de la invención. Para respaldar el análisis, se ha realizado una nueva simulación con el número de Schmidt Se = 103, y con la configuración geométrica del Punto A (el candidato óptimo con la mayor eficiencia de mezcla y consumo de energía de bombeo). Si el presente estudio se compara con trabajos anteriores con números de Schmidt de 104 (ver Tabla 12), la eficiencia del diseño propuesto supera a todas las de estos dispositivos. En cuanto a la caída de presión (Ap se proporciona en los trabajos anteriores y se hace adimensional en este documento con pU2), se puede ver que también presenta una baja caída de presión para lograr una alta eficiencia. Aunque la caída de presión en [30] presenta valores más bajos, la eficiencia también es mucho menor. En realidad, si realmente se necesita reducir la Ap, los puntos B y C tienen únicamente valores de 1.85 y 0.891, respectivamente, y aun así superan en eficiencia a los diseños comparables. Los trabajos en la literatura relacionados con Se = 103, a pesar de lograr buenas capacidades de mezclado, reportan grandes valores de caída de presión. Es importante considerar la caída de presión (o la potencia de bombeo) en la comparación, porque cada diseño es diferente en términos base a sus elementos y configuración, por lo que centrarse solo en la mezcla no sería fiable, ya que los requerimientos de energía pueden ser grandes. Por lo tanto, en términos de relativos entre mezclado y bombeo, el nuevo microdispositivo aquí diseñado parece superarlos a todos, debido a que logra una buena eficiencia de mezclado con una potencia de bombeo muy baja. Sin embargo, para aclarar esta afirmación, la comparación debe contextualizarse aún más. Aunque una eficiencia de mezcla de r¡ ~ 50% puede parecer pequeña en comparación con otras eficiencias de mezcla en la literatura, el dispositivo óptimo propuesto en esta invención está diseñado para un fluido con un Peclet muy alto, un microcanal muy corto y con la interacción de un único obstáculo. Si se aumenta la longitud del microcanal y/o el número de obstáculos, la eficacia de la mezcla también aumentará notablemente, pero la caída de presión aumentará excesivamente, como se ve por ejemplo en [90, 112]. En estos trabajos cobra especial relevancia el fuerte impacto en la caída de presión por la agregación de numerosos obstáculos. Como agregar más obstáculos a nuestro diseño modificaría la mecánica de desprendimiento de vórtices y aumentaría notablemente la caída de presión, extender la longitud del microcanal sería una opción sencilla para aumentar la eficiencia de la mezcla. Para dar evidencia de esto, se ha realizado un cálculo adicional: se ha probado el impacto del aumento de la longitud del microcanal en una unidad de longitud para la configuración óptima del Punto A con Se = 104. Los resultados de esta longitud L = 6 se comparan con el diseño base óptimo de longitud L = 5. En la comparación se observa un aumento importante de r¡ = 49% a r¡ = 58.6%, con un aumento casi inapreciable de la caída de presión adimensional de 3,284 a 3,348. Así, el mezclador óptimo diseñado en este trabajo es una opción muy eficiente, especialmente cuando el microcanal debe ser corto y la caída de presión (potencia de bombeo) debe ser baja. Por lo tanto, si se extiende la longitud del microcanal, se espera que la eficiencia del micromezclador basado en desprendimiento de vórtices mejore notablemente bajo un coste de bombeo decente, ya que no se colocan objetos adicionales a lo largo del microcanal. Diseños futuros pueden estar orientados a desarrollar la presente invención como un problema de optimización multiobjetivo que incluya L como parámetro de diseño, ya que no debe ser excesivamente largo para evitar un aumento considerable de la caída de presión.

[90], 50, 103 [113], 0.29, 104 [112], 200, 700 [30], 200, 104 PI, ~200, 104 PI, ~200, 103

_{n [ % \} 13-90 ~90 12.8-17.8 ~49 ~52 Ap [-] 9.3-35.9

24.6-470.5

0.74-1.96

3.284

3.284

Tabla 12: Comparación de la eficiencia de mezcla obtenida por diferentes autores. Cada encabezado de columna significa: Referencia, Re, Se. PI son las siglas de Presente Invención (PI).

Referencias

[30] J Ortega-Casanova. On the onset of vortex shedding from 2D confined rectangular cylinders having different aspect ratios: Application to promote mixing fluids. Chemical Engineering and Processing-Process Intensification, 120:81-92, 2017.

[90] Sourav Sarkar, KK Singh, V Shankar, and KT Shenoy. Numerical simulation of mixing at 1-1 and 1-2 microfluidic junctions. Chemical Engineering and Processing: Process Intensification, 85:227-240, 2014.

[95] S Turki, H Abbassi, and S Ben Nasrallah. Effect of the blockage ratio on the flow in a channel with a built-in square cylinder. Computational Mechanics, 33(1):22-29, 2003.

[96] Pratish P Patil and Shaligram Tiwari. Effect of blockage ratio on wake transition for flow past square cylinder. Fluid Dynamics Research, 40(11-12):753, 2008.

[97] Atul Sharma and V Eswaran. Heat and fluid flow across a square cylinder in the twodimensional laminar flow regime. Numerical Heat Transfer, Part A: Applications, 45(3):247-269, 2004.

[108] S0 ren N Lophaven, Hans Bruun Nielsen, Jacob Sondergaard, and A Dace. DACE. A matlab kriging toolbox. Technical University of Denmark, Lyngby, Technical Report No. IMMTR-2002,12, 2002.

[112] Nita Solehati, Joonsoo Bae, and Agus P Sasmito. Numerical investigation of mixing per formance in microchannel T-junction with wavy structure. Computers & Fluids, 96:10-19, 2014.

[113] J Ortega-Casanova. Application of CFD on the optimization by response surface method- ology of a micromixing unit and its use as a chemical microreactor. Chemical Engineering and Processing: Process Intensification, 117:18-26, 2017.

Claims

REIVINDICACIONES

1. Método de optimización de diseño asistido por aprendizaje automático caracterizado por que comprende:

a. Crear un predictor usando datos previos disponibles o generados mediante simulaciones poco precisas o de bajo coste, y que permita identificar la aparición o no de un efecto o situación de interés;

b. Usar el predictor para identificar qué simulaciones (combinaciones de los parámetros de diseño) deben realizarse con mayor precisión, por corresponder a regiones del espacio de diseño en las que se produce con una mayor frecuencia la aparición o no del efecto o situación de interés;

c. Realizar simulaciones de mayor precisión para generar o refinar un conjunto de datos fiable para crear el al menos un modelo sustituto;

d. Crear el al menos un modelo sustituto relativo a la al menos una variable objetivo de interés; y

e. Determinar el al menos un diseño óptimo mediante la aplicación de un algoritmo de optimización sobre el modelo sustituto.

2. Método de optimización según la reivindicación anterior caracterizado por que comprende la adquisición o generación de los datos iniciales necesarios para crear el predictor empleando un procedimiento anidado de construcción de la muestra de entrenamiento y validación del predictor.

3. Método de optimización según cualquiera de las reivindicaciones anteriores caracterizado por que el predictor es un algoritmo de inteligencia artificial.

4. Método de optimización según la reivindicación anterior caracterizado por que el predictor se selecciona entre Regresión Logística (LR), Bosque Aleatorio (RF), Redes Neuronales en Aprendizaje Profundo (Deep Leaming) o métodos de boosting como AdaBoost (Adapting Boosting), GBM (Gradient Boosting Machine) o XGBoost (eXtreme Gradient Boosting).

5. Método de optimización según la reivindicación anterior caracterizado por que el predictor se selecciona entre Regresión Logística (LR) y Bosque Aleatorio (RF).

6. Método de optimización según cualquiera de las reivindicaciones anteriores caracterizado por que el algoritmo de optimización es un algoritmo de optimización metaheurístico.

7. Método de optimización según la reivindicación anterior caracterizado por que el algoritmo de optimización es un algoritmo de optimización evolutivo (por ejemplo, genético).

8. Método de optimización según la reivindicación anterior caracterizado por que el algoritmo de optimización empleado es el Algoritmo Genético de Ordenación No dominada II (NSGA-II).

9. Método de optimización según cualquiera de las reivindicaciones anteriores caracterizado por que el al menos un modelo sustituto creado debe ser refinado antes de poder aplicarse sobre él el algoritmo de optimización, y dicho refinamiento comprende un proceso de transformación de malla (mesh morphing) para probar diferentes geometrías rápidamente dentro del bucle de optimización, facilitar el análisis en tiempo real o posibilitar la definición de métodos automáticos de detección de ciertas características del modelo sustituto.

10. Micromezclador mecánico diseñado empleando un método de optimización conforme cualquiera de las reivindicaciones anteriores para, mediante desprendimiento de vórtices, conseguir la máxima eficiencia de la mezcla a la par que la mínima caída presión, dicho mezclador comprendiendo:

• Un pilar rectangular bidimensional que mide h y l metros de largo en la dirección x e y, respectivamente, dicho pilar ubicado en la línea central de un canal recto centrado en el origen de las coordenadas (x, y), estando la entrada del canal a una distancia Lu metros aguas arriba de la cara frontal del pilar;

• el canal tiene L metros de longitud, la entrada del canal consiste en un flujo laminar completamente desarrollado formado por dos fluidos idénticos con diferentes valores de la magnitud escalar Y a mezclar, coeficiente de difusividad escalar Dc, y con densidad y viscosidad conocidas, p y p, respectivamente;

• el flujo ingresa a la entrada del canal con una velocidad media U, por la mitad superior de la sección de entrada entra el fluido con valor nulo de la magnitud escalar Y, es decir, Y = 0 en x = - (Lu l / 2), 0 < y < H / 2, mientras que el otro fluido entra por la mitad inferior con un valor unidad de la magnitud escalar, es decir, Y = 1 en x = - (Lu l / 2), - H / 2 < y < 0; • los campos de velocidad v = (u ,v) y presión p del flujo 2D en la geometría descrita anteriormente están gobernados por las ecuaciones de conservación de masa y cantidad de movimiento; y

• la mezcla de Y se rige por la ecuación escalar de conveccióndifusión;

dicho micromezclador caracterizado por que presenta:

• Una relación de bloqueo BR para h comprendida en el rango 0,2 0,5;

• una relación de aspecto AR para l en el rango 0,125-0,5; y

• un número de Reynolds Re de 200.

11. Micromezclador mecánico según la reivindicación anterior diseñado para priorizar el incremento en la eficiencia de mezcla respecto de la potencia de bombeo adimensional / caída de presión, dicho micromezclador caracterizado por que presenta:

• una relación de bloqueo BR para h en el rango 0,37-0,5;

• una relación de aspecto AR para l de 0,125-0,131; y

• un número de Reynolds Re de 200.

12. Micromezclador mecánico según la reivindicación anterior caracterizado por que presenta:

• una relación de bloqueo BR para h en el rango 0,37;

• una relación de aspecto AR para l de 0,13; y

• un número de Reynolds Re de 200.

13. Micromezclador mecánico según la reivindicación 11 caracterizado porque presenta:

• una relación de bloqueo BR para h en el rango 0,5;

• una relación de aspecto AR para l de 0,13; y

• un número de Reynolds Re de 200.

14. Micromezclador mecánico según la reivindicación 10 diseñado para priorizar la reducción de la potencia de bombeo adimensional / caída de presión respecto de la eficiencia de mezcla, dicho micromezclador caracterizado por que presenta:

• una relación de bloqueo BR para h de 0,2;

• una relación de aspecto AR para l de 0,25; y

• un número de Reynolds Re de 200.