Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip

Equation Chapter 1 Section 1 Trabajo Fin de Máster Máster Universitario en Ingeniería Industrial Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip Autora: Eva Aguilera Mazuela Tutor: José Manuel Quero Reboul Dpto. Ingeniería Electrónica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 iii Trabajo Fin de Máster Máster Universitario en Ingeniería Industrial Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip Autora: Eva Aguilera Mazuela Tutor: José Manuel Quero Reboul Dpto. de Ingeniería Electrónica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 v Trabajo Fin de Máster: Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip Autor: Eva Aguilera Mazuela Tutor: José Manuel Quero Reboul El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2021 El Secretario del Tribunal vii A mi familia y amigos ix Resumen Para diversas aplicaciones se desea obtener una zona de un Lab-on-Chip a una temperatura constante y, para ello, se emplean pistas de cobre a modo de calentadores. Sin embargo, no es tarea fácil obtener una temperatura constante, ya que las zonas internas de los calentadores tienden a liberar menos calor y, por tanto, se encuentran a temperaturas más elevadas. Para lograr una temperatura constante se puede diseñar una pista de ancho variable, pero, a día de hoy, la decisión de qué anchos imponer y en qué zonas, solo puede realizarse de forma manual. Así que, para simplificar este trabajo, se desea elaborar un programa de cálculo capaz de definir los anchos de una pista de cobre para obtener una temperatura constante en una determinada zona. En este documento se explica, no solo el contenido de dicho programa, sino también los cálculos previos necesarios, el procedimiento a seguir para su correcto uso y algún ejemplo que muestre su funcionamiento de manera ilustrativa. xvii ÍNDICE DE FIGURAS Figura 1-1. Modelo en COMSOL de un ejemplo de LoC con un serpentín de cobre. 1 Figura 2-1. Representación de las resistencias térmicas laterales y verticales en una celda. 6 Figura 2-2. Perspectiva del modelo construido en COMSOL. Bloque de 7x7 celdas en el centro, rodeado de un bloque de 20mm de lado y atravesado por la pista de cobre en línea recta. 7 Figura 2-3. Planta del bloque central del modelo, donde están señalizadas las celdas en las que se ha medido cada temperatura para obtener las resistencias laterales. Siendo T1 y T2 usadas para el cálculo de la resistencia térmica entre dos celdas sin pista de cobre, T3 y T4 para el de dos celdas con pista de cobre, y T3 y T5 para el de una celda con pista de cobre y otra sin ella. 8 Figura 2-4. Ecuación del coeficiente entre una celda con cobre y otra sin cobre, en función de la relación entre el ancho de la pista de cobre y el lado de la celda. 11 Figura 2-5. Planta del bloque central del modelo, donde están señalizadas la distancias a una celda con cobre de cada celda en la que se han medido temperaturas para obtener las resistencias verticales. Siendo D=0 una celda que contiene pista de cobre, D=1 una celda contigua a una celda con pista de cobre y así sucesivamente. 11 Figura 2-6. Corte de temperatura obtenido en COMSOL, en el que se distinguen las variaciones de temperatura en als superficies en las que se han tomadado medidas para el cálculo de las resistencias térmicas verticales. 13 Figura 2-7. Diagrama de flujo que sigue el programa de cálculo para obtener la resistencia de ancho variable definitiva. 15 Figura 2-8. Imagen de una línea que define el recorrido de la pista de cobre, en la que se ha marcado, en rojo, el borde de la zona activa y, en azul, el de la zona ocupada por el plano de cobre. 17 Figura 2-9. Gráficas mostradas en Matlab para el seguimiento de parámetros. 27 Figura 2-10. Diagarama de flujo para el correcto uso del código en el programa de cálculo. 29 Figura 3-1. Planta obtenida de COMSOL del LoC del ejemplo en el que se desea colocar una pista de cobre. Los tres círculos que se muestran son las zonas en la que se desea obtener una temperatura de 65ºC y por encima de los cuales no puede pasar la pista de cobre. 31 Figura 3-2. Primer diseño de la línea que define el recorrido de la pista de cobre con el contorno de la zona activa marcado en rojo y el de la zona ocupada por el plano de cobre marcado en azul. Las tres zonas en las que no hay línea de cobre se corresponden con las zonas ocupadas por los tres círculos de la Figura 3-1. 32 Figura 3-3. Pista de cobre obtenida en L-Edit tras deshacer la discretización en el resultado obtenido en Matlab para una línea como la mostrada en la Figura 3-2. 32 Figura 3-4. Corte de temperatura en el plano de la pista de cobre, obtenido en COMSOL, para una pista de cobre como la mostrada en la Figura 3-3. 33 Figura 3-5. Distintos diseños de la línea que define el recorrido de la pista de cobre, numerados según el orden cronológico en el que fueron creados. 34 Figura 3-6. Mapa de temperatura del corte realizado en el plano de la pista de cobre, obtenido con las distintas líneas que definen la pista de cobre de la Figura 3-5. 35 Figura 3-7. Diseño de la línea que define el recorrido de la pista de cobre correspondiente a R5 en la Figura 3-5 y 3-6, con el contorno de la zona activa marcado en rojo y el de la zona ocupada por el plano de cobre marcado en azul. Zona activa definida tras comprobar que no se puede alcanzar la misma temperatura en los círculos que en la pista de cobre. 36 Figura 3-8. Pista de cobre obtenida en L-Edit tras deshacer la discretización en el resultado obtenido en Matlab para una línea correspondiente a R5 en la Figura 3-5 y 3-6. 36 Figura 3-9. Corte de temperatura en el plano de la pista de cobre, obtenido en COMSOL, para la pista de cobre de la Figura 3-8. 37 Figura 3-10. Pista de cobre de la Figura 3-8 fabricada en un PCB de doble cara. Los puntos azules se corresponden con los distintos puntos en los que se ha obtenido una medida de la temperatura. 38 Figura 3-11. Mapa de temperaturas de la superficie del metacrilato. La zona roja se corresponde con la de mayor temperatura, mientras que la azul se corresponde con la de menor temperatura. La temperatura media de la superficie es de 54.9ºC. 39 Figura 3-12. Puesto de trabajo en el laboratorio. 39 Figura AA-0-1. Recorrido de la pista de cobre diseñado en L-Edit. 43 Figura AA-0-2. Imágenes devueltas por leelinea.m tras convertir el fichero .DXF a una matriz. En la imagen de la izquierda se muestra una línea continua según el recorrido de la pista de cobre, y en la de la derecha se marcan los puntos discretos que conforman el recorrido. 44 Figura AA-0-3. Imagen del recorrido de la pista de cobre con el contorno de la zona activa marcado en rojo y el del plano de cobre marcado en azul. 45 Figura AA-0-4. Gráficas mostradas en Matlab para el seguimiento de parámetros. 46 Figura AA-0-5. Perfil de temperaturas de la línea media del LoC en una iteración. 47 Figura AA-0-6. Pista de cobre obtenido tras el cálculo de los anchos. A la izquierda se muestran los anchos de la pista mediante escala de grises, siendo el negro un ancho 0. A la derecha se muestran los anchos de manera discretizada, de manera que se ve la pista de cobre formada por pequeños cuadrados de distinto tamaño. 47 Figura AA-0-7. Pista de cobre obtenida en Matlab importada en L-Edit. Se encuentra seleccionada la línea completa y marcado el icono de “merge”. 48 Figura AA-0-8. Mapa de temperaturas en el plano de la pista de cobre obtenido en COMSOL. A la izquierda con el rango completo de temperaturas del LoC. A la derecha con el rango de temperaturas permitido, es decir Tobj±1ºC. 49 1 1 INTRODUCCIÓN AL PROBLEMA e quiere emplear una línea de cobre como resistencia, para obtener una temperatura constante en una determinada zona de un Lab-on-Chip (LoC). Para ello es necesaria una definición detallada del ancho de la línea de cobre. Actualmente se utiliza un simulador basado en el método de los elementos finitos, pero este no soporta la simulación de una geometría variable, por lo que se necesita un proceso iterativo manual para la modificación del ancho de la pista de cobre a lo largo de su recorrido. Es decir, se puede establecer una forma concreta para la resistencia con un ancho de cobre variable, realizar una simulación en COMSOL con dicha resistencia y a partir del resultado obtener conclusiones. Pero esas conclusiones se verán muy limitadas, ya que se puede intuir en qué zona puede ser necesario un aumento o disminución del ancho de la pista de cobre, pero esas zonas son muy amplias y un cambio en ellas puede cambiar el resultado por completo, alejándose del objetivo deseado. Incluso si se hicieran cambios muy localizados, estos cambios también afectan a las zonas cercanas, debido a la transmisión de calor, y dicho cambio no es fácil de intuir al hacerlo de forma manual para intetar compensarlo si fuera necesario. Por todo ello, este proceso iterativo manual se vuelve muy largo y tedioso, por lo que se hace prácticamente imposible llegar al objetivo a través de él. Y para agilizar este proceso, se pretende realizar un programa que permita calcular de forma automática los anchos de la pista de cobre necesarios. Figura 1-1. Modelo en COMSOL de un ejemplo de LoC con un serpentín de cobre. Un posible ejemplo de lo que podría ser un LoC con una pista de cobre que hace de calentador se muestra en la Figura 1-1. En ella se distinguen las capas de PCB (inferior) y metacrilato (superior), así como la figura serpenteante que forma la pista de cobre. En este ejemplo se ve un canal que pasa por encima del calentador, el cual serviría para llevar un líquido que se desea poner a una determinada temperatura. Esta podría ser una de las aplicaciones para las que se desea tener una zona a una temperatura constante. S 3 2 DESCRIPCIÓN DE LA SOLUCIÓN omo se ha dicho anteriormente, se necesita agilizar el proceso de definición de los anchos de la pista de cobre y, para ello, se pretende realizar un programa de cálculo que permita calcular los anchos necesarios de forma automática. Dicho programa debe ser capaz:  de leer el diseño de resistencia hecho previamente, en el que se define la forma y recorrido de la línea de cobre,  de estimar matemáticamente la generación y transferencia de calor,  de modificar de manera localizada el ancho de la pista de cobre, así como realizar cambios en la corriente para lograr la temperatura deseada,  de tener en cuenta las limitaciones tecnológicas y  de generar la geometría de la línea de resistencia para que pueda ser posteriormente leída desde la herramienta CAD. Estas funciones se explican con mayor grado de detalle en el apartado 2.2. Además, al realizar distintas simulaciones con dicho programa, se ha visto que es importante para lograr el objetivo deseado, colocar un plano de cobre en la zona ocupada por la resistencia, ya que esto contribuye en gran medida a la homogeneización de la temperatura. Este efecto es debido a que el cobre es el material, de los empleados, con mayor conductividad eléctrica, lo cual hace que tenga una resistencia térmica menor, por lo que en el intercambio entre una zona con cobre y otra sin cobre, hay una gran diferencia en el valor de las resistencias, siendo la zona con cobre la dominante en el intercambio de calor. Sin embargo, al colocar el plano de cobre, ya no existe esa gran diferencia en las resistencias, ya que toda la zona contiene cobre, por lo que no hay unas zonas dominantes frente a otras y el calor se distribuye de forma más homogénea. 2.1 Transmisión de calor Para poder elaborar el código que genere automáticamente una resistencia de ancho variable, es necesario establecer las ecuaciones que definen la transmisión de calor y el cálculo de resistencias térmicas. El modelo de transferencia de calor en sólidos por conducción es [1]: 𝜌𝐶𝑝𝜕𝑇 𝜕𝑡 +∇𝑞 =𝑄 (2–1) donde ρ es la densidad, Cp es la capacidad calorífica específica, T es la temperatura, t es el tiempo, q es el flujo de calor por conducción y Q es la fuente de calor. Simplificando para un sistema en el que no hay variación de temperatura con el tiempo, la ecuación de transferencia de calor en sólidos por conducción es: 𝑞 =−𝑘∇𝑇 (2–2) C Descripción de la solución 4 4 donde q es el flujo de calor por conducción, k es la conductividad térmica y T es temperatura. El modelo de transferencia de calor por convección es: 𝑞 =ℎ(𝑇𝑠−𝑇𝑎𝑚𝑏) (2–3) donde q es el flujo de calor por convección, h es el coeficiente de película convectivo, Ts es la temperatura en la superficie del sólido y Tamb es la temperatura ambiente. Discretizando el LoC y considerando que este está formado por celdas de lado L, cuyas superficies se encuentran a temperatura constante, las ecuaciones que se obtienen para una celda son las siguientes:  Generación de calor por efecto Joule [1] 𝑄𝑐𝑢 =𝑅𝑒_𝑐𝑢𝐼2 (2–4) donde Qcu es el calor generado, I es la corriente que circula por la pista de cobre y 𝑅𝑒_𝑐𝑢 = 𝜌𝐿 𝐻𝑐𝑢𝑊𝑐𝑢 es la resistencia eléctrica del cobre [2], donde ρ es la resistividad del cobre, L es la longitud del lado de la celda, Hcu es el espesor del cobre y Wcu es el ancho de la pista de cobre.  Transmisión de calor vertical, suponiendo que las temperaturas en las superficies de cada celda son constantes [1] 𝑇𝑐𝑢 −𝑇𝑡=𝑄𝑢𝑝𝑅𝑚𝑒𝑡 𝑇𝑐𝑢 −𝑇𝑏=𝑄𝑑𝑤𝑅𝑝𝑐𝑏 (2–5) donde Tcu es la temperatura en el centro de la superficie de la celda que separa el metacrilato y el PCB, Tt es la temperatura en el centro de la superficie superior de la celda, Tb es la temperatura en el centro de la superficie inferior de la celda, Qup es el calor disipado hacia arriba, Qdw es el calor disipado hacia abajo, 𝑅𝑚𝑒𝑡 =𝐻𝑚𝑒𝑡 𝐾𝑚𝑒𝑡𝐿2 es la resistencia térmica del metacrilato y 𝑅𝑝𝑐𝑏 =𝐻𝑝𝑐𝑏 𝐾𝑝𝑐𝑏𝐿2 es la resistencia térmica del PCB, 5 5 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip donde L es la longitud del lado de la celda, Hmet es el espesor del metacrilato, Hpcb es el espesor del PCB, Kmet es la conductividad térmica del metacrilato y Kpcb es la conductividad térmica del PCB.  Transmisión de calor lateral [1] 𝑇𝑐𝑢 −𝑇𝑐𝑢(±1,±1)=𝑄𝑙𝑎𝑡𝑅𝑙𝑎𝑡 (2–6) donde Tcu es la temperatura en el centro de la superficie de la celda que separa el metacrilato y el PCB, Tcu(±1, ±1) es la temperatura en el centro de la superficie que separa el metacrilato del PCB en cada una de las celdas contiguas a la primera, Qlat es el calor disipado hacia los laterales y Rlat es la resistencia térmica entre dos celdas contiguas, que puede calcularse de tres formas: 𝑅𝑙𝑎𝑡_𝑠𝑖𝑛𝑐𝑢 =𝐿 𝐾𝑝𝑐𝑏𝐻𝑝𝑐𝑏𝐿+𝐾𝑚𝑒𝑡𝐻𝑚𝑒𝑡𝐿 entre dos celdas sin pista de cobre, 𝑅𝑙𝑎𝑡_𝑐𝑜𝑛𝑐𝑢 =𝐿 𝐾𝑝𝑐𝑏𝐻𝑝𝑐𝑏𝐿+𝐾𝑚𝑒𝑡𝐻𝑚𝑒𝑡𝐿+𝐾𝑐𝑢𝐻𝑐𝑢𝑊𝑐𝑢 entre dos celdas con pista de cobre y 𝑅𝑙𝑎𝑡_𝑐𝑜𝑛𝑐𝑢_𝑠𝑖𝑛𝑐𝑢 =𝑊𝑐𝑢/2 𝐿(𝐾𝑚𝑒𝑡(𝐻𝑚𝑒𝑡−𝐻𝑐𝑢)+𝐾𝑐𝑢𝐻𝑐𝑢+𝐾𝑝𝑐𝑏𝐻𝑝𝑐𝑏)+𝐿−𝑊𝑐𝑢/2 𝐿(𝐾𝑚𝑒𝑡𝐻𝑚𝑒𝑡+𝐾𝑝𝑐𝑏𝐻𝑝𝑐𝑏) entre una celda con pista de cobre y otra sin ella, donde L es la longitud del lado de la celda, Wcu es el ancho de la pista de cobre, Hmet es el espesor del metacrilato, Hpcb es el espesor del PCB, Hcu es el espesor del cobre, Kmet es la conductividad térmica del metacrilato, Kpcb es la conductividad térmica del PCB y Kcu es la conductividad térmica del cobre.  Disipación de calor vertical [1] 𝑇𝑡−𝑇𝑎𝑚𝑏 =𝑄𝑢𝑝𝑅ℎ 𝑇𝑏−𝑇𝑎𝑚𝑏 =𝑄𝑑𝑤𝑅ℎ (2–7) donde Tt es la temperatura en el centro de la superficie superior de la celda, Tb es la temperatura en el centro de la superficie inferior de la celda, Tamb es la temperatura ambiente, Qup es el calor disipado hacia arriba, Qdw es el calor disipado hacia abajo y 𝑅ℎ=1 ℎ𝐿2 es la resistencia térmica convectiva, donde L es la longitud del lado de la celda y Descripción de la solución 6 6 h es el coeficiente de película convectivo. Haciendo balance de energía con las ecuaciones anteriores se tiene que 𝑄𝑢𝑝 +𝑄𝑑𝑤 +∑𝑄𝑙𝑎𝑡 4 1=𝑄𝑐𝑢 (2–8) de donde se puede obtener la temperatura de cada celda (Tcu(i, j)) que, al compararla con la temperatura deseada, permite tomar decisiones en el cambio del ancho de la pista y en la corriente a aplicar. Figura 2-1. Representación de las resistencias térmicas laterales y verticales en una celda. 7 7 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip 2.1.1 Cálculo de coeficientes La obtención de las resistencias térmicas en el apartado anterior se ha hecho de manera teórica, pero en la realidad la transmisión de calor no es homogénea, sino que se deforma a medida que se aleja del foco de emisión de calor (la pista de cobre). Para no tener en cuenta este efecto en el cálculo teórico, pero que el resultado se asemeje a la realidad, se aplican unos coeficientes sobre las resistencias térmicas obtenidas de manera teórica. Para obtener dichos coeficientes se han realizado unas simulaciones en COMSOL. Se ha tomado un sistema de 7x7 celdas de lado L=1mm, donde cada celda está compuesta de una capa superior de metacrilato de altura Hmet = 5mm, cuya conductividad térmica es Kmet = 0.19W/mK, y una capa inferior de PCB de altura Hpcb = 1,6 mm, cuya conductividad térmica es Kpcb = 0.3W/mK. Este sistema es más pequeño que el real, e incluso las celdas más internas estarían muy cerca del exterior y, por tanto, se verían afectadas por las condiciones del ambiente. Así que para alejarlas, y que el resultado se parezca más a lo que realmente se va a implementar, se ha añadido un bloque de 20mm de lado que rodea al anteriormente descrito de lado 7mm, con la misma estructura respecto a las capas de metacrilato y PCB que este. Por último, se añade una pista de cobre de espesor Hcu = 18µm , cuya conductividad térmica es Kcu = 400W/mK, entre la capa de metacrilato y la de PCB por la zona central del bloque. Se hace pasar por dicha pista una corriente de 1.5A (tras comprobar previamente que variar la corriente no afecta al resultado) y se va modificando su ancho desde 200µm hasta 600μm. No se han hecho pruebas con anchos mayores a este, porque al acercarse la pista a la zona en la que se toman las medidas (pared que separa una celda de otra) se distorsionan los resultados. Figura 2-2. Perspectiva del modelo construido en COMSOL. Bloque de 7x7 celdas en el centro, rodeado de un bloque de 20mm de lado y atravesado por la pista de cobre en línea recta. Se han obtenido coeficientes tanto para las resistencias laterales, como para las verticales. Y en el primer caso, se ha diferenciado si la transmisión de calor era entre dos celdas sin pista de cobre, dos celdas con pista de cobre o entre una de cada tipo. Descripción de la solución 14 14 2.2 Descripción del programa de cálculo En este apartado se trata el contenido del programa de cálculo que, tras leer la forma de la pista de cobre es capaz de modificar el ancho de esta para lograr una temperatura objetivo en una zona determinada del LoC. Se explica el proceso que sigue, así como las rutinas que lo forman. En la Figura 2-6 se muestra el diagrama de flujo que sigue el programa de cálculo para obtener la resistencia de ancho variable definitiva. Para entender mejor este diagrama de flujo se explican a continuación de manera detallada cada uno de los pasos a seguir. Paso 1. En primer lugar se lee el fichero .DXF elaborado en un programa de diseño, que contiene la línea que marca el recorrido de la pista de cobre. Una vez realizado esto ir al paso 2. Paso 2. A continuación, se transforma el fichero leído en el apartado anterior en una matriz que contiene el valor del ancho en cada punto de la pista de cobre, para que Matlab pueda operar con ella. Tras esto ir al paso 3. Paso 3. En este paso se definen las constantes que son necesarias en el cálculo, así como los valores iniciales de los parámetros que lo requieren. Estos valores pueden ser definidos directamente en el código principal, u obtenerse de otro fichero. Entre los valores de las constantes se encuentra el cálculo de las resistencias que no dependen del ancho de la pista de cobre. Otro valor relevante para seguir el diagrama de flujo, que se define en este punto, es el índice “Pasos”, que se inicializa a 0 e indica el número de iteraciones realizadas. Una vez definido esto, ir al paso 4. Paso 4. Se calculan las resistencias restantes, las que dependen del ancho de la pista de cobre. Y, tras ello, se va al paso 5. Paso 5. Se elaboran las matrices de conductancias térmicas, 𝐺 = 1 𝑅, y de fuentes (o sumideros) de energía, Q. Y se va al paso 6. Paso 6. Obtenidas las matrices G y Q se procede al cálculo de las temperaturas en cada punto del LoC, mediante la ecuación 𝑄 =𝐺𝑇, y se va al paso 7. Paso 7. Si en una celda Tobj – Tcu > máximo error local permitido ir al paso 8, si no ir al paso 11. Paso 8. Se disminuye el ancho de la pista de cobre en la celda en cuestión y se va al paso 9. Paso 9. Si, tras diminuir el ancho de la pista de cobre en una celda, este es menor que el mínimo permitido ir al paso 10, si no ir al paso15. Paso 10. Se iguala el ancho que es menor que el ancho mínimo permitido a este y se aumenta la corriente. Tras esto, ir al paso 15. 15 15 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip Figura 2-7. Diagrama de flujo que sigue el programa de cálculo para obtener la resistencia de ancho variable definitiva. Descripción de la solución 16 16 Paso 11. Si en una celda Tcu – Tobj > máximo error local permitido ir al paso 12, si no ir al paso 16. Paso 12. Se aumenta el ancho de la pista de cobre en la celda en cuestión y se va al paso 13. Paso 13. Si, tras aumentar el ancho de la pista de cobre en una celda, este es mayor que el máximo permitido ir al paso 14, si no ir al paso 15. Paso 14. Se iguala el ancho que es mayor que el ancho máximo permitido a este y se disminuye la corriente. Tras esto, ir al paso 15. Paso 15. Se compara el error medio obtenido en una iteración con la anterior y, en base a ello, se aumenta o disminuye el paso de modificación del ancho de la pista de cobre. Tras esto, ir al paso 17. Paso 16. Si el error medio obtenido es mayor que el error medio máximo permitido ir al paso 17, si no ir al paso 19. Paso 17. Se aumenta el índice de iteración, Pasos, en una unidad y se va al paso 18. Paso 18. Si Pasos < Pasos máximo permitido ir al paso 4, si no ir al paso 19. Paso 19. Bien porque se ha convergido y se ha obtenido la resistencia de anchos variables necesaria para obtener el objetivo, o bien porque se ha llegado al límite de iteración establecido, se finaliza el proceso de cálculo y se escribe un fichero .DXF que contiene la última figura de la pista de cobre obtenida con los anchos calculados. En los siguientes subapartados, se explica y se muestra parte del código elaborado en Matlab. 2.2.1 Constantes En el fichero MisCons.m se encuentran algunas de las constantes necesarias, que serán llamadas desde otros de los ficheros Matlab. Tales como:  la temperatura objetivo (Tcu),  las dimensiones del LoC (DimLoC),  la longitud del lado de las celdas en la que se divide el LoC para discretizarlo (Paso),  el ancho inicial de la pista de cobre (anchoCu_inicial),  el ancho mínimo que se le puede dar a la pista de cobre debido a las limitaciones técnicas (Emin),  el espesor de la pista de cobre (EspesorCu),  el espesor del plano de cobre (EspesorplanoCu),  el coeficiente que se necesita aplicar sobre la resistencia lateral teórica entre dos celdas sin pista de cobre (coefCOMSOL_Resistencia_lateral_sincu),  la corriente inicial suministrada (I0), 17 17 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip  las coordenadas de la zona activa (Xmin,Ymin, Xmax, Ymax), es decir, la zona en la que se desea conseguir la temperatura objetivo, ya que si la zona en la que se desea una temperatura constante es muy amplia, sería imposible alcanzar dicha temperatura solo cambiando los anchos de la pista; y  las coordenadas en las que se encuentra el plano de cobre (XminCu,YminCu, XmaxCu, YmaxCu). Este plano de cobre, es importante colocarlo en la zona en la que se dispone la resistencia, porque cumple un papel decisivo en la homogeneización de la temperatura. classdef MisCons properties (Constant = true) Tcu=95; %degC ConConveccion=1; DimLoC=64; %en mm, tamaño del LoC Paso=0.5; %en mm, tamaño de la cuadrícula de diseño anchoCu_inicial=500; %en micras, ancho pista cobre inicial Emin=200; %en micras, ancho mínimo de pista EspesorCu=18e-6; %en m, espesor pista cobre EspesorplanoCu=100e-6; %en m, espesor plano Cu coefCOMSOL_Resistencia_lateral_sincu=1.85; I0=1.02; %en A, corriente inicial %área activa Xmin=15; Ymin=32; Xmax=25; Ymax=65; %área de la placa de cobre XminCu=3; YminCu=4; XmaxCu=46; YmaxCu=100; end end Figura 2-8. Imagen de una línea que define el recorrido de la pista de cobre, en la que se ha marcado, en rojo, el borde de la zona activa y, en azul, el de la zona ocupada por el plano de cobre. Descripción de la solución 18 18 2.2.2 Lectura de la línea que define la forma de la pista de cobre El código encargado de la lectura de la línea que define la forma de la pista de cobre es leelinea.m. Esta parte del código recibe como entrada:  Resistencia_variable.dxf, que es el fichero en el que se encuentra el diseño de la línea que define el recorrido de la pista de cobre, realizado previamente con un programa de diseño como, por ejemplo, L-Edit. Y se obtiene como salida:  linea.dat, que es el archivo que contiene los valores de los anchos de la pista de cobre, obtenido a partir del fichero de entrada, para que Matlab pueda operar con los datos que contiene. Esta conversión se realiza gracias a la función DXFtoolq.m, obtenida de modificar la librería DXFtool.m. %mis variables global primervertex dimmat paso primervertex=1; global milinea tamano_PCB tamano_PCB=MisCons.DimLoC; %LoC cuadrado de 64mm de lado, en unidades de 1 mm paso=MisCons.Paso; % 0.5 mm/cuadrado dimmat=round(tamano_PCB/paso); %64/0.5=128 cuadrados/lado milinea=255*ones(dimmat,dimmat); %un pixel es 1mm %matriz de dimensiones 128x128 con un 255 en cada elemento de la matriz % read file and plot dxf = DXFtoolq('Resistencia_variable.dxf'); 2.2.3 Cálculo y definición de los anchos de la pista de cobre El fichero resistencia_automatica.m contiene el código principal, en el que se realizan los cálculos necesarios para definer los anchos de la pista de cobre. Esta parte del código recibe como entrada:  linea.dat, que es el archivo que contiene los valores de los anchos de la pista de cobre, obtenido de leelinea.m. Y genera a la salida:  matrizTODO.dat y  parametros.dat, que son, respectivamente, el archivo que contiene los valores de los anchos de la pista de cobre en todo el LoC tras haber sido calculados para conseguir el objetivo deseado, y el archivo con los valores obtenidos de los parámetros más relevantes. Una vez cargados los datos del archivo de entrada en una matriz (Me), se calculan las resistencias que no dependen del ancho de la pista de cobre, es decir, las que mantienen un valor constante durante todo el proceso de cálculo. Para calcularlas se emplean:  la longitud del lado de la celda (L),  los espesores de cada capa (Hpcb, Hmet, Hplanocu) ,  las conductividades térmicas (Kpcb, Kmet, Kcu),  el coeficietne de película de convección (h) y 19 19 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip  los coeficientes necesarios calculados en el apartado 2.1.1 (coefCOMSOLlateral_sincu, coef_vertical_met =4.78, coef_vertical_pcb). Las resistencias finales que se usarán en el cálculo de los anchos son:  Rlsincu. Es la resistencia térmica lateral por conducción entre dos celdas sin cobre.  Rlsincu_planoCu. Es la resistencia térmica lateral por conducción entre dos celdas sin pista de cobre, pero que sí contienen al plano de cobre.  Rhl. Es la resistencia térmica lateral por convección entre las celdas que se encuentran en el contormo del LoC y el ambiente.  Ramb. Es la resistencia térmica vertical total. Aúna la resistencia vertical por conducción debida al metacrilato y la debida al PCB, y también incluye las resistencias convectivas de la superfecie superior e inferior de cada celda. %resistencias térmicas coefCOMSOLlateral_sincu=MisCons.coefCOMSOL_Resistencia_lateral_sincu; Rpcbl=L/(Kpcb*L*Hpcb); %K/W resistencia térmica lateral pcb Rml= L/(Kmet*L*Hmet); %K/W resistencia térmica lateral metacrilato RplanoCU= L/(Kcu*L*Hplanocu); %K/W resistencia termica lateral plano cobre Rlsincu_i=Rpcbl*Rml/(Rpcbl+Rml); %resistencia total lateral sin cobre Rlsincu_planoCu_i=Rpcbl*Rml*RplanoCU/(Rpcbl*Rml+Rpcbl*RplanoCU+RplanoCU*Rml); %resistencia total lateral sin pista de cobre Rlsincu=Rlsincu_i*coefCOMSOLlateral_sincu; Rlsincu_planoCu=Rlsincu_planoCu_i*coefCOMSOLlateral_sincu; h=5; %coeficiente de película convectivo Rhl=1/(h*L*(Hpcb+Hmet)); %resistencia térmica lateral convección %resistencia vertical Rmet=Hmet/(Kmet*L*L); %K/W %resistencia térmica vertical metacrilato Rpcb=Hpcb/(Kpcb*L*L); %K/W %resistencia térmica vertical PCB Ramb=zeros(N,N); %resistencia térmica vertical total for i=1:N for j=1:N if Me(i,j)>0 %si hay cobre coef_vertical_met =4.78; coef_vertical_pcb =7.29; else %si no hay cobre if (i<N && Me(i+1,j)>0) || (i>1 && Me(i-1,j)>0) || (j<N && Me(i,j+1)>0) || (j>1 && Me(i,j-1)>0) %D=1 coef_vertical_met =2.54; coef_vertical_pcb =1.42; elseif (i<(N-1) && Me(i+2,j)>0) || (i>2 && Me(i-2,j)>0) || (j<(N-1) && Me(i,j+1)>0) || (j>2 && Me(i,j-2)>0) %D=2 coef_vertical_met =1.64; coef_vertical_pcb =0.39; elseif (i<(N-2) && Me(i+3,j)>0) || (i>3 && Me(i-3,j)>0) || (j<(N-2) && Me(i,j+3)>0) || (j>3 && Me(i,j-3)>0) %D=3 coef_vertical_met =1.19; coef_vertical_pcb =0.32; else %D>3 coef_vertical_met =1; coef_vertical_pcb =1; end end Rsup=coef_vertical_met*Rmet; %resistencia térmica vertical superior Rinf=coef_vertical_pcb*Rpcb; %resistencia térmica vertical inferior if MisCons.ConConveccion==1 %conveccion Rh=1/(h*L*L); Rsup=Rsup+Rh; Rinf=Rinf+Rh; end Ramb(i,j)=Rsup*Rinf/(Rsup+Rinf); %resistencia al ambiente, paralelo vertical PCB y metacrilato end %j end %i Descripción de la solución 20 20 Tras esto, se inicia el bucle while en el que se calcula la temperatura de cada celda y se realizan las modificaciones sobre el ancho de la pista de cobre y/o sobre la corriente aplicada. De dicho bucle se sale si:  se alcanza un error considerado razonable, es decir, si el error medio de temperatura es menor de 2ºC y en ninguna celda hay un error mayor de 1ºC (ErrorTMedio < ErrorTMedioMax y ErrorLocal = 0), o si  se ha repetido este bucle 30 veces (Pasos = PasosMax). Dentro de este bucle, lo primero que se hace es calcular las resistencias que dependen del ancho de la pista de cobre y que, por tanto, varían en cada iteración del código. Para calcularlas se emplean:  la longitud del lado de la celda (L),  los espesores de cada capa (Hpcb, Hmet,Hcu) ,  las conductividades térmicas (Kpcb, Kmet, Kcu),  el ancho de la pista de cobre (Wcu),  la resistividad del cobre (ro_Cu) y  los coeficientes necesarios calculados en el apartado 2.1.1 (coefCOMSOLlateral_concu_sincu). De aquí se obtienen:  ReCu. Es la resistencia eléctrica en cada celda debida a la circulación de corriente.  RCu. Es la resistencia térmica lateral por conducción de la pista de cobre.  Rlconcu_sincu. Es la resistencia térmica lateral por conducción entre una celda que contiene pista de cobre y otra que no. Puesto que el plano de cobre se encuentre en la zona ocupada por la pista de cobre, se entiende que cualquier intercambio entre una celda con pista de cobre y otra sin ella va a contener parte del plano de cobre.  Rlconcu. Es la resistencia térmica lateral entre dos celdas con pista de cobre y que, por tanto, contienen al plano de cobre. . . . ReCu=zeros(N,N); %resistencia eléctrica Rcu=zeros(N,N); %resistencia térmica de cobre Rlconcu_sincu=zeros(N,N); %resistencia térmica entre celda con cobre y sin Rconcu=zeros(N,N); %resistencia térmica entre celdas con cobre %se itera while(Pasos < PasosMax) && (ErrorTMedio > ErrorTMedioMax || ErrorLocal==1 ) Resist_Total=0; %resistencia eléctrica total Wcu=Me*1e-6; %se pasa a metros %matriz con los espesores del cobre en cada cuadrado (i,j) %matriz de resistencias electricas y térmicas for i=1:N for j=1:N if Me(i,j)>0 %si hay cobre %resistencia eléctrica ReCu(i,j)=ro_Cu*L/(Hcu*Wcu(i,j)); %ohm Resist_Total=Resist_Total+ReCu(i,j); %resistencia térmica lateral pista cobre Rcu(i,j)=L/(Kcu*Wcu(i,j)*Hcu); %K/W %resistencia térmica entre celda con Cu y celda sin Cu coefCOMSOLlateral_concu_sincu(i,j)=4.3089*((Wcu(i,j)/L)^(-0.214)); 21 21 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip R1(i,j)=(Wcu(i,j)/2)/(L*(Kcu*Hcu+Kpcb*Hpcb+Kmet*(Hmet-Hcu))); R2(i,j)=(L-AnchoCu(i,j)/2)/(L*(Kpcb*Hpcb+Kmet*Hmet)); Rlconcu_sincu(i,j)=R1(i,j)+R2(i,j); Rlconcu_sincu(i,j)=Rlconcu_sincu(i,j)*RplanoCU/(Rlconcu_sincu(i,j)+RplanoCU); %se añade plano de cobre Rlconcu_sincu(i,j)=Rlconcu_sincu(i,j)*coefCOMSOLlateral_concu_sincu(i,j); %resistencia térmica entre dos celdas con cobre Rlconcu(i,j)=Rlsincu_planoCu_i*Rcu(i,j)/(Rlsincu_planoCu_i+Rcu(i,j)); else %si no hay cobre ReCu(i,j)=0; Rcu(i,j)=1e6; Rlconcu_sincu(i,j)=1e6; Rlconcu(i,j)=1e6; end end %j end %i . . . end %while Tras el cálculo de las resistencias, y a partir de:  la corriente suministrada (I),  la temperatura ambiente (Tamb) y  las resistencias calculadas previamente (ReCu, Ramb, Rhl, Rlsincu, Rlsincu_planoCu, Rlconcu_sincu,Rlconcu). se generan  la matriz de conductancias térmicas (G) y  la matriz de fuentes (o sumideros) de energía (Q). Dichas matrices, a partir de la ecuación 𝑄 =𝐺𝑇, permiten calcular  las temperaturas en cada celda (T_cu), para, en base a ellas, realizar modificaciones sobre el ancho de la pista de cobre o sobre el valor de la corriente. En primer lugar, se calcula la generación térmica, así como la contribución al ambiente de cada celda. Para esto ultimo, se tiene en cuenta que 𝑄𝑘=𝑇𝑘−𝑇𝑎𝑚𝑏 𝑅𝑎𝑚𝑏 =𝑇𝑘 𝑅𝑎𝑚𝑏 −𝑇𝑎𝑚𝑏 𝑅𝑎𝑚𝑏 . Por lo que 1 𝑅𝑎𝑚𝑏 pertenece a la matriz G y 𝑇𝑎𝑚𝑏 𝑅𝑎𝑚𝑏 al ser un valor constante, que no depende de Tk, pasa sumando a la matriz Q. for i=1:N for j=1:N %para cada nodo i,j %generación térmica Q((i-1)*N+j)=Q((i-1)*N+j)+I*I*ReCu(i,j); %contribucion al ambiente Q((i-1)*N+j)=Q((i-1)*N+j)+Tamb/Ramb(i,j); %transmisión con el ambiente G((i-1)*N+j,(i-1)*N+j)=G((i-1)*N+j,(i-1)*N+j)+1/Ramb(i,j); . . . Descripción de la solución 22 22 Después, se calcula la transmisión lateral de una celda que no contiene pista de cobre con las cuatro celdas (también llamadas nodos) vecinas. Al igual que para el caso anterior, hay que tener en cuenta que la transmisión lateral con el ambiente es 𝑄𝑘=𝑇𝑘−𝑇𝑎𝑚𝑏 (𝑅𝑙𝑠𝑖𝑛𝑐𝑢+𝑅ℎ𝑙)=𝑇𝑘 (𝑅𝑙𝑠𝑖𝑛𝑐𝑢+𝑅ℎ𝑙)−𝑇𝑎𝑚𝑏 (𝑅𝑙𝑠𝑖𝑛𝑐𝑢+𝑅ℎ𝑙) , por lo que se opera de la misma manera. Y para el intercambio entre dos celdas contiguas se tiene que 𝑄𝑘=𝑇𝑘−𝑇𝑘+𝑁 𝑅=𝑇𝑘 𝑅− 𝑇𝑘+𝑁 𝑅, donde R puede ser Rlsincu_planoCu, Rlconcu_sincu o Rlsincu. En este caso 1 𝑅 se sumará a Gk y −1 𝑅 se sumará a Gk+N. A continuación se muestra el código para la transmisión con una de las celdas vecinas, para las celdas vecinas restantes se usa el mismo código cambiando solo los índices según corresponda. if Me(i,j)==0 %nodo sin cobre %nodo i+1,j %transmisión de calor con el nodo de abajo(i+1) if i==N %si es nodo de un canto (Rhl), no hay cobre al lado (Rlsincu) G((i-1)*N+j,(i-1)*N+j)=G((i-1)*N+j,(i-1)*N+j)+1/(Rlsincu+Rhl); %contribucion al ambiente Q((i-1)*N+j)=Q((i-1)*N+j)+Tamb/(Rlsincu+Rhl); else %si es nodo interior if i>=MisCons.YminCu && i<=MisCons.YmaxCu && j>=MisCons.XminCu && j<=MisCons.XminCu %si hay plano de cobre if Me(i+1,j)==0 %si no hay cobre en el adyacente (Rlsincu_planoCu) G((i-1)*N+j,(i-1)*N+j) =G((i-1)*N+j,(i-1)*N+j) +1/Rlsincu_planoCu; G((i-1)*N+j,(i-1+1)*N+j)=G((i-1)*N+j,(i-1+1)*N+j)-1/Rlsincu_planoCu; else %si hay cobre en el adyascente (Rlconcu_sincu) G((i-1)*N+j,(i-1)*N+j) =G((i-1)*N+j,(i-1)*N+j) +1/Rlconcu_sincu(i+1,j); G((i-1)*N+j,(i-1+1)*N+j)=G((i-1)*N+j,(i-1+1)*N+j)-1/Rlconcu_sincu(i+1,j); end else %si no hay plano de cobre if Me(i+1,j)==0 %si no hay cobre en el adyacente (Rlsincu) G((i-1)*N+j,(i-1)*N+j) =G((i-1)*N+j,(i-1)*N+j) +1/Rlsincu; G((i-1)*N+j,(i-1+1)*N+j)=G((i-1)*N+j,(i-1+1)*N+j)-1/Rlsincu; else %si hay cobre en el adyascente (Rlconcu_sincu) G((i-1)*N+j,(i-1)*N+j) =G((i-1)*N+j,(i-1)*N+j) +1/Rlconcu_sincu(i+1,j); G((i-1)*N+j,(i-1+1)*N+j)=G((i-1)*N+j,(i-1+1)*N+j)-1/Rlconcu_sincu(i+1,j); end end end . . . Por ultimo, se calcula la transmisión lateral de una celda con pista de cobre con sus cuatro celdas vecinas. Tras esto, se obtienen los valores de las temperaturas en cada celda. Se opera igual que en el caso anterior, pero variando las resistencias empleadas y sin comprobar si hay plano de cobre, ya que, al tratarse de una celda con pista de cobre, siempre va a contener parte del plano de cobre. Se muestra el código para la transmisión con una de las celdas vecinas, para las celdas vecinas restantes se usa el mismo código cambiando solo los índices según corresponda. else %Me(i,j)>0 nodo con cobre %nodo i+1,j %transmisión de calor con el nodo de abajo(i+1) if i==N %si es nodo de un canto (Rhl), no hay cobre al lado (Rlsincu) G((i-1)*N+j,(i-1)*N+j)=G((i-1)*N+j,(i-1)*N+j)+1/(Rlsincu+Rhl); %contribucion al ambiente Q((i-1)*N+j)=Q((i-1)*N+j)+Tamb/(Rlsincu+Rhl); else %si es nodo interior if Me(i+1,j)==0 %si no hay cobre en el adyacente (Rlconcu_sincu) G((i-1)*N+j,(i-1)*N+j) =G((i-1)*N+j,(i-1)*N+j)+1/Rlconcu_sincu(i,j); G((i-1)*N+j,(i-1+1)*N+j)=G((i-1)*N+j,(i-1+1)*N+j)-1/Rlconcu_sincu(i,j); else %si hay cobre en el adyacente (Rlconcu) G((i-1)*N+j,(i-1)*N+j) =G((i-1)*N+j,(i-1)*N+j)+1/Rlconcu(i,j); G((i-1)*N+j,(i-1+1)*N+j)=G((i-1)*N+j,(i-1+1)*N+j)-1/Rlconcu(i,j); end end 23 23 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip . . . end %Me end %j end %i %fin generacion ecuaciones %cálculo temperaturas por matriz inversa T=G\Q; T_cu = reshape(T(1:N*N),[N,N]); %para volver a pasar del vector de T(k) a la matriz de posiciones T(i,j) T_cu=T_cu'; Pasos=Pasos+1; Tras el cálculo de temperaturas, se procede a ver el error, para actuar modificando el ancho de la pista (Me) y/o la corriente aplicada (I). Se compara la diferencia entre:  la temperatuar objetivo (Tobjetivo) y  la temperatura calculada en cada celda (T_cu), con  el máximo error local de temperatura permitido (MaxErrorLocalT). Si Tobjetivo-T_cu(i,j)>MaxErrorLocalT,  se indica que hay un error local (ErrorLocal=1) y  se reduce el ancho de la pista de cobre (Me(i,j)= Me(i,j)-PasoAncho).  Pero si al reducir el ancho, este es menor que el mínimo permitido, se sube la corriente (I=I+deltaI). Si, por el contrario, T_cu(i,j)-Tobjetivo>MaxErrorLocalT,  se indica que hay un error local (ErrorLocal=1) y  se aumenta el ancho de la pista de cobre (Me(i,j)= Me(i,j)+PasoAncho).  Pero si al aumentar el ancho, este es mayor que el máximo permitido, se disminuye la corriente (I=I-deltaI). for i=MisCons.Ymin:MisCons.Ymax for j=MisCons.Xmin:MisCons.Xmax %en zona activa if Me(i,j)>0 %nodo con cobre ErrorTMedio=ErrorTMedio+(Tobjetivo-T_cu(i,j))^2; NumNodos=NumNodos+1; %contabilizamos el nodo porque tiene pista cobre %ajuste resistencia if Tobjetivo-T_cu(i,j)>MaxErrorLocalT ErrorLocal=1; Me(i,j)= Me(i,j)-PasoAncho; %en um if Me(i,j)<emin Me(i,j)=emin; subir_corriente=1; %se aumenta la corriente porque no se puede disminuir la pista end elseif T_cu(i,j)-Tobjetivo>MaxErrorLocalT ErrorLocal=1; Me(i,j)= Me(i,j)+PasoAncho; if Me(i,j)>emax Me(i,j)=emax; bajar_corriente=1; %se disminuye la corriente porque no se puede aumentar la pista end end end %Me end %j Descripción de la solución 30 30 Paso 8. Si converge y se finaliza la simulación con éxito ir al paso 8, si no, ir al paso 10. Paso 9. Si el tamaño de la zona activa definida es el deseado ir al paso 12, si no, ir al paso 9. Paso 10. Se incrementa el tamaño de la zona activa, se simula de nuevo el código y se vuelve al paso 3. Paso 11. Si la mayor zona activa para la que se ha conseguido que converja anteriormente resulta suficiente para el objetivo buscado, ir al paso 12, si no, ir al paso 11. Paso 12. Se vuelve a L-Edit (o equivalente) para modificar el recorrido inicial de la línea separando las líneas que la conforman. Inicialmente se separarán al menos dos de las líneas centrales de la resistencia al doble de la separación inicial. Según en qué situación se llegue a este punto será necesario separar más líneas o separarlas más espacio. En cualquier caso se operará de manera gradual. Una vez realizada una modificación en el recorrido de la línea ir al paso 2. Paso 13. Una vez conseguida una forma de pista de cobre válida, se introduce el archivo dado por el programa de cálculo en L-Edit (o equivalente) y se deshace la discretización de la línea (en L-Edit hacer “merge”). El resultado se introduce en el programa de simulación (COMSOL) en el que previamente se ha diseñado el escenario sobre el que se quiere trabajar y se va al paso 13. Paso 14. Si en COMSOL se consigue el resultado esperado ir al paso 14, si no ir, al paso 11. Paso 15. Una vez obtenida una simulación en COMSOL satisfactoria, se finaliza el proceso y se procede a implementar la pista obtenida y a comprobar su validez en laboratorio. 31 3 CASO PRÁCTICO e ha llevado a cabo un ejemplo en el que se desea obtener una temperatura de 65ºC en tres puntos concretos de un LoC de 68mm de lado, compuesto por una capa inferior de PCB de 1.6mm de espesor y una capa superior de metacrilato de 5mm de espesor. La línea de cobre que conforma la resistencia se colocorá entre ambas capas. Los puntos en los que se quiere obtener la temperatura objetivo son tres círculos, colocados en las posiciones que se muestran en la Figura 3-1, que tienen una composición distinta al resto del LoC, no tienen una capa de PCB y otra de metacrilato, sino que están compuestos en su totalidad de metacrilato. Además, por encima de esos tres puntos no puede pasar la pista de cobre, ya que se necesita que estén libres. Figura 3-1. Planta obtenida de COMSOL del LoC del ejemplo en el que se desea colocar una pista de cobre. Los tres círculos que se muestran son las zonas en la que se desea obtener una temperatura de 65ºC y por encima de los cuales no puede pasar la pista de cobre. Inicialmente se dibujó como recorrido de la pista de cobre una línea de separación uniforme de 1mm que bordeaba los tres puntos de interés, tal y como se muestra en la Figura 3-2. Tras esto se introdujo la resistencia obtenida en Matlab, se definió una zona activa que abarcaba una gran parte de la pista de cobre, se añadió un plano de cobre en la zona ocupada por la resistencia y se simuló. Este resultado no convergió, pero se decidió introducir la pista de cobre resultante, mostrada en la Figura 3-3, en COMSOL, para comprobar su comportamiento y así tener una idea de cómo sería el mapa de temperaturas con una pista de cobre con esta separación entre líneas. S Caso práctico 32 32 Figura 3-2. Primer diseño de la línea que define el recorrido de la pista de cobre con el contorno de la zona activa marcado en rojo y el de la zona ocupada por el plano de cobre marcado en azul. Las tres zonas en las que no hay línea de cobre se corresponden con las zonas ocupadas por los tres círculos de la Figura 3-1. Figura 3-3. Pista de cobre obtenida en L-Edit tras deshacer la discretización en el resultado obtenido en Matlab para una línea como la mostrada en la Figura 3-2. 33 33 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip Figura 3-4. Corte de temperatura en el plano de la pista de cobre, obtenido en COMSOL, para una pista de cobre como la mostrada en la Figura 3-3. El resultado obtenido en COMSOL se muestra en la Figura 3-4, donde se ve que la temperatura alcanzada en los puntos de interés es algunos grados menor que la que se obtiene en la zona por la que pasa la resistencia. Además, el gradiente de temperatura en los círculos era alto, variando 4ºC desde el centro al perímetro, y diferente en los tres puntos, siendo el de la izquierda el que se encontraba a menor temperatura y el central el de mayor temperatura. Por ello se decidió separar las líneas centrales de la resistencia. Se hicieron varias pruebas, separando las líneas en zonas progresivamente más amplias e incluyendo pista de cobre por los huecos que se habían dejado inicialmente junto a los círculos (Figura 3-5). Todas estas formas se simularon en COMSOL con la misma zona activa que la que aparece en la Figura 3-2 y, a pesar de que ninguna convergió, se decidió introducir el resultado en COMSOL y ver cómo variaba el mapa de temperaturas en el plano de la pista de cobre en función del recorrido de esta (Figura 3-6). R1 R2 Caso práctico 34 34 R3 R4 R5 R6 Figura 3-5. Distintos diseños de la línea que define el recorrido de la pista de cobre, numerados según el orden cronológico en el que fueron creados. R1 R2 R3 R4 35 35 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip R5 R6 Figura 3-6. Mapa de temperatura del corte realizado en el plano de la pista de cobre, obtenido con las distintas líneas que definen la pista de cobre de la Figura 3-5. El primer diseño que se realizó fue el de R1, separando algunas de las pistas centrales 2mm entre sí, y se vio una clara mejora con respecto al resultado obtenido en la Figura 3-4. Los tres círculos tenían unas temperaturas similares y además, el gradiente de temperatura dentro de ellos era ligeramente menor. Así que se decidió mantener esta forma y añadir pista de cobre por debajo de uno de los círculos, para comprobar qué efecto tendría ese cambio en el recorrido. Por lo tanto, se rodeó uno de los círculos por completo con pista de cobre y, por simplificar el trabajo, se bordeó solo uno para ver el resultado y decidir si era necesario bordear los demás. Primero se añadió una parte de la pista de cobre con una distancia de 1mm entre las líneas, obteniendo R2, y después con una distancia de 2mm, obteniendo R3, y, aunque el resultado de esta segunda era algo mejor, realmente ninguno de los dos era válido. Ya que, aunque, como se observa en la Figura 3-6, el hecho de incluir parte de la pista de cobre por debajo de uno de los círculos hace que la diferencia de temperaturas en toda la zona que ocupa la resistencia sea menor, al estudiar el círculo que ha sido rodeado por la pista de cobre, se observa que la distribución de temperatuars dentro de él es menos uniforme que en los otros círculos y que el gradiente de temperatura es mayor. Por lo que estas dos propuestas fueron rechazadas y no se probó cuál sería el resultado si se bordearan todos los círculos con la pista de cobre. Desechada la idea de bordear todos los círculos, se decidió seguir separando líneas de la pista de cobre, obteniendo así R4, R5 y R6 sucesivamente. En todas ellas el resultado era parecido, siendo la temperatura en todos los círculos bastante similar. Sin embargo, en R5 la diferencia de temperatura entre el centro del círculo y su perímetro era ligeramente inferior (2ºC), por lo que se decidió tomar esta forma, aunque cualquiera de las otras dos podría también haber sido válida. Tras ver los distintos resultados y comprobar que la temperatura en los círculos siempre es menor que en la pista de cobre, se llegó a la conclusión de que no se podía coger una zona activa que incluyese una zona con pista de cobre y otra sin ella, ya que, al ser la resistencia del cobre más pequeña que la del resto de los materiales, este se convierte en dominante en la transmisión de calor y, por tanto, las zonas que contienen a la pista de cobre siempre tendrán una temperatura mayor a la de los círculos en los que nos interesa obtener la temperatura objetivo. Por ello, se decidió tomar una zona activa que no incluyese a los círculos (Figura 3-7) y aumentar la temperatura objetivo en el código de Matlab a 68ºC, para así obtener una temperatura algo menor en los círculos. La zona activa seleccionada es tal que no incluya a los círculos, que esté próxima a todos ellos, y, por tanto, centrada, y que permita que la simulación converja. Una vez obtenida dicha zona activa e introducida la línea en el programa de cálculo, este devuelve el resultado de la Figura 3-8, y dicho resultado se dio para una corriente de 0.855A. Tras ser introducido en COMSOL ,se obtiene el corte de temperatura de la Figura 3-9. Caso práctico 36 36 Figura 3-7. Diseño de la línea que define el recorrido de la pista de cobre correspondiente a R5 en la Figura 3-5 y 3-6, con el contorno de la zona activa marcado en rojo y el de la zona ocupada por el plano de cobre marcado en azul. Zona activa definida tras comprobar que no se puede alcanzar la misma temperatura en los círculos que en la pista de cobre. Figura 3-8. Pista de cobre obtenida en L-Edit tras deshacer la discretización en el resultado obtenido en Matlab para una línea correspondiente a R5 en la Figura 3-5 y 3-6. 37 37 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip Figura 3-9. Corte de temperatura en el plano de la pista de cobre, obtenido en COMSOL, para la pista de cobre de la Figura 3-8. En este último resultado de la Figura 3-9, se aprecia cómo el salto de temperatura es menor, no solo en los puntos en los que se desea obtener la temperatura objetivo, sino también en comparación con la temperatura en la resistencia. Además, la distribución de temperatura en los tres puntos es bastante similar. También es importante destacar, que el ancho de la resistencia ha sufrido pocas modificaciones, por lo que queda de manifiesto que es más relevante la forma dada inicialmente a la pista de cobre, que la modificación del ancho realizada por Matlab. En el resultado de COMSOL se obutvo una tensión de 1.63V para una corriente de 0.855A, por lo que la resistencia de la pista de cobre es 𝑅 =𝑉 𝐼=1.63 0.855 =1.91Ω. Por último, se fabricó esta resistencia en una de las caras de un PCB de doble cara (Figura 3-10) y se hicieron varias comprobaciones en el laboratorio. En primer lugar se colocó un bloque de metacrilato sobre el PCB para semejarlo lo máximo posible al modelo sobre el que se habían realizado las simulaciones previas. Y se colocó un sensor de temperatura entre el PCB y el metacrilato para poder medir la temperatura en distintos puntos de dicho plano. Inicialmente dicho sensor se colocó sobre la pista de cobre en el centro de la zona inferior (T1), y se sometió la pista a una corriente, que se fue aumentando hasta obtener en el sensor una temperatura similar a 65ºC. La corriente alcanzada fue de 1.57A y la tensión obtenida para dicha corriente fue de 3V, por lo que la resistencia eléctrica de la pista de cobre es 𝑅 =𝑉 𝐼=3 1.57 =1.91Ω, coincidiendo con el resultado obtenido en COMSOL. Se usó un termómetro que muestra un mapa de temperaturas, para ver la temperatura en la superficie del metacrilato, de lo cual se muestra una imagen en la Figura 3-11, y se colocó el sensor de temperatura ubicado entre el PCB y el metacrilato en distintos puntos de ese plano para medir sus temperaturas. Se obtuvieron los siguientes resultados, para cada uno de los puntos en los que se midió la temperatura (Figura 3-10): Caso práctico 38 38  T1 = 65.3ºC. Punto inferior central sobre la pista de cobre.  T2= 66.2 ºC. Punto perteneciente a la pista de cobre, en la zona inferior izquierda.  T3 = 57.8ºC. Punto correspondiente al hueco inferior izquierdo.  T4 = 59.6 ºC. Punto correspondiente al hueco superior.  T5 = 59.6ºC. Punto correspondiente al hueco inferior derecho. T4 T3 T2 T1 T5 Figura 3-10. Pista de cobre de la Figura 3-8 fabricada en un PCB de doble cara. Los puntos azules se corresponden con los distintos puntos en los que se ha obtenido una medida de la temperatura. La primera temperatura se ha tomado para definir la corriente que se le suministraba y la segunda, para poder comparar la temperatura en dos puntos pertenecientes a la pista de cobre. Se comprueba que se obtiene una diferencia de temperatura de aproximadamente 1ºC, lo cual es aceptable. En las temperaturas tomadas en los huecos correspondientes a las posiciones que ocuparían los círculos, se observa que el punto central y el de la derecha se encuentran a la misma temperatura (T4 = T5), sin embargo, no hay simetría y el punto de la izquierda (T3) se encuentra a una temperatura que es inferior en 2ºC. Además, estas tres zonas se encuentran a unos 6ºC menos que las zonas que se encuentran sobre la pista de cobre. Sin ser el obtenido un resultado plenamente satisfactorio, puede ser aceptado, ya que la diferencia de temperatura entre el hueco de la izquierda y los otros no es excesivamente grande y las condiciones en las que se están realizando estas medidas no son exactamente las mismas para las que se ha simulado. Entre las diferencias de las condiciones entre la experimentación física y la simulación se encuentra que en el primer caso no se han realizado los agujeros, luego la zona en la que se está midiendo no está formada solo por metacrilato, sino que también hay PCB y cobre. También en la experimentación, el bloque de metacrilato se ha movido para poder colocar el sensor de temperatura en distintas posiciones, dejando en ocasiones parte del PCB en contacto con el aire. Por todo ello, es lógico, que los resultados obtenidos en la experimentación física no coincidan por completo con los obtenidos en la simulación en COMSOL. 39 39 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip Figura 3-11. Mapa de temperaturas de la superficie del metacrilato. La zona roja se corresponde con la de mayor temperatura, mientras que la azul se corresponde con la de menor temperatura. La temperatura media de la superficie es de 54.9ºC. Figura 3-12. Puesto de trabajo en el laboratorio. En la imagen se distingue la fuente de energía al fondo, marcando la corriente (rojo) y tensión (verde) a la que se somete la pista de cobre. Esta está conectada con los cables rojo y negro a la placa de PCB, la cual tiene un bloque de metacrilato sobre ella sujeto por dos pinzas. Se observa también, un cable verde y blanco que sale de la placa por su parte superior, en cuyo extremo se encuentra el sensor de temperatura, y que se conecta al aparato que aparece en primer plano, por el que se muestra el valor obtenido con el sensor. Por último, delante de la placa, se encuentra el termómetro, mostrando por pantalla el mapa de temperaturas de la superficie del metacrilato. Anexo A: Manual de usuario 46 46  Durante la simulación, muestra o la evolución de (de izquierda a derecha y de arriba a abajo):  el error medio de temperatura en ºC,  el paso de modificación del ancho de la pista de cobre en µm,  el ancho medio de la pista de cobre (cruz azul), su ancho máximo (círculo turquesa) y su ancho mínimo (círculo rojo) en µm y dentro de la zona activa,  el mapa de temperaturas en ºC, así como los puntos de máxima (cruz roja) y mínima (círculo azul) temperatura dentro de la zona activa,  la corriente aplicada en A y  la temperatura máxima (cruz roja) y mínima (círculo azul) alcanzadas en la zona activa; y o el perfil de temperaturas de la línea media del LoC en cada iteración. Figura AA-0-4. Gráficas mostradas en Matlab para el seguimiento de parámetros. De izquierda a derecha y de arriba abajo se muestran: En la gráfica 1, el error medio de temperatura en ºC. En la gráfica 2, el paso de modificación del ancho de la pista de cobre en µm. En la gráfica 3, el ancho medio de la pista de cobre (cruz azul), su ancho máximo (círculo turquesa) y su ancho mínimo (círculo rojo) en µm y dentro de la zona activa. En la figura 4, el mapa de temperaturas en ºC, así como los puntos de máxima (cruz roja) y mínima (círculo azul) temperatura dentro de la zona activa. En la gráfica 5, la corriente aplicada en A. En la figura 6, la temperatura máxima (cruz roja) y mínima (círculo azul) alcanzadas en la zona activa. 47 47 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip Figura AA-0-5. Perfil de temperaturas de la línea media del LoC en una iteración.  Al finalizar la simulación, muestra dos imágen con los espesores de la línea de cobre discretizados. Figura AA-0-6. Pista de cobre obtenido tras el cálculo de los anchos. A la izquierda se muestran los anchos de la pista mediante escala de grises, siendo el negro un ancho 0. A la derecha se muestran los anchos de manera discretizada, de manera que se ve la pista de cobre formada por pequeños cuadrados de distinto tamaño. Anexo A: Manual de usuario 48 48 5. Importar la máscara en L-Edit desde mascara_cu.dxf. Al importarla, indicar 0.001 en la escala (1dxf = 1µm). Una vez importada, hacer visible la capa “21”, en la que se encuentra la máscara, seleccionar toda la pista de cobre y hacer “merge”. El resultado se exporta en DXF en el archivo MiResistencia.dxf. Figura AA-0-7. Pista de cobre obtenida en Matlab importada en L-Edit. Se encuentra seleccionada la línea completa y marcado el icono de “merge”. 6. Importar la máscara MiResistencia.dxf generada en L-Edit en el fichero resistencia2.0.mph de COMSOL para poderla simular. Este fichero ha sido previamente generado con las condiciones del LoC que se desea simular. Para cargar la máscara, ir al apartado Componentes → Geometría → Plano de trabajo cobre → Geometría de plano → Importar MiResistencia.dxf → Seleccionar archivo en el navegador → Construir todo. Introducir en el apartado Parámetros → Ipcr95 el valor de la corriente obtenida en Matlab y simular en COMSOL. 49 49 Sistema de cálculo automático de resistencias PCB para calentadores en Lab-on-chip Figura AA-0-8. Mapa de temperaturas en el plano de la pista de cobre obtenido en COMSOL. A la izquierda con el rango completo de temperaturas del LoC. A la derecha con el rango de temperaturas permitido, es decir Tobj±1ºC. 51 REFERENCIAS [1] P. Streit, J. Nestler, R. Schulze, A.Shaporin and T. Otto, «Investigation on the Temperature Distribution of Integrated Heater Configuration in a Lab-on-a-Chip System,» 2017 18th International Conference on Thermal, Mechanical and Multi-Physics Simulation and Experiments in Microelectronics and Microsystems (EuroSimE), pp. 1-8, 2017, doi: 10.1109/EuroSimE.2017.7926229 [2] J. Martinez-Quijada et al., «Deterministic Design of Thin-Film Heaters for Precise Spatial Temperature Control in Lab-on-Chip Systems, ,» Journal of Microelectromechanical Systems, vol. 25, no. 3, pp. 508516, 2016, doi: 10.1109/JMEMS.2016.2536561.