Repositorio Institucional de Documentos

Trabajo Fin de M´aster M´aster en Mec´anica Aplicada Programa Oficial de Posgrado en Mec´anica Computacional Curso 2010-2011 Modelos de comportamiento para suspensiones de nanotubos de carbono. Aplicaci´on a la simulaci´on de procesos de fabricaci´on. Rosa Mar´ıa Monge Prieto Septiembre de 2011 Director: El´ıas Cueto Prendes Co-Director: David Gonz´alez Ib´a˜nez Departamento de Ingenier´ıa Mec´anica Escuela de Ingenier´ıa y Arquitectura Universidad de Zaragoza Resumen. Los nanotubos de carbono son estructuras de escala nanom´etrica formadas a partir de l´aminas de grafito enrolladas sobre s´ı mismas. Dependiendo del n´umero de capas, los nanotubos (CNTS por sus siglas en ingl´es) pueden ser single-walled (una s´ola capa) o multiwalled (multicapa). Tanto unos como otros presentan unas propiedades mec´anicas, t´ermicas y el´ectricas que los hacen muy atractivos para el desarrollo de materiales compuestos. Para el buen aprovechamiento de estas caracter´ısticas es necesario conocer la orientaci´on de los CNTS en el seno de la matriz polim´erica, ya que el comportamiento mec´anico depende de esta orientaci´on y el comportamiento el´ectrico y t´ermico es notablemente mejor cuanto menos aglomeraciones presenten los nanotubos. Para este objetivo es de suma importancia disponer de t´ecnicas num´ericas de simulaci´on que permitan una evaluaci´on r´apida de condiciones de fabricaci´on de estos compuestos, as´ı como para la validaci´on de hip´otesis realizadas a partir de los ensayos experimentales realizados en el laboratorio. El proceso de fabricaci´on que se ha simulado en el trabajo ha sido el de spin-coating, presentando ´este unas particularidades (grades deformaciones, superficie libre, etc.) que hacen que la simulaci´on num´erica mediante los m´etodos tradicionales con malla no sea f´acil, siendo m´as adecuado un m´etodo sin malla, m´as concretamente, el M´etodo de los Elementos Naturales (MEN), ya que este m´etodo nos permite la descripci´on Lagrangiana actualizada de la cinem´atica del flujo. Otro de los aspectos que merece la pena ser estudiado en estos compuestos es su particular reolog´ıa. Para ellos hay que recurrir a modelos microestructurales si se quiere, mediante Din´amica Browniana, simular lo que ocurre a escala macrosc´opica. Los nanotubos de carbono (tanto los single-walled como los multi-walled) pueden ser tratados qu´ımicamente, variando de esta manera su comportamiento reol´ogico, por lo que el mismo modelo no sirve para explicar ambos comportamientos (el de los que est´an tratados y el de los que no). Se tiene as´ı el modelo de Orientaci´on (nanotubos tratados qu´ımicamente) y el de Agregaci´on/Orientaci´on (nanotubos no tratados). Los resultados de las simulaciones correspondientes a ambos modelos son presentados en este trabajo. ´ Indice general I Memoria 1 1. Introducci´on 3 1.1. Reolog´ıa de las suspensiones de nanotubos de carbono . . . . . . . . . . . . 4 1.1.1. La formaci´on de bandas helicoidales en suspensiones de CNTS no tratados qu´ımicamente. . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.2. Reolog´ıa experimental de las suspensiones de CNTS. . . . . . . . . 4 2. Un modelo de tipo FENE 7 2.1. Introducci´on. .................................. 8 2.2. Formulaci´on del modelo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.1. Din´amica del pol´ımero. . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.2. Tensor de tensiones. . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3. Funciones del material para un flujo a cortadura. . . . . . . . . . . . . . . 10 2.4. Elm´etodonum´erico. .............................. 11 2.4.1. Resultados de las simulaciones. . . . . . . . . . . . . . . . . . . . . 12 3. Proceso de spin-coating 15 4. M´etodo de los Elementos Naturales 17 4.1. Interpolaci´on por vecinos naturales. . . . . . . . . . . . . . . . . . . . . . . 17 5. Dos modelos basados en la teor´ıa cin´etica y la orientaci´on de las fibras 21 5.1. Modelo de Fokker-Plack . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 I 5.1.1. El modelo num´erico. . . . . . . . . . . . . . . . . . . . . . . . . . . 21 5.1.2. Relaci´on entre modelo y reolog´ıa . . . . . . . . . . . . . . . . . . . 24 5.1.3. Discretizaci´on del problema. . . . . . . . . . . . . . . . . . . . . . . 25 5.1.4. Resultados de las simulaciones. . . . . . . . . . . . . . . . . . . . . 27 5.2. El modelo de Agregaci´on/Orientaci´on. . . . . . . . . . . . . . . . . . . . . 29 5.2.1. El modelo num´erico. . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.2.2. Resultados de las simulaciones. . . . . . . . . . . . . . . . . . . . . 30 6. Conclusiones 35 6.1. Trabajofuturo. ................................. 35 II Anexos 39 A. Tecnolog´ıa y compuestos de nanotubos de carbono. 41 A.1. Subestructura y morfolog´ıa de los nanotubos de carbono. . . . . . . . . . . 42 A.1.1. Estructura de los nanotubos de carbono. . . . . . . . . . . . . . . . 42 A.1.2. Morfolog´ıa de los nanotubos de carbono. . . . . . . . . . . . . . . . 43 A.2. Procesado de los nanotubos de carbono para materiales compuestos. . . . . 43 A.2.1.Ablaci´onl´aser. ............................. 44 A.2.2.Descargadearco............................. 44 A.2.3.CVD................................... 45 A.3. Caracterizaci´on de los nanotubos de carbono. . . . . . . . . . . . . . . . . 46 A.4. Composites de nanotubos de carbono . . . . . . . . . . . . . . . . . . . . . 48 A.4.1. Procesado y caracterizaci´on de compuestos polim´ericos con nanotubosdecarbono.............................. 48 A.4.2. Compuestos cer´amicos y met´alicos. . . . . . . . . . . . . . . . . . . 48 B. Resultados 49 II C. Spin-coating 53 C.1. Descripci´on del proceso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 C.2. Velocidad de rotaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 C.3.Aceleraci´on ................................... 55 C.4.Secado...................................... 55 C.5.Gr´aficosdelproceso............................... 55 C.6. Problemas del proceso de spin-coating .................... 56 D. Formas-α. 61 D.1.Introducci´on. .................................. 61 D.2. Teor´ıa de las formas-α.............................. 62 D.2.1. C´omo escoger el valor de α. ...................... 63 D.2.2. Problemas con la t´ecnica de las formas-α. .............. 64 E. Ampliaci´on del modelo de Agregaci´on/Orientaci´on. 67 E.0.3. El modelo num´erico. . . . . . . . . . . . . . . . . . . . . . . . . . . 67 E.0.4. Resultados de las simulaciones. . . . . . . . . . . . . . . . . . . . . 68 III IV ´ Indice de figuras 1.1. Nanotubos de carbono single-walled (SWCNTS) y multi-walled (MWCNTS). ...... 3 1.2. Evoluci´on temporal de la microestructura ´optica de una disoluci´on de 0.03 %de MWCNT0s para distintas distancias de placas y un ritmo de cizalladura constante de 0.5 s−1.. . . 5 1.3. Viscosidad aparente (ηa) para (a) epoxi, (b) disoluci´on con 0.25 %de CNTS tratados y (c) disoluci´on con 0.25 %de CNTS no tratados qu´ımicamente. ............ 6 1.4. Viscosidad aparente (ηa) en funci´on de la velocidad de deformaci´on para disoluciones con (a) 0.25 %, (b) 0.1 %, (c) 0.05 %, (d) 0.025 %de CNTS y (e) epoxi. ........ 6 2.1. Dos modelos de dumbbells: (a) FENE dumbbell con menor movilidad para CNTS agregados y (b) FENE dumbbell con mayor movilidad para CNTS que no presentan agregaci´on. . . 8 2.2. Viscosidad aparente (ηa) para (a) epoxi, (b) disoluci´on con 0.25 %de CNTS tratados qu´ımicamente y (c) el resultado de la simulaci´on. .................. 13 3.1. Representaci´on esquem´atica del proceso de spin-coating. ............... 15 4.1. Triangulaci´on de Delaunay y diagrama de Voronoi de una nube de puntos. ...... 18 4.2. Definici´on de las coordenadas de un Vecino Natural de un punto x........... 19 4.3. T´ıpica funci´on φ(x)................................ 19 5.1. Representaci´on de una fibra con su vector orientaci´on ρ................ 22 5.2. Ajuste del modelo de orientaci´on para suspensiones con 0.05 %. 0.2 % y 0.33 % de CNTS tratados qu´ımicamente en una resina epoxi. Npy Drson los par´ametros de ajuste. . . . 25 5.3. Campo de orientaciones para la suspensi´on de CNTS en el plano r-z, para distintos pasos de simulaci´on. Pasos t=1, 10, 20, 30 y 40. ..................... 27 5.4. Detalle del campo de orientaciones para la suspensi´on de CNTS en el plano r-z para el paso de simulaci´on t=1. .............................. 28 V Tanto los nanotubos single-walled como los multi-walled pueden ser tratados qu´ımicamente para que exista repulsi´on electrost´atica entre ellos. Este hecho es muy importante a la hora del comportamiento reol´ogico que pueden presentar las suspensiones de nanotubos de carbono y que ser´a explicado a continuaci´on. 1.1. Reolog´ıa de las suspensiones de nanotubos de carbono Los nanotubos de carbono (tanto los single-walled como los multi-walled) pueden ser tratados qu´ımicamente con agentes sulfatantes para conseguir una repulsi´on electrost´atica entre ellos. Esta repulsi´on electrost´atica va a influir a la hora del comportamiento reol´ogico de la suspensi´on. 1.1.1. La formaci´on de bandas helicoidales en suspensiones de CNTS no tratados qu´ımicamente. En suspensiones de resina epoxi con baja viscosidad con CNTS no tratados qu´ımicamente y mediante la aplicaci´on de un flujo a cortadura se ha podido observar un fen´omeno mediante el cual los nanotubos forman unas estructuras anis´otropas llamadas “bandas helicoidales” (“helical bands”, HBS en ingl´es). En la Figura 1.2 se puede observar la evoluci´on de la microestructura de la suspensi´on para diferentes tiempos y distancia entre placas de un re´ometro. El ritmo de cizalladura y la distancia entre placas son los par´ametros que m´as influyen en la formaci´on de las bandas helicoidales. Los CNTS, que en un primer momento son estructuras is´otropas, bajo el flujo a cortadura al que est´an sometidos poseer´an distinto gradiente de velocidades dependiendo de la profundidad a la que se encuentren en el fluido. Esto har´a que colisionen entre ellos formando unos agregados. A medida que estos agregados choquen unos con otros, se formar´an unos nuevos m´as grandes que adem´as presentan las caracter´ısticas de poseer un di´ametro m´as o menos constante y de mostrar una alineaci´on preferencial en la direcci´on perpendicular a la velocidad del fluido. Se ha comprobado que la formaci´on de bandas helicoidales tambi´en est´a presente en otro tipo de suspensiones. 1.1.2. Reolog´ıa experimental de las suspensiones de CNTS. Cuando a un resina epoxi se le a˜nade nanotubos, tanto tratados como no tratados qu´ımicamente, se observa un comportamiento pseudopl´astico del compuesto. Como puede observarse en la Figura 1.3, la matriz de epoxi posee una viscosidad constante de 10 Pa·s. En el caso de los compuestos con CNTS esta viscosidad es mayor cuando los ritmos de 4 Figura 1.2: Evoluci´on temporal de la microestructura ´optica de una disoluci´on de 0.03 %de MWCNT0s para distintas distancias de placas y un ritmo de cizalladura constante de 0.5 s−1. cizalladura son bajos, disminuyendo a medida que se aumenta el ritmo. Tambi´en se observa como para un mismo porcentaje de CNTS disueltos en la suspensi´on, la viscosidad aparente del compuesto es muy superior para el caso con nanotubos no tratados qu´ımicamente que para el caso de CNTS tratados qu´ımicamente. Este hecho se debe a la formaci´on de bandas helicoidales. La reolog´ıa para las suspensiones con CNTS no tratados qu´ımicamente ha sido estudiada m´as detalladamente y en la Figura 1.4 se muestra la evoluci´on de la viscosidad en funci´on de la velocidad de deformaci´on para cuatro suspensiones con distinta concentraci´on de nanotubos. Las suspensiones con nanotubos no tratados qu´ımicamente muestran un aumento de la viscosidad para bajos ritmos de cizalladura, lleg´andose a incrementar en un orden de magnitud con la adici´on de un 0.1 % de CNTS (Figura 1.4 b). Se puede observar como para ritmos de cizalladura muy altos la viscosidad de las suspensiones tiende asint´oticamente al valor de la viscosidad de la matriz de epoxi. 5 Figura 1.3: Viscosidad aparente (ηa) para (a) epoxi, (b) disoluci´on con 0.25 %de CNTS tratados y (c) disoluci´on con 0.25 %de CNTS no tratados qu´ımicamente. Figura 1.4: Viscosidad aparente (ηa) en funci´on de la velocidad de deformaci´on para disoluciones con (a) 0.25 %, (b) 0.1 %, (c) 0.05 %, (d) 0.025 %de CNTS y (e) epoxi. 6 Cap´ıtulo 2 Un modelo de tipo FENE Las soluciones polim´ericas con CNTS presentan un comportamiento no-Newtoniano, lo que significa que la relaci´on entre el tensor de tensiones y el gradiente de velocidades no puede ser descrito por las ecuaciones de Navier-Stokes. Una alternativa al estudio de estos fluidos a escala macrosc´opica es el desarrollo de modelos basados en la teor´ıa cin´etica que describan los movimientos, a nivel molecular, en los fluidos con base polim´erica y CNTS. Es habitual que muchos de estos modelos microestructurales lleven a simulaciones de Din´amica Browniana, que no s´olo explican la din´amica de las soluciones en estado estacionario sino tambi´en en estados transitorios. El modelo de teor´ıa cin´etica m´as simple para una soluci´on diluida de pol´ımeros lineamente flexibles consiste en dos masas unidas por un muelle lineal (que sigue la ley de Hook) suspendido en un fluido Newtoniano incompresible. De ahora en adelante a este conjunto de masas y muelle pasar´a a denominarse “dumbbell” (por su nombre en ingl´es)[9]. Las masas representan segmentos moleculares de varios mon´omeros y el muelle los efectos de la entrop´ıa. La consideraci´on de que el muelle es lineal (Hookeano) s´olo es realista para peque˜nas deformaciones cerca del equilibrio (distribuci´on Gaussiana) y no pone restricci´on en la elongaci´on que la dumbbell puede tener. La idea de introducir un l´ımite en la elongaci´on responde por un lado a la correcci´on de este comportamiento f´ısicamente imposible y a la importancia que tienen los fen´omenos no lineales en la reolog´ıa. Se ha visto anteriormente que los CNTS pueden presentar distinto grado de agregaci´on lo cual tambi´en debe ser tenido en cuenta en el modelo, es por ello que cuanto mayor sea el estado de agregaci´on de los nanotubos, menor ser´a la movilidad que se le confiera a la dumbbell (Figura 2.1). En este trabajo se supone que los CNTS no presentan agregaci´on ninguna, es decir, que han sido tratados qu´ımicamente. A continuaci´on se presenta una modificaci´on del modelo FENE cl´asico (del ingl´es Finitely Extensible Non-linear Elastic) donde esta limitaci´on de la elongaci´on es tenida en cuenta. 7 Figura 2.1: Dos modelos de dumbbells: (a) FENE dumbbell con menor movilidad para CNTS agregados y (b) FENE dumbbell con mayor movilidad para CNTS que no presentan agregaci´on. 2.1. Introducci´on. La fuerza del muelle en el modelo FENE ser´a: Fc=H 1−Q2 Q2 0 Q(2.1) donde Hes la constante del muelle, Qes el vector tridimensional que une las masas y Q0 es la elongaci´on m´axima permitida del muelle. Se observa en la Ecuaci´on 2.1, que para peque˜nas extensiones se mantiene el comportamiento lineal del muelle y que en el l´ımite de elongaci´on obtendr´ıamos una fuerza infinita. El precio a pagar por esta limitaci´on es la no existencia de una ecuaci´on constitutiva cerrada para el tensor de tensiones del pol´ımero y por lo tanto, no es posible una soluci´on anal´ıtica. 2.2. Formulaci´on del modelo. Se va a considerar un flujo homog´eneo e incompresible, por lo que el campo de velocidades del fluido puede ser escrito como v=κ·rdonde κ= (Ov)T(tensor gradiente de velocidades) es independiente del vector posici´on r[10][11][12]. 2.2.1. Din´amica del pol´ımero. La ecuaci´on de Fokker-Planck (o ecuaci´on de difusi´on) para la funci´on de distribuci´on ψ(Q, t) para el modelo de dumbbell mostrado en la Ecuaci´on 2.1 es: ∂ψ ∂Q=−∂ ∂tκ·Qψ−2 ζFcψ+2kBT ζ ∂ ∂Q ∂ ∂Qψ, (2.2) 8 donde ζes el coeficiente de fricci´on, kBla constante de Boltzmann y Tla temperatura. El primer t´ermino de la derecha de la ecuaci´on proviene de la fuerza de arrastre hidrodin´amica (debida al movimiento de las dumbbells en el solvente), el segundo se refiere a la fuerza del muelle y el ´ultimo se debe a las fuerza Brownianas en las masas (causado por las fluctuaciones t´ermicas en el solvente). El coeficiente de difusi´on kBT/2ζse supone constante. Para una funci´on B(Q) del vector conector, la media de la configuraci´on espacial con dependecia temporal, est´a dada por: hBi=ZB(Q)ψ(Q, t)d3Q(2.3) Por lo tanto, de la ecuaci´on de Fokker-Planck (Ecuaci´on 2.2) se obtiene la siguiente ecuaci´on de movimiento: d dthQQi − κhQQi−hQQiκT=4kBT ζI−4 ζhQFci(2.4) donde Ies la matriz identidad. En el equilibrio (estado estacionario con κ=0), se obtiene: hQFcieq =kBTI(2.5) reflejando que en el equilibrio, la solucion de la Ecuaci´on 2.3 es una distribuci´on de Boltzmann ψeq ∼eφc/kBT, con un potencial φcdefinido por Fc=∂φc ∂Q. Para este modelo de dumbbell se pueden definir dos constantes de tiempo. La primera, λH=ζ 4H, ser´ıa la misma que para un modelo con un muelle Hookeano y la segunda, λQ=ζQ2 0 12kBT, ser´ıa an´aloga a la anterior pero para dumbbells completamente r´ıgidas. Se ha trabajado con la primera de las constantes de tiempo y con un par´ametro b, adimensional, definido como: b=3λQ λH =HQ2 0 kBT. Este par´ametro es llamado coeficiente adimensional de extensi´on del muelle y representa, de manera aproximada, el n´umero de cadenas monom´ericas modelado por la dumbbell en cuesti´on. Esto quiere decir que el valor de bno puede ser cualquiera. Habitualmente, se toman b= 20 como m´ınimo (para que tenga significado f´ısico) y b= 100 como m´aximo (equivalente al modelo de Oldroyd-B, que no ser´a explicado aqu´ı, pero que supone un comportamiento lineal del muelle). 9 2.2.2. Tensor de tensiones. El tensor de tensiones τde la soluci´on polim´erica se supone sim´etrico y resulta de la suma de la contribuci´on del solvente Newtoniano y del pol´ımero, τ=τs+τp=ηs˙γ+τp, donde ηses la viscosidad del solvente y ˙γ=κ+κTes el tensor velocidad de deformaci´on. 2.3. Funciones del material para un flujo a cortadura. Para un flujo a cortadura (que es el caso que nos ocupa) el campo de velocidades viene dado por: v=   ˙γ·y 0 0   donde ˙γes el ritmo de cizalladura y puede ser dependiente del tiempo. La matriz traspuesta del gradiente de velocidades ser´a: κ(t) = ˙γ(t)  0 1 0 0 0 0 0 0 0   y el tensor de tensiones tendr´a la forma: τ=  τxx τxy 0 τyx τyy 0 0 0 τzz   con τxy =τyx, es decir, que el tensor τes sim´etrico. Para un flujo a cortadura, las funciones del material que tienen inter´es de estudio son la viscosidad η, el primer coeficiente normal de tensiones Ψ1y el segundo coeficiente normal de tensiones Ψ2, definidos (de forma adimensional) como sigue: η(˙γ)−ηs nkBTλH =−1 λH˙γ τp,xy nkBT,(2.6) Ψ1(˙γ) nkBTλ2 H =−1 (λH˙γ)2 τp,xx −τp,yy nkBT,(2.7) Ψ2(˙γ) nkBTλ2 H =−1 (λH˙γ)2 τp,yy −τp,zz nkBT.(2.8) Para t= 0, suponemos que el fluido est´a en equilibrio y que el tensor de tensiones τes igual a cero. Para tiempos t > 0 se aplica un ritmo de cortadura ˙γconstante que hace que el valor de las tensiones aumente hasta sus valores estacionarios. Para ritmos de cizalladura suficientemente altos, las funciones del material pueden alcanzar un m´aximo para despu´es aproximarse a un valor constante. 10 2.4. El m´etodo num´erico. Como se ha visto anteriormente, la ecuaci´on de Fokker-Planck del modelo FENE (Ecuaci´on 2.2) no es lineal, por lo que no puede ser resuelta de manera anal´ıtica. La evaluaci´on de las medias en el tensor de tensiones se har´a mediante Din´amica Browniana [19]. Mediante la introducci´on de unidades de longitud (kBT/H)1/2y de tiempo λH= ζ/4H, obtenemos la forma adimensional de: - el vector que une las dos masas de la dumbbell: ˆ Q=Q (kBt/H)1/2, - el tiempo: ˆ t=t λH , - el gradiente de velocidades traspuesto: ˆ κ=κλH. De esta manera, podemos escribir el tensor de tensiones de forma adimensional como: τp nkBT=−*ˆ Qˆ Q 1−ˆ Q2/b++I. Con todos estos par´ametros adimensionales y definiendo una funci´on de distribuci´on, tambi´en adimensional, ˆ ψ(ˆ Q, t) = ψ(Q, t)(kBT/H)3/2, la ecuaci´on de Fokker-Planck (Ecuaci´on 2.2) adimensional ser´ıa: ∂ˆ ψ ∂ˆ t=−∂ ∂ˆ Q· ˆ κ·ˆ Qˆ ψ−1 2 1 1−ˆ Q2/b ˆ Qˆ ψ!+1 2 ∂ ∂ˆ Q ∂ ∂ˆ Q ˆ ψ, (2.9) que es equivalente a la ecuaci´on diferencial estoc´astica de Itˆo para un proceso de Markov tridimensional ˆ Qˆ t: dˆ Qˆ t= ˆ κ·ˆ Qˆ t−1 2 1 1−ˆ Q2 ˆ t/b ˆ Qˆ t!dˆ t+dWˆ t,(2.10) donde Wˆ tes un proceso de Wiener tridimensional. Esta equivalencia significa que la evoluci´on de la densidad de probabilidad, que caracteriza la distribuci´on continua del proceso ˆ Qˆ t (proceso de Markov tridimensional) obtenida resolviendo la Ecuaci´on 2.10, est´a gobernada por la ecuaci´on diferencial descrita en la Ecuaci´on 2.9. Debido a la no linealidad de la ecuaci´on diferencial estoc´astica de Itˆo (Ecuaci´on 2.10) necesitaremos de un m´etodo num´erico 11 para su resoluci´on. El elegido es un algoritmo predictor-corrector semi-impl´ıcito de segundo orden [13]: ˆ Q 0 ˆ tj+1 =ˆ Qˆ tj+ ˆ κˆ tj·ˆ Qˆ tj−1 2 1 1−ˆ Q2 ˆ tj/b ˆ Qˆ tj ∆ˆ t+ ∆Wj(2.11) y  1 + 1 4 ∆ˆ t 1−ˆ Q2 ˆ tj+1 /b ˆ Qˆ tj+1 =ˆ Qˆ tj+1 2 ˆ κˆ tj+1 ·ˆ Q 0 ˆ tj+1 +ˆ κˆ tj·ˆ Qˆ tj−1 2 1 1−ˆ Q2 ˆ tj/b ˆ Qˆ tj ∆ˆ t+∆Wj, (2.12) donde ˆ tj=j∆ˆ t(j= 1,2,3, ...)y∆Wj=Wˆ tj+1 −Wˆ tj, cuyas componentes son variables aleatorias Gaussianas independientes con media cero y varianza ∆ˆ t. 2.4.1. Resultados de las simulaciones. A continuaci´on se muestran algunos de los resultados obtenidos de las simulaciones realizadas mediante la implementaci´on de un c´odigo en MATLAB. Los par´ametros elegidos para la obtenci´on de estos resultados han sido: - n´umero de dumbbels = 50000 - n´umero de pasos de simulaci´on = 15000 - incremento de tiempo = 0.001 -b=20 -λ=2.5 En la Figura 2.2 se muestra la aproximaci´on (punteado rojo) de la evoluci´on de la viscosidad aparente de una disoluci´on con 0.25 % de nanotubos de carbono tratados qu´ımicamente en epoxi. Se puede observar c´omo los resultados se ajustan perfectamente a los datos obtenidos experimentamente tanto en valores como en tendencia. Para otras proporciones de disoluciones de nanotubos de carbono habr´ıa que realizar una nueva b´usqueda de los par´ametros byλque se ajustasen a la correspondiente curva de viscosidad aparente. En el Anexo B se pueden ver m´as resultados. 12 Figura 2.2: Viscosidad aparente (ηa) para (a) epoxi, (b) disoluci´on con 0.25 %de CNTS tratados qu´ımicamente y (c) el resultado de la simulaci´on. 13 20 Cap´ıtulo 5 Dos modelos basados en la teor´ıa cin´etica y la orientaci´on de las fibras El comportamiento de las suspensiones de nanotubos de carbono depende fundamentalmente de si estos han sufrido un tratamiento qu´ımico o no, por lo que los modelos num´ericos no ser´an los mismos. 5.1. Modelo de Fokker-Plack Esta secci´on se centra en suspensiones con CNTS tratados qu´ımicamente, es decir, aquellos que sufren repulsi´on electrost´atica y no muestran tendencia a formar agregados (bandas helicoidales) con el tiempo. A continuaci´on se explica un modelo basado en la teor´ıa cin´etica para pol´ımeros reforzados con fibras cortas. 5.1.1. El modelo num´erico. Para el modelado de una suspensi´on de nanotubos de carbono (CNTS), se ha considerado que se comporta como una suspensi´on de fibras elipsoidales (Figura 5.1) en el seno de un fluido newtoniano. El flujo arrastra a las fibras, que se suponen indeformables. Las ecuaciones que describen el modelo son las siguientes: -Balance de momentos, donde s´olo se consideran las fuerzas centr´ıfugas: Divσ=−ρω×(ω×r) = fr donde σdenota la matriz de tensiones, ρes la densidad,×es el producto tensorial y ωes la velocidad de rotaci´on del plato (spin-coating). 21 Figura 5.1: Representaci´on de una fibra con su vector orientaci´on ρ. -Condici´on de incompresibilidad: Divv= 0 donde ves el campo de velocidades. - La ecuaci´on de comportamiento, con una relaci´on cuadr´atica de cierre para la matriz de orientaci´on de cuarto orden: σ=−pI+ 2η{D+Nptr(a·D)a} donde pdenota la presi´on, Ies la matriz identidad, ηes la viscosidad equivalente de la suspensi´on, Des la matriz de deformaci´on, Npes un par´ametro escalar que depende de la forma y la concentraci´on de los CNTS, tr es la traza y aes la matriz de orientaci´on de segundo orden, definida por: a=Iρ⊗ρΨ(ρ)dρ donde ρdenota el vector unitario de direcci´on de los CNTS, ⊗es el producto vectorial yΨ(ρ) es la funci´on de distribuci´on de la orientaci´on, que debe satisfacer la condici´on de normalidad: IΨ(ρ)dρ= 1. Si Ψ(ρ) = δ(ρ−ˆρ) , con δ() la funci´on Delta de Dirac, toda la probabilidad de orientaci´on se concentra en la direcci´on definida por ˆρ, y la correspondiente matriz de orientaci´on ser´a ˆa=ˆρ⊗ˆρ. En este caso la relaci´on cuadr´atica de cierre resulta exacta. Cuando los nanotubos no est´an perfectamente alineados en una cierta direcci´on, la relaci´on cuadr´atica de cierre se convierte en una aproximaci´on aunque es ampliamente aceptada seg´un la bibliograf´ıa. Desde un punto de vista f´ısico, los valores propios de la matriz de orientaci´on de segundo orden (a) representan la probabilidad de encontrar nanotubos de carbono en la direcci´on correspondiente al vector propio. 22 - Con una relaci´on cuadr´atica de cierre, la ecuaci´on de orientaci´on se expresa como: da dt =Ω·a−a·Ω+k(D·a+a·D−2tr(a·D)a)−6Dra−I 3 DyΩson, respectivamente, las matrices sim´etrica y antisim´etrica de las componentes del Gradv,kes una constante que depende de la relaci´on de forma de los nanotubos (longitud y di´ametro), siendo r=l/ρ, el par´ametro kse define como: k= (r2−1)/(r2+ 1) con k≈1, en la mayor´ıa de los casos. Dres un coeficiente rotacional de difusi´on que tiene en cuenta sucesos como la interacci´on de las fibras en la suspensi´on. El modelo del flujo es definido en el volumen que ocupa el fluido en el tiempo t, Ω(t). En el contorno, Γf(t)≡∂Ωf(t) o la velocidad o la tensi´on (presi´on) son impuestas: v(x∈Γ1) = vg o bien σn(x∈Γ2) = Fg con Γ1∪Γ2=Γf(t),Γ1∩Γ2=∅, y donde n(x) es el vector normal unitario, definido en el contorno en el punto x. En t´erminos de condiciones iniciales, se asume que los nanotubos se encuentran orientados de manera perfectamente is´otropa (es decir, a=I/3). En este modelo, el vector de orientaci´on utilizado para describir la suspensi´on de los CNTS es: ρ=  r θ z  con θ∈[0,2π]. El campo de velocidades es expresado en coordenadas cil´ındricas, despu´es la imposici´on de la condici´on de simetr´ıa axial: v=  vr ω·r vz (5.1) Si se supone que el sistema de referencia rota con el disco, se tendr´a que: v=  vr 0 vz  En este caso la expresi´on de los tensores DyΩen coordenadas cil´ındricas ser´a: 23 D=  ∂vr ∂r 1 2(1 r ∂vr ∂θ +∂vθ ∂r−vθ r)1 2(∂vr ∂z+∂vz ∂r) 1 2(1 r ∂vr ∂θ +∂vθ ∂r−vθ r)1 r ∂vθ ∂θ +vr r 1 2(1 r ∂vz ∂θ +∂vθ ∂z) 1 2(1 r ∂vr ∂z+∂vz ∂r)1 2(∂vz ∂θ +∂vθ ∂z)∂vz ∂z   y siempre que el flujo sea axialmente sim´etrico (vθ= 0,vr=vr(r,z),vz=vz(r,z)), y combin´andolo con la Ecuaci´on 5.1, la matriz Dqueda reducida a lo siguiente: D=  ∂vr ∂r01 2(∂vr ∂z+∂vz ∂r) 0vr r0 1 2(∂vr ∂z+∂vz ∂r) 0 ∂vz ∂z   Similarmente, la matriz de vorticidad tiene la siguiente forma: Ω=  0 0 1 2(∂vr ∂z+∂vz ∂r) 0 0 0 −1 2(∂vr ∂z+∂vz ∂r) 0 0   Como condiciones iniciales en el primer paso de simulaci´on se tomar´a: - La orientaci´on is´otropa de las fibras es definida en el modelo mediante la distribuci´on uniforme: ψ(ρ) = 1 4π - El vector orientaci´on para distribuciones is´otropas: a=1 3I 5.1.2. Relaci´on entre modelo y reolog´ıa El modelo utilizado para la descripci´on del comportamiento de las suspensiones de nanotubos de carbono tratados qu´ımicamente ha sido el de Fokker-Planck (FP). En este modelo aparecen dos par´ametros que son escogidos de manera experimental para que este modelo se ajuste lo mejor posible a los resultados experimentales (DryNp). As´ı, en la siguiente gr´afica (Figura 5.2), podemos observar, para tres concentraciones distintas de nanotubos, la evoluci´on de la viscosidad aparente. Tomando como ejemplo la suspensi´on con 0.33 % de CNTS se comprueba que la pareja de constantes que mejor ajusta la curva es Dr= 0,005s−1 y Np= 7. Para suspensiones con nanotubos que no est´en tratados qu´ımicamente, este modelo de Fokker-Planck no es capaz de explicar su comportamiento, por lo que se necesita un nuevo modelo llamado de Agregaci´on/Orientaci´on (AO), como una ampliaci´on del modelo de Fokker-Planck. Este modelo ser´a explicado en la siguiente secci´on del trabajo. 24 Figura 5.2: Ajuste del modelo de orientaci´on para suspensiones con 0.05 %. 0.2 % y 0.33 % de CNTS tratados qu´ımicamente en una resina epoxi. Npy Drson los par´ametros de ajuste. 5.1.3. Discretizaci´on del problema. Anteriormente se han visto las bases del M´etodo de los Elementos Naturales y se ha comentado c´omo, debido a las caracter´ısticas del problema, este era el m´etodo de interpolaci´on m´as adecuado. Para la discretizaci´on de la forma fuerte del problema, se ha utilizado una interpolaci´on de elementos naturales C0− C−1, m´as concretamente se ha utilizado una interpolaci´on C0(interpolante de Sibson, suave en todo punto excepto en los nodos, φIen la Ecuaci´on 5.2) para la aproximaci´on del campo de velocidades, mientras que para el caso de la aproximaci´on de la presi´on, la interpolaci´on utilizada ha sido una discontinua C−1(ψI en la Ecuaci´on 5.3) : vh(x) = n X I=1 φI(x)vI(5.2) ph(x) = n X I=1 ψI(x)pI= n X I=1 1 npI(5.3) donde vIypIson las velocidades y presiones nodales, respectivamente, y nes el n´umero de vecinos naturales del punto xconsiderados para la interpolaci´on. Este tipo de aproximaci´on no verifica la condici´on LBB (Ladyzhenskaya-Babuˆska-Brezzi). De esta manera tenemos que, con el dominio fluido Ωfextra´ıdo de la nube de nodos mediante el m´etodo de las formas α(Anexo D) y la velocidad y presi´on (definidas por las ecuaciones anteriores), se puede proceder a la discretizaci´on est´andar de la formulaci´on variacional de las ecuaciones del flujo: 25 ZΩf(t) σ:D∗dΩ = ZΩf(t) frv∗dΩ (5.4) ZΩf(t) Divvp∗dΩ = 0 (5.5) y σ=−pI+ 2η{D+NpTr(aD)a}(5.6) donde se asume una tracci´on nula en el frente del fluido y la velocidad impuesta para los nodos en contacto con el plato rotatorio es cero. La ecuaci´on de orientaci´on es resuelta para cada incremento de tiempo. Con la cin´etica conocida en el tiempo t,vt(x), la posici´on de los nodos puede ser actualizada al mismo tiempo que la evoluci´on de la orientaci´on de las fibras. Las ecuaciones de actualizaci´on de la posici´on y de la matriz de orientaci´on de las fibras, de manera expl´ıcita, son las siguientes: xt+∆t I=xt I+vt I∆t,∀I(5.7) y at+∆t I=at I+Ωt Iat I−at IΩt I+kDt Iat I+kat IDt I−2kTr(at IDt I)at I−6Drat I−I 3∆t,∀I (5.8) donde DtIyΩtIson las componentes sim´etrica y antisim´etrica de la matriz del gradiente de velocidades, en el tiempo ten el nodo xI, respectivamente. En la etapa de resoluci´on de la cinem´atica, la matriz de orientaci´on de las fibras se supone conocida en los nodos, para el tiempo t(at I). El valor de aen los puntos de integraci´on para evaluar las Ecuaciones 5.4 y 5.5, se obtiene por interpolaci´on (elementos naturales): at(x) = n X I=1 φI(x)at I La ´unica dificultad a la hora de resolver la Ecuaci´on 5.8 para actualizar la orientaci´on es la no derivabilidad de las funciones de forma en sus nodos. As´ı, el tensor del gradiente de velocidades se puede evaluar mediante la siguiente expresi´on: vt ,k(x) = n X I=1 φI,k(x)vt I(5.9) donde la kcon la coma denota la derivada espacial con respecto a la coordenada k-´esima. El gradiente de velocidades puede ser calculado en todo punto excepto en los nodos, ya que φI,k(xI) no est´a definida. Una posible soluci´on a este problema es el uso de t´ecnicas que realicen una proyecci´on desde los puntos de integraci´on y despu´es se realice un promedio. 26 5.1.4. Resultados de las simulaciones. Para la obtenci´on de los resultados num´ericos que simulan el comportamiento de suspensiones de nanotubos de carbono tratados qu´ımicamente en una resina epoxi, se ha implementado un programa en C basado en el M´etodo de los Elementos Naturales explicado anteriormente y en el modelo de Fokker-Planck. El modelo est´a compuesto por 977 nodos bajo condiciones de simetr´ıa axisim´etrica para el flujo (pero no para el campo de orientaciones), condiciones de no deslizamiento en el plano horizontal, simetr´ıa en condiciones de contorno respecto del eje de simetr´ıa y una distribuci´on de orientaciones inicial is´otropa. Para entender los resultados que a continuaci´on se presentan hay que tener en cuenta que la orientaci´on de los nanotubos se representa mediante elipses: cuando m´as alejada de la isotrop´ıa sea la distribuci´on de probabilidad de encontrar un CNT en una determinada direcci´on, m´as estirada estar´a la elipse en la direcci´on de la orientaci´on. Figura 5.3: Campo de orientaciones para la suspensi´on de CNTS en el plano r-z, para distintos pasos de simulaci´on. Pasos t=1, 10, 20, 30 y 40. En la Figura 5.3 se puede ver la evoluci´on tanto del fluido como de la orientaci´on de los nanotubos para los pasos de simulaci´on t=1, 10, 20, 30 y40. Como se puede observar 27 para el tiempo t=1, la orientaci´on est´a representada por circunferencias, esto quiere decir que para ese paso de tiempo la orientaci´on es nula, convirti´endose estas circunferencias en elipses a medida que pasa el tiempo, es decir los nanotubos se van orientando en el plano r-z. En las Figuras 5.4 y 5.5 se pueden ver con m´as detalle los pasos t=1 yt=30, respectivamente. La diferencia de orientaci´on que existe entre ambos pasos de tiempo se ve aqu´ı claramente. En el primer caso los nanotubos no tienen ning´un tipo de direcci´on preferente, mientras que tras varios pasos de tiempo, los nanotubos se encuentran fuertemente orientados en la direcci´on perpendicular al flujo, como nos indicaba la teor´ıa. Esta orientaci´on queda representada por la gran excentricidad que presentan las elipses (cuanta m´as excentricidad, m´as orientaci´on). Figura 5.4: Detalle del campo de orientaciones para la suspensi´on de CNTS en el plano r-z para el paso de simulaci´on t=1. Figura 5.5: Detalle del campo de orientaciones para la suspensi´on de CNTS en el plano r-z para el paso de simulaci´on t=30. Al igual que se ha representado el plano r-z, tambi´en se puede obtener una representaci´on de lo que ocurre en el plano r-θpara distintos pasos de tiempo (Figura 5.6). Lo que se puede ver en la Figura 5.6, no es el plano r-θdirectamente, si no el abatimiento de este plano, es por eso que, aparentemente, las elipses aparecen orientadas hacia arriba. Como ya se vio en la representaci´on anterior, a medida que pasa el tiempo, la orientaci´on de los nanotubos es m´as predominante en la direcci´on perpendicular a la direcci´on del flujo. En la Figura 5.7 se observa en m´as detalle la orientaci´on en los pasos de tiempo t=1 y en la Figura 5.8, para t=30. 28 Figura 5.6: Campo de orientaciones para la suspensi´on de CNTS en el plano r-θ, para distintos pasos de simulaci´on. Pasos t=1, 10, 20, 30 y 40. 5.2. El modelo de Agregaci´on/Orientaci´on. Debido a la tendencia que presentan los nanotubos no tratados qu´ımicamente a formar agregados bajo las condiciones de flujo a cortadura, el modelo de Fokker-Plank no sirve para explicar la cin´etica de lo que sucede en las suspensiones en este caso. Es por ello que se necesita un nuevo modelo en el cual se incluya un nuevo par´ametro que tenga en cuenta el grado de agregaci´on que hay presente en el fluido. De esta manera se presenta el modelo de Agregaci´on/Orientaci´on (A/O, Aggregation/Orientation, en ingl´es) como una modificaci´on del modelo de Fokker-Plank en el cual se introduce un par´ametro λque modela el grado de agregaci´on presente en el fluido [2]. 5.2.1. El modelo num´erico. Al tratarse este nuevo modelo de una modificaci´on del anterior, las suposiciones presentadas en el apartado anterior no cambian (fibras cortas, flujo incompresible y presencia 29 con las ya existentes. Otra de las simplificaciones que se ha hecho es la de aceptar que los nanotubos forman agregados desde un primer momento, pero no se ha tenido en cuenta el hecho de que a altas velocidades, estos agregados pueden deshacerse, es decir, llegado el agregado a un tama˜no dado, puede irse desprendiendo y formar as´ı agregados m´as peque˜nos. Esta modificaci´on podr´ıa sumarse a la de establecer una variaci´on de la velocidad a lo largo de la simulaci´on de spin-coating, c´omo ocurre en los procesos industriales de esta t´ecnica. Teniendo en cuenta las aplicaciones industriales que podr´ıan tener los nanotubos de carbono, parece interesante tambi´en tener en cuenta otras t´ecnicas, como la del ink-jet oroll-casting. Esto supone un trabajo mucho m´as amplio y duro, ya que no supondr´ıa simplemente una modificaci´on como en los casos anteriores, sino que habr´ıa que crear un modelo nuevo, as´ı como nuevos estudios reol´ogicos para conocer el comportamiento que tienen las suspensiones bajo estas t´ecnicas. 36 Glosario σ:matriz de tensiones. ρ:densidad. ×:producto tensorial. ω:velocidad de rotaci´on de la placa (spin-coating). fr:fuerzas centr´ıfugas. v:campo de velocidades. I:matriz identidad. D:matriz de deformaci´on. Np:escalar que depende de la forma y la concentraci´on de los nanotubos. a:matriz de orientaci´on de segundo orden. ρ:vector unitario de direcci´on de los nanotubos. ⊗:producto vectorial. Ψ(ρ) : fuci´on de distribuci´on de la orientaci´on. δ() : funci´on Delta de Dirac. D:matriz sim´etrica del Gradv. Ω:matriz antisim´etrica del Gradv. k:constante dependiente de la relaci´on de forma de los nanotubos. Dr:coeficiente rotacional de difusi´on. Fc:fuerza del muelle en el modelo FENE. H:constante del muelle en el modelo FENE. Q:vector tridimensional que une las masas de la dumbbell. 37 Q0:elongaci´on m´axima permitida en el muelle de la dumbbell. κ:tensor gradiente de velocidades. r:vector posici´on. ψ(Q, t) : funci´on de distribuci´on para el modelo FENE. ζ:coeficiente de fricci´on. kB:constante de Boltzmann. T:temperatura. λH:constante de tiempo. λQ:constante de tiempo. b:coeficiente adimensional del muelle. τ:tensor de tensiones. τs:tensor de tensiones del solvente. τp:tensor de tensiones del pol´ımero. ηs:viscosidad del solvente. ˙γ:tensor velocidad de deformaci´on. ˙γ(t) : ritmo de cizalladura. η:viscosidad. Ψ1:primer coeficiente normal de tensiones. Ψ2:segundo coeficiente normal de tensiones. Wˆ t:proceso de Wiener tridimensional. λ:estado de agregaci´on de los nanotubos de carbono.