Dinámica de gases diluidos

Grado en F´ısica TRABAJO DE FIN DE GRADO Din´amica de gases diluidos Javier Caro Ruiz 4 de julio de 2024 Dirigido por Mar´ıa Isabel Garc´ıa de Soria Lucena y Pablo Maynar Blanco Facultad de F´ısica de la Universidad de Sevilla ´ Indice 1. Introducci´on 1 2. Modelo de esferas duras 3 3. Funci´on de distribuci´on 7 3.1. Ecuaci´on de balance para la funci´on de distribuci´on . . . . . . . . . . . . . . 8 3.1.1. T´erminodeflujo ............................. 9 3.1.2. T´ermino de colisiones . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2. Teorema HdeBoltzmann ............................ 15 3.3. Tendenciaalequilibrio .............................. 17 4. Simulaci´on 23 4.1. Par´ametros de entrada e inicializaci´on . . . . . . . . . . . . . . . . . . . . . . 24 4.2. Evoluci´ontemporal................................ 26 4.2.1. Subrutina1:Flujo ............................ 27 4.2.2. Subrutina 2: Tiempo m´ınimo de colisi´on . . . . . . . . . . . . . . . . 27 4.3. Resoluci´on de la colisi´on y exportaci´on de variables . . . . . . . . . . . . . . 29 5. An´alisis de la simulaci´on 31 5.1. Distribuci´onplana ................................ 32 5.2. Distribuci´on sim´etrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.3. Distribuci´on asim´etrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 6. Conclusiones 44 6.1. Modelo de esferas duras y funci´on de distribuci´on . . . . . . . . . . . . . . . 44 6.2. Simulaci´on y an´alisis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 1. Introducci´on Todo juego consta de dos factores fundamentales: dos o m´as elementos bien definidos y una serie de reglas que rigen su comportamiento y c´omo se relacionan. Por ejemplo, en Polic´ıas y ladrones existen dos tipos de elementos bien diferenciados, los polic´ıas y los ladrones, y una serie de reglas que los relacionan, como que los polic´ıas deben atrapar a todos los ladrones para ganar el juego, que los ladrones pueden ayudar a sus compa˜neros a salir del calabozo, etc. La F´ısica trata de dar una explicaci´on a nivel fundamental de todo lo que nos rodea con un esquema id´entico al de Polic´ıas y ladrones: una serie de elementos y un conjunto de reglas que los relacionan. La necesidad de estos elementos bien definidos que conforman nuestra realidad es la que, tras siglos de confrontaci´on de ideas, llega al concepto de ´atomo (y, m´as tarde, a las part´ıculas fundamentales), mientras que las leyes de la F´ısica toman el papel de las reglas. Suponiendo que lleg´aramos a coleccionar el conjunto completo de reglas, lo m´as razonable ser´ıa pensar que, efectivamente, somos capaces de comprender todo lo que nos rodea. A esta corriente de pensamiento se la denomina Determinismo, y era aceptada por la gran mayor´ıa de cient´ıficos contempor´aneos con la figura central de este trabajo: Ludwig Boltzmann. Sin embargo, aunque l´ogica, esta idea presenta una dificultad pr´actica fundamental. Las ecuaciones de Newton describen adecuadamente el movimiento de dos cuerpos. Cuando a˜nadimos un tercero, sin embargo, obtener una soluci´on anal´ıtica de las ecuaciones de movimiento se vuelve casi imposible. Si tratamos un gas con decenas de ´ordenes de magnitud de cuerpos, las ecuaciones son, sencillamente, irresolubles anal´ıticamente. Es justamente este el problema que presentan los gases (entre otros sistemas f´ısicos) en la descripci´on determinista. Se puede describir la din´amica de un gas con las ecuaciones de Newton, pero la gran cantidad de ecuaciones diferenciales que se deben resolver se traduce en una inviabilidad pr´actica. En este contexto es en el que se desarrolla la Termodin´amica, que, en lugar de considerar un sistema como un conjunto de mol´eculas, crea un formalismo matem´atico en el que se incluyen cantidades del sistema como la energ´ıa libre, la entrop´ıa, la temperatura, etc. y se determina c´omo deben actuar en el equilibrio. Boltzmann, por otro lado, decide partir del sistema como un conjunto muy numeroso de part´ıculas, pero con la gran idea de aplicar una maquinaria 1 estad´ıstica en la que condensa todas las ecuaciones de Newton de las part´ıculas en una sola ecuaci´on integrodiferencial, como veremos en este trabajo. Su estudio, sin embargo fue completamente rechazado por sus contempor´aneos, que no terminaron de comprender c´omo la descripci´on estad´ıstica no romp´ıa con el determinismo. A´un as´ı, Boltzmann defendi´o sus ideas hasta el d´ıa de su muerte y, algunos a˜nos m´as tarde, uno de sus mayores detractores, Max Planck, utilizar´ıa su trabajo para comenzar a desarrollar lo que m´as tarde ser´ıa la Mec´anica Cu´antica. En este trabajo deduciremos la ecuaci´on de Boltzmann que describe la evoluci´on temporal de la funci´on de distribuci´on, as´ı como su conocido teorema H, que demuestra de forma rigurosa que cierta funcional siempre crece y que podemos identificarla con la entrop´ıa. Adem´as, tras desarrollar un programa que permite simular un gas de esferas duras, probaremos diferentes configuraciones iniciales y analizaremos c´omo se ajustan a los resultados de Boltzmann. 2 2. Modelo de esferas duras Consideremos un sistema de N esferas duras de di´ametro σy masa mque se mueven bajo la acci´on de una fuerza externa F=−∇ϕy que colisionan el´asticamente. Analicemos en detalle c´omo es la colisi´on entre dos part´ıculas con posiciones ryr1y velocidades vyv1, respectivamente. Conviene reescribir estas posiciones y velocidades en funci´on de un par de variables que nos faciliten los c´alculos. Generalmente, se trabaja en el sistema centro de masas y usando las posiciones y velocidades relativas, R=r−r1;V=v−v1; (1) rcm =r+r1 2;vcm =v+v1 2.(2) Una vez definidas estas variables, nos interesa buscar magnitudes que se conserven para el estudio de las colisiones. Partiremos de dos leyes de conservaci´on fundamentales en las colisiones el´asticas: la del momento lineal y la de la energ´ıa cin´etica. La formulaci´on de la primera es p=p′;  mv+  mv1=  mv′+  mv′ 1,(3) donde, aplic´andola a nuestro caso en el que la masa de todas las part´ıculas es id´entica, ves la velocidad inicial de la part´ıcula 1, v1la de la part´ıcula 2, v′es la velocidad final de la part´ıcula 1 y v′ 1es la velocidad final de la part´ıcula 2. La expresi´on de la ley de conservaci´on de la energ´ıa cin´etica ser´a T=T′;  1 2mv2+  1 2mv2 1=  1 2mv′2+  1 2mv′ 1 2.(4) Ahora bien, podr´ıamos reescribir los vectores posici´on y velocidad de las part´ıculas en funci´on de estas nuevas cantidades, esto es, cambiar las coordenadas del espacio de fases. Se obtienen como sigue, r=rcm +1 2R,v=vcm +1 2V; (5) r1=rcm −1 2R,v1=vcm −1 2V.(6) Con esta transformaci´on obtenemos un resultado interesante, debido a que estas nuevas variables no dependen del origen de coordenadas, sino ´unicamente de las caracter´ısticas de la 3 colisi´on en particular. Por tanto, las leyes de conservaci´on enunciadas anteriormente deben imponer alg´un tipo de ley de conservaci´on en estas nuevas variables. Veamos qu´e obtenemos si aplicamos el nuevo cambio de variables, vcm +  1 2V+vcm −  1 2V=v′ cm +  1 2V′+vcm ′−  1 2V′.(7) De la conservaci´on del momento lineal obtenemos que la velocidad del centro de masas se conserva en las colisiones, algo que ya esper´abamos por la propia definici´on de la cantidad. Veamos qu´e obtenemos al aplicar el cambio de variables en la ecuaci´on de conservaci´on de la energ´ıa, vcm +1 2V2 +vcm −1 2V2 =vcm ′+1 2V′2 +vcm ′−1 2V′2 , 2vcm2+ 21 2V 2 = 2vcm ′2+ 21 2V′ 2(8) Y, a consecuencia de la ec. (7), llegamos a |V|=|V′|,(9) esto es, que el m´odulo de la velocidad relativa de las part´ıculas no cambia tras la colisi´on. No obstante, debemos subrayar que la orientaci´on de la velocidad relativa s´ı que debe cambiar tras el choque en general. Dado que las colisiones entre esferas duras se producen en un ´unico plano, denotaremos con θel ´angulo de desviaci´on de la velocidad relativa. Teniendo esto en cuenta, la colisi´on de dos esferas duras se traducir´a en un cambio instant´aneo de las velocidades de las mismas. El ´unico factor por determinar, por tanto, ser´a el ´angulo de giro de la velocidad relativa, aunque es posible hallar el cambio en las velocidades sin determinar directamente este ´angulo. En primer lugar, definiremos como σel vector que une los centros de las esferas en el instante del choque, cuyo m´odulo ser´a el di´ametro de las mismas y, su vector unitario, ˆ σ=σ σ. Asumiendo una colisi´on el´astica, el cambio de velocidad, si lo hay, debe producirse en la direcci´on de ˆ σ, pues si no hay rozamiento la ´unica fuerza que act´ua sobre las part´ıculas es la propia de la colisi´on, que s´olo se ejerce en la direcci´on mencionada. Es decir, que el cambio de las velocidades de las esferas debe ser proporcional a la proyecci´on de la velocidad relativa 4 Ven la direcci´on ˆ σ, v′=v+k1(V·ˆ σ)ˆ σ, v′ 1=v1+k2(V·ˆ σ)ˆ σ. (10) Aplicando la conservaci´on del momento lineal, v′+v′ 1=v+v1, v+k1(V·ˆ σ)ˆ σ+v1+k2(V·ˆ σ)ˆ σ=v+v1. (11) Como el producto escalar no es nulo de forma general, esta igualdad implica que k1+k2= 0, es decir, k1=−k2=k. Incluyendo la conservaci´on de la energ´ıa, que es cin´etica en su totalidad, v′2+v2 1=v2+v2 1, v2+ 2k(V·ˆ σ)ˆ σ·v+k2(V·ˆ σ)2+v2 i−2k(V·ˆ σ)ˆ σ·v1+ (−k)2(V·ˆ σ)2=v2+v2 1 (12) Despejando lleva a la condici´on, k2+k= 0,(13) por lo que kdebe ser nula o igual a −1. La soluci´on nula es trivial: corresponde a que no se transfiera energ´ıa en la colisi´on, es decir, que esta no ocurra y las esferas se atraviesen. Por tanto, llegamos a que k=−1. De esta forma, el cambio de velocidades debido a una colisi´on de dos esferas duras ser´a, v′=v−(V·ˆ σ)ˆ σ, v′ 1=v1+ (V·ˆ σ)ˆ σ. (14) Podemos reescribir este cambio de variables con un operador bσ, v′≡bσv, v′ 1≡bσv1. (15) Es decir, que para cualquier funci´on ϕ(v,v1), este operador act´ua de tal forma que bσϕ(v,v1) = ϕ(bσv, bσv1) = ϕ(v′,v′ 1). Nos ser´a ´util calcular el jacobiano de esta transformaci´on para c´alculos posteriores. Este se escribe como J=∂(v′,v′ 1) ∂(v,v1)=∂(v′,v′ 1) ∂(V′,v′ cm) ∂(V′,v′ cm) ∂(V,vcm) ∂(V,vcm) ∂(v,v1).(16) 5 Comencemos calculando la inversa del tercer t´ermino, cuyas relaciones entre variables propusimos al comienzo de esta secci´on, ∂(v,v1) ∂(V,vcm)=∂(vx, vx1, vy, vy1, vz, vz1) ∂(Vx, vcm,x, Vy, vcm,y, Vz, vcm,z)= =  −1 21 0 0 0 0 1 21 0 0 0 0 0 0 −1 21 0 0 0 0 1 21 0 0 0 0 0 0 −1 21 0 0 0 0 1 21  =−1. (17) Ahora bien, el primer t´ermino de la ec. (16) es id´entico al calculado, ya que el cambio de variables de las velocidades finales a los vectores velocidad centro de masa y relativa finales es el mismo. Por tanto, obtenemos que, ∂(v′,v′ 1) ∂(V′,v′ cm )=−1. Para calcular el segundo t´ermino, el cambio de la velocidad del centro de masas y la relativa se puede obtener aplicando sus correspondientes f´ormulas a la ec. (14) obteniendo v′ cm =vcm V′=V−2(V·ˆ σ)ˆ σ. (18) Y, efectivamente, nos queda ∂(V′,v′ cm) ∂(V,vcm)=∂V′ ∂V ∂v′ cm ∂vcm = 1 ·1=1.(19) Finalmente, tenemos que J=∂(v′,v′ 1) ∂(v,v1)= 1.(20) 6 3. Funci´on de distribuci´on Cuando nos enfrentamos a un problema f´ısico en el que tenemos un gran n´umero de part´ıculas interaccionando entre ellas, el estudio de sus colisiones individuales resulta impracticable anal´ıticamente. A´un con la ayuda de ordenadores, la complejidad de los c´alculos escala de forma desmedida con el n´umero de part´ıculas, por lo que los m´etodos estad´ısticos pasan a ser el centro de atenci´on. No obstante, los m´etodos num´ericos nos ser´an muy ´utiles para poner de manifiesto las conclusiones obtenidas en el estudio estad´ıstico. En el trabajo veremos que se puede obtener una ecuaci´on de evoluci´on cerrada para la funci´on de distribuci´on de velocidades, f(r,v;t), definida de forma que f(r,v;t)drdves el n´umero medio de part´ıculas que, en un determinado instante t, se encuentran en el elemento de volumen drcentrado en r, con una velocidad comprendida en el elemento de volumen dv centrado en v. Cabe destacar que, cuando nos refiramos al valor medio, estaremos hablando de un promedio sobre condiciones iniciales distintas. Adem´as, para poder obtener una descripci´on cerrada en t´erminos de fser´a necesario admitir que esta var´ıa en una escala tal que es posible utilizar el c´alculo diferencial en elementos de volumen que contienen multitud de part´ıculas. Teniendo esto en cuenta, relacionamos la funci´on de distribuci´on con el n´umero total de part´ıculas del sistemas, N, como, N=Zdrdvf(r,v;t).(21) Si trabajamos con un gas de esferas duras (esto es, que las interacciones son instant´aneas), las propiedades referidas al comportamiento del mismo s´olo van a requerir informaci´on patente en la funci´on de distribuci´on, suponiendo que las part´ıculas son indistinguibles. Por ejemplo, una cantidad interesante es la densidad local de part´ıculas, n(r;t), es decir, el n´umero medio de part´ıculas que se encuentran entre ry drpara un instante t. Como cabe esperar, basta contar todas las part´ıculas que se encuentren en dicho intervalo sin considerar la velocidad de las mismas, n(r;t) = Zdvf(r,v;t).(22) 7 Figura 2: Representaci´on de unas posibles paredes el´asticas (l´ınea curva) del recipiente. vrepresenta la velocidad inicial de la part´ıcula, v′la velocidad de la part´ıcula tras la colisi´on con la pared, y N(r) el vector unitario normal a la pared en el punto de colisi´on. Si consideramos el´astica la pared de este recipiente, las part´ıculas interact´uan con ella reflejando su velocidad. Teniendo en cuenta las cantidades mostradas en la Figura 2, v′≡bNv−2(v·N(r))N(r).(45) La funci´on de distribuci´on, por tanto, experimentar´a un cambio de la forma f(r,v;t)dv= f(r, bNv;t)dv′. Sin embargo, al ser una reflexi´on ∂v′ ∂v= 1. Nos queda, entonces, que la condici´on de contorno para paredes el´asticas de la ecuaci´on de Boltzmann es, f(r,v;t) = f(r, bNv;t)∀r∈∂Υ; ∀v v v∈R,(46) siendo ∂Υ toda la superficie de las paredes. Evidentemente, podr´ıamos establecer diferentes condiciones de contorno para modelar otro tipo de paredes, pero en este trabajo ´unicamente trataremos las el´asticas. Es f´acil ver que la distribuci´on de equilibrio introducida anteriormente en la ec. (26) satisface la ecuaci´on de Boltzmann con estas condiciones de contorno. Aunque la hayamos determinado para esferas duras, la ecuaci´on de Boltzmann es mucho m´as general. No s´olo es aplicable a potenciales de interacci´on entre part´ıculas arbitrarios (siempre que sean de corto alcance), sino que tambi´en es capaz de modelar diferentes tipos de part´ıculas en un mismo gas (mezclas), esferas duras inel´asticas, etc. Particularmente, puede deducirse de forma an´aloga para discos duros, sistema que simularemos en la siguiente secci´on. Cabe destacar adem´as que se ha realizado una gran simplificaci´on a nivel de la descripci´on del sistema. Hemos pasado de un modelo para un gas con 3N ecuaciones diferenciales ordinarias 14 (leyes de Newton para cada part´ıcula), condiciones iniciales e interacci´on con las paredes a una ´unica ecuaci´on integrodiferencial con condiciones iniciales y de contorno. Aunque no es f´acil obtener soluciones expl´ıcitas de la ecuaci´on, veremos que se pueden extraer importantes conclusiones f´ısicas. 3.2. Teorema Hde Boltzmann Ahora que tenemos una descripci´on cerrada que rige la din´amica de f(r,v;t) para bajas densidades, nos planteamos si esta se adec´ua a la realidad y, en particular, a los resultados de la Termodin´amica. Cabe destacar que esto no es algo trivial: durante la obtenci´on de la ecuaci´on de balance, hemos realizado una serie de suposiciones que, en principio, no tendr´ıan por qu´e ser v´alidas. Existen dos resultados fundamentales de la Termodin´amica que deben estar reflejados en la ecuaci´on de balance: la obtenci´on de un estado de equilibrio para tiempos largos y, suponiendo un sistema aislado, el aumento de la entrop´ıa hasta un cierto valor m´aximo en el propio estado de equilibrio. El reflejo de ambas propiedades se materializa en el teorema Hde Boltzmann, que pasamos a analizar. Definimos la siguiente funcional de f H(t)≡ZdrZdvf(r,v;t) ln f(r,v;t).(47) Vamos a demostrar que, si fverifica la ecuaci´on de Boltzmann, entonces dH(t) dt≤0. Esto es el denominado teorema Hde Boltzmann. Efectivamente, dH(t) dt=ZdrZdv∂f ∂t ln f+∂f ∂t .(48) Sin embargo, acudiendo a la ecuaci´on 21, ZdrZdv∂f ∂t =d dtZdrdvf=dN dt= 0,(49) lo que nos deja un s´olo t´ermino. Aplicando la ecuaci´on de balance (27), podemos dividir el c´alculo de este en dos, uno de flujo y otro de colisiones. Comencemos con el de flujo, Zdrdv∂f ∂t flujo ln f=−Zdrdvv·∂f ∂r+F m·∂f ∂vln f. (50) 15 Dividimos esta integral en dos t´erminos. A˜nadiendo que ∂ ∂r vfln f=v∂f ∂rln f+∂f ∂r, podemos reescribir lo anterior como, Zdrdv v ·∂f ∂rln f=Zdrdv∂ ∂r·[vfln f−f].(51) Aplicando el teorema de Gauss, y teniendo en cuenta que dado un sistema aislado, el flujo de cualquier cantidad en sus fronteras debe ser nulo por las condiciones de contorno, llegamos a Zdrdv v ·∂f ∂rln f=ZdvZS dS·vfln f−f= 0.(52) Para el segundo t´ermino, utilizamos ∂ ∂v·hF(r) mfln fi=F(r) m ∂f ∂vln f+F(r) m ∂f ∂v, ZdrdvF m·∂f ∂vln f=ZdrZdv∂ ∂vF(r) mfln f−F(r) mf= 0,(53) donde hemos aplicado el teorema de Gauss en el ´ultimo paso. Por tanto, s´olo el t´ermino de colisiones ser´a relevante para la dependencia temporal de la cantidad H. Desarroll´andolo, dH(t) dt=σ2 2ZdrZdvZdv1Zdˆ σ|V·ˆ σ|(bσ−1)ff1ln f, (54) Entendiendo fcomo f(r,v;t) y f1como f(r,v1;t), de la misma forma que m´as adelante utilizaremos an´alogamente f′yf′ 1. Podemos reescribir esta expresi´on de tres formas distintas. En primer lugar, intercambiemos las variables vyv1, dH(t) dt=σ2 2ZdrZdvZdv1Zdˆ σ|V·ˆ σ|(bσ−1)ff1ln f1,(55) que resulta equivalente a (54) salvo por el t´ermino del logaritmo. En segundo lugar, trataremos la colisi´on inversa, es decir, (v′,v′ 1)→(v,v1). La cantidad Hpara esta se definir´a como dH(t) dt=σ2 2ZdrZdv′Zdv′ 1Zdˆ σ|V′·ˆ σ|(1 −bσ)ff1ln f′.(56) Pero, si tenemos en cuenta la ec. (20), podemos volver a integrar sobre las variables anteriores, dv′,dv′ 1= dvdv1. Adem´as, por la ec. (18), sabemos que |V′·ˆ σ|=|−V·ˆ σ|, por lo que podemos reescribir lo anterior como dH(t) dt=−σ2 2ZdrZdvZdv1Zdˆ σ|V·ˆ σ|(bσ−1)ff1ln f′.(57) 16 Por ´ultimo, haciendo el cambio de variable v→v1,v1→ven la colisi´on inversa, tenemos que dH(t) dt=−σ2 2ZdrZdvZdv1Zdˆ σ|V·ˆ σ|(bσ−1)ff1ln f′ 1.(58) Sumando estas cuatro expresiones y despejando la derivada temporal de H, obtenemos dH(t) dt=−σ2 2ZdrZdvZdv1Zdˆ σ|V·ˆ σ|(f′f′ 1−ff1) ln f′f′ 1 ff1 .(59) Ahora bien, el producto [f′f′ 1−ff1] ln f′f′ 1 ff1da necesariamente un valor positivo o nulo, ya que, si uno de los t´erminos es positivo, el otro es tambi´en lo es y viceversa. Como las funciones de distribuci´on fson positivas o nulas, dH(t) dt≤0.(60) Precisamente, este es el nexo de uni´on que busc´abamos con la Termodin´amica. Al encontrar una cantidad que disminuye necesariamente con el transcurso del tiempo podemos relacionarla con la entrop´ıa, aunque a´un no podamos obtener una relaci´on clara. M´as adelante las relacionaremos rigurosamente. Adem´as, la ec. (59) nos dice que ∂H ∂t = 0 si y s´olo si ff1=f′f′ 1 ∀r,v,v1,ˆ σ. 3.3. Tendencia al equilibrio La necesidad de un estado de equilibrio impl´ıcita en el teorema Hes ligeramente m´as compleja de encontrar. Desde luego, si existe un estado de equilibrio, las cantidades que promedien sobre todo el sistema, como ser´ıa H(t), deben ser constantes una vez se alcancen tiempos largos, es decir, dH(t→ ∞)/dt= 0. Ya vimos que esto implica directamente f′ 1f′=f1f, (61) es decir, que las colisiones entre part´ıculas no afecten a la funci´on de distribuci´on. En otras palabras, estamos buscando una cota inferior para el valor de Hpues, si conocemos que siempre disminuye, necesitamos asegurar que s´olo puede hacerlo hasta un cierto valor y no de forma indefinida. Trataremos esto para un sistema con paredes el´asticas sobre el que act´ua una fuerza F(r)=−∇ϕ(r). 17 Para demostrar que, efectivamente, existe esta cota, partiremos de la desigualdad matem´atica xln x≥x−1∀x≥0.(62) Si tomamos x=f/feq , definiendo feq como la funci´on de distribuci´on en el caso en que realmente el sistema llegue al equilibrio, es decir, la distribuci´on de Maxwell-Boltzmann feq =neq(r)β π3/2exp{−βv2}con neq =Aexp−2β mϕyβ=m 2kT , se cumple fln f≥fln feq +f−feq (63) Si integramos ahora los miembros de cada lado y recordamos que tanto fcomo feq est´an normalizadas a N, ZdrZdvfln f≥ZdrZdvfln feq = =ZdrZdvf3 2lnβ π+ ln A−2β mhm 2v2+ϕ(r)i>−∞, (64) ya que la energ´ıa asociada a fyfeq es equivalente, se cumple ZdrZdvfln feq =ZdrZdvfeq ln feq,(65) o, lo que es lo mismo, H[f]≥H[feq], como quer´ıamos demostrar. Adem´as, si Hs´olo disminuye y est´a acotado, resulta evidente que dH(t→ ∞)/dt= 0, es decir, que existe un estado f∗(r,v;t) que verifica esta igualdad. Como propusimos anteriormente, esta distribuci´on debe cumplir f∗(r,v;t)f∗(r,v1;t) = f∗(r,v′;t)f∗(r,v′ 1;t)∀r,v,v1,ˆ σ(66) Tomando logaritmos, ln f∗(r,v;t) + ln f∗(r,v1;t) = ln f∗(r,v′;t) + ln f∗(r,v′ 1;t).(67) De esta ecuaci´on inferimos que ln f∗(r,v;t) es un invariante aditivo en las colisiones y, por tanto, debe ser combinaci´on lineal de los invariantes aditivos estudiados, es decir, ln f∗(r,v;t) = −α(r;t)−β(r;t)v2−γ(r;t)·v,(68) donde αdebe estar relacionado con el n´umero de part´ıculas, βcon la energ´ıa y αcon el momento lineal. Si deshacemos el logaritmo, f∗(r,v;t) = exp−α(r;t)−β(r;t)v2−γ(r;t)·v.(69) 18 Para determinar las funciones α,βyγvamos, en primer lugar, a reescribir f∗como f∗= exp−α+γ2 4βexp(−βv+γ 2β2).(70) Calculemos ahora la densidad local de part´ıculas, en la que aplicaremos el cambio w=v+γ 2β, n(r;t) = Zdvf∗(r,v;t) = exp−α+γ2 4βZdwexp−βw2= =π β3/2 exp−α+γ2 4β. (71) Por otro lado, la velocidad media local, se obtiene por n(r;t)u(r;t) = Zdv vf∗(r,v;t) = exp−α+γ2 4βZdww−γ 2βexp−βw2= =−n(r;t)γ 2β, (72) de tal forma que u=−γ 2β. Y, de nuevo, podemos reescribir f∗con estas nuevas cantidades, f∗(r,v;t) = n(r;t)β(r;t) π3/2 expβ(r;t)[v−u(r;t)]2.(73) Esta funci´on resulta similar a la funci´on de equilibrio si denotamos β(r;t) = m 2kT(r;t), con la particularidad de que la densidad, velocidad media y temperatura depende de r. Se conoce como funci´on de equilibrio local. Por ´ultimo, nos falta imponer que f∗cumple la ecuaci´on de Boltzmann, ya que ´unicamente la hemos definido como funci´on que sigue la condici´on dH(t→ ∞)/dt= 0. Efectivamente, ya que ∂f∗ ∂t colisiones = 0 por definici´on, ∂f∗ ∂t +v·∂f∗ ∂r+F m·∂f∗ ∂v= 0.(74) Por simplicidad, volveremos a f∗(r,v;t) = exp{−α(r;t)−β(r;t)v2−γ(r;t)·v}. las deri- 19 vadas involucradas pueden describirse de la forma ∂f∗ ∂t =−∂α ∂t +∂β ∂t v2+∂γ ∂t ·vf∗, ∂f∗ ∂xi =−∂α ∂xi +∂β ∂xi v2+∂γ ∂xi·vf∗, v·∂f∗ ∂r=−f∗vi ∂α ∂xi +vi ∂β ∂xi v2+vi ∂γj ∂xi vj, ∂f∗ ∂vi =−2βvi+γif∗, F m·∂f∗ ∂v=−2βv·F m+γ·F mf∗, (75) donde se usa el convenio de suma para ´ındices repetidos. Volviendo a la ecuaci´on de Boltzmann y agrupando en t´erminos de v, se obtienen ∂α ∂t γ+F m+v·∂γ ∂t +∂α ∂r+ 2βF m+v2v·∂β ∂r+v2∂β ∂t +∂γi ∂xi+vivj ∂γj ∂xi = 0,(76) donde el ´ultimo t´ermino es una suma en la que los ´ındices no pueden repetirse (i=j), ya que hemos incluido los elementos en los que i=jen el pen´ultimo. Ahora bien, cada uno de estos t´erminos debe ser igual nulo por separado, ya que la velocidad es arbitraria. Del tercer t´ermino de la ec. (76), ∂β ∂r= 0, obtenemos que β(t) es independiente de r. Si escribimos el pen´ultimo t´ermino para dos dimensiones, tenemos ˙ β(t) + ∂γx ∂x = 0 →∂2γx ∂x∂y = 0 (77) ˙ β(t) + ∂γy ∂y = 0 →∂2γy ∂x∂y = 0.(78) Por la ec. (77) sabemos que ∂γx ∂y es independiente de xy, por (78), que ∂γy ∂x es independiente de y. Adem´as, seg´un el ´ultimo t´ermino de la ec. (76), ∂γx ∂y =−∂γy ∂x =−J(t),(79) que, unido a lo anterior, nos da γ(r;t) = γ0(t)−˙ β(t)r+J(t)×r,(80) donde γ0(t) y J(t) son funciones del tiempo (en particular, Jest´a relacionado con una velocidad angular). Ya que u(r;t) = −1 2β(t)γ(r;t), el campo de velocidades debe ser una 20 combinaci´on lineal de un campo de velocidades constante, uno radial y uno rotacional, sin embargo, al considerar al gas encerrado en una superficie arbitraria, estos tres deben ser nulos por condiciones de contorno (el flujo en las paredes debe ser nulo, ya que las esferas o discos no pueden salir). Este an´alisis se basa en la referencia [3]. Es decir, γ0= 0, J(t) = 0 y˙ β(t) = 0. Esto implica tanto u(r;t) = 0como β= cte. Volviendo a la ec. (76), introduciendo γ= 0 en el primero llegamos a que ∂α ∂t = 0, es decir, α=α(r), mientras que el segundo queda ∂α ∂r+ 2βF(r) m=0→∂ ∂rα−2β mϕ=0,(81) es decir, α(r) = 2β mϕ(r) + K, siendo Kuna constante. Sustituyendo todo lo obtenido, llegamos a f∗(r,v;t) = exp−α(r) + βv2=Aexp−1 kT hm 2v2+ϕ(r)i=feq(r,v),(82) que es justamente lo que quer´ıamos probar: f∗(r,v;t) es soluci´on de la ecuaci´on de Boltzmann si y s´olo si coincide con la con figuraci´on de equilibrio. Hemos demostrado, por tanto, que partiendo de cualquier configuraci´on inicial, la funci´on de distribuci´on tiende a la funci´on de equilibrio correspondiente al NyEdado para la condici´on inicial. De esta forma, tendremos que, para un sistema homog´eneo, aislado (F= 0) y cerrado, la funci´on de distribuci´on de equilibrio ser´a feq(v) = nm 2πkT 3/2exp−mv2 2kT (83) Como propusimos anteriormente, debe existir cierta relaci´on entre la cantidad H(t) y la entrop´ıa S(t), ya que esta crece en todo proceso irreversible, tal y como es la evoluci´on de un gas diluido. Propongamos una definici´on de la entrop´ıa, S(t) = −k H(T) + S0,(84) siendo S0una constante y, por tanto, dS(t) dt=−kdH(t) dt≥0.(85) 21 Ahora bien, si obtenemos la expresi´on de H(t) integrando en (47) con la expresi´on calculada para la ecuaci´on de balance, llegamos a S=−kN ln n−3 2ln T+3 2lnm 2πk−3 2+S0,(86) que, tomando adecuadamente S0, llegamos a la expresi´on de la entrop´ıa para un gas ideal en equilibrio, S=kN 3 2ln N−ln n+3 2ln(2πmk) + 5 2.(87) El argumento para la identificaci´on de H(t) con la entrop´ıa en una y dos dimensiones, aunque la expresi´on sea distinta, es completamente an´alogo: ´unicamente ser´ıa necesario integrar las funciones de distribuci´on pertinentes y compararlas con la entrop´ıa para un sistema unidimensional o bidimensional, seg´un corresponda. 22 4. Simulaci´on Ahora que tenemos un marco te´orico completo para el estudio de un gas diluido de esferas duras, lo pondremos a prueba con una simulaci´on en la que podremos partir de diferentes condiciones iniciales y hacer mediciones en el equilibrio. Partiendo ´unicamente de la ecuaci´on para el cambio de velocidades en una colisi´on (Secci´on 2), el programa hace evolucionar el sistema temporalmente y modifica instant´aneamente las velocidades en caso de colisi´on. Por motivos de tiempo de simulaci´on y capacidad de procesamiento del ordenador a disposici´on, las simulaciones ´unicamente se han realizado para dos dimensiones, aunque el c´odigo es f´acilmente extrapolable a tres dimensiones. Aunque en un principio se escribieron en MATLAB, se pas´o a Fortran 90, ya que MATLAB es mucho menos eficiente para los c´alculos pertinentes. No obstante, el tratamiento de datos y trazado de gr´aficas se ha realizado en MATLAB. El esquema general de la simulaci´on est´a fundamentado en [4]. Para llevar acabo la simulaci´on, hemos supuesto, al igual que en el estudio te´orico, que nuestro gas es homog´eneo e is´otropo a escala mesosc´opica. Por tanto, podremos encontrar un recinto razonablemente grande como para asumir que el gas en general es un conglomerado de copias id´enticas de este (Figura 3(a)). (a) (b) Figura 3: (a) Interpretaci´on f´ısica de las condiciones de contorno peri´odicas del sistema. (b) Diagrama del recinto en el que realizaremos la simulaci´on, aunque con n´umero de part´ıculas simplificado. Las simulaciones, sin embargo, trabajar´an exclusivamente en el recinto de la Figura 3 (b). Este recinto tendr´a condiciones de contorno peri´odicas, de tal forma que funciona tal y como si se tratara de una porci´on de un gas de volumen infinito. Podemos hacer esto gracias a las suposiciones anteriores de homogeneidad e isotrop´ıa. 23 coordenada seg´un corresponda. Tras este apunte se calcula la colisi´on, modific´andose ´unicamente las componentes de las velocidades de las part´ıculas, ya que la interacci´on es instant´anea. Este proceso se repite para el n´umero de colisiones totales a excepci´on de la inicializaci´on y, a su vez, todos los pasos se repiten tantas veces como trayectorias se requieran. Durante y tras el transcurso del programa, seg´un corresponda, calculamos numerosos conjuntos de datos que, posteriormente, se exportar´an en archivos de texto para ser analizados en un programa externo como MATLAB. Entre ellos, se encuentran las distribuciones de velocidades iniciales y finales, as´ı como varias distribuciones intermedias para visualizar la evoluci´on del sistema. Tambi´en se calculan y exportan los momentos de la velocidad (desde ⟨v⟩hasta ⟨v6⟩), ya que es la mejor forma de comprobar que se ha llegado al estado de equilibrio: una vez alcanzan los valores te´oricos. Los momentos tender´an a estos valores, de tal forma que esperamos que coincida aproximadamente la funci´on de distribuci´on correspondiente a la colisi´on en la que los momentos se han estabilizado con una gaussiana. De nuevo, mientras se promedie en m´as trayectorias, mayor suavidad tendr´a el histograma y m´as se asemejar´a a la gaussiana. 30 5. An´alisis de la simulaci´on Una vez desarrollado con ´exito, el programa permite crear simulaciones de gases diluidos con cualesquiera variables de entrada en dos dimensiones. En esta secci´on veremos c´omo, si introducimos valores de densidad bajos (la fijaremos en 0.02 part´ıculas por elemento de superficie), llegamos al estado de equilibrio tal y como predice el teorema H. Adem´as, mediremos la funci´on de distribuci´on de velocidades en este estado y compararemos lo obtenido con los resultados esperados por Boltzmann. Para ello, como comentamos anteriormente, hemos partido de tres distribuciones de velocidades iniciales diferentes, de tal forma que esperamos diferentes ritmos de alcanzar el estado de equilibrio. Para verificar que hemos llegado al mismo nos hemos valido de los momentos de las velocidades, ya que son cantidades que deben tender a un cierto valor de equilibrio y que son sencillas de calcular teniendo los datos de las velocidades de todas las part´ıculas. Utilizaremos el concepto de colisiones por part´ıcula del sistema, s. Los valores te´oricos de los momentos se obtienen adaptando la funci´on de distribuci´on de equilibrio de la secci´on anterior para el caso de dos dimensiones, quedando f(v)eq,2D=nm 2πkT exp−mv2 2kT .(91) La obtenci´on de los momentos de la velocidad de las part´ıculas en el equilibrio se calcula por, ⟨vn⟩eq =1 nZdvvnf(v)2D.(92) Durante toda la simulaci´on, fijamos los valores n= 0,02σ−2,m= 1kg,kT = 1J. Resolviendo las integrales, obtenemos que los valores te´oricos de los momentos son, ⟨v⟩=rπ 2ms−1v2= 2m2s−2v3= 3rπ 2m3s−3 v4= 8m4s−4v5= 15rπ 2m5s−5v6= 48m6s−6 (93) Principalmente, trabajamos con dos formas de visualizar las distribuciones de velocidades: Distribuciones marginales (ec. (25)): Una vez obtenemos las componentes de las velocidades de todas las part´ıculas del sistema en un instante concreto (ya sea inicial, 31 final o tras un cierto n´umero de colisiones) para cada una de las trayectorias realizadas, hacemos un promedio de todas ellas y las representamos en un histograma. No s´olo promediaremos en todas las trayectorias, sino tambi´en en las componentes vxyvy. Podemos hacer esto debido a que el programa construye las condiciones iniciales de las dos componentes de forma sim´etrica, por lo que son independientes salvando la condici´on de que ⟨v2⟩debe ser constante. Distribuciones completas: Por otro lado, representaremos todas las distribuciones obtenidas en diagramas de dos variables, una para cada componente de la velocidad. Estas representaciones son mucho m´as restrictivas que las marginales en el sentido de que, evidentemente, tendremos menos conteos para cada par de coordenadas. Adem´as, como veremos m´as adelante, esta representaci´on nos da una informaci´on m´as precisa de la evoluci´on real del sistema, pues habr´a casos en los que la distribuci´on marginal se asemeje a una gaussiana y la completa diste en gran medida de ella. Como norma general, normalizaremos las distribuciones marginales a N. Las distribuciones completas, ya que tienen un papel meramente representativo, recogen los datos de todas las trayectorias sin ning´un tipo de normalizaci´on. 5.1. Distribuci´on plana Comenzaremos con la distribuci´on m´as sencilla de las tres, en las que las componentes de las velocidades iniciales toman valores aleatorios entre dos extremos de forma independiente. Como podemos observar en la Figura 6, la representaci´on completa (b) muestra una desviaci´on de una distribuci´on plana al uso de forma mucho m´as evidente que la marginal (a). Esto se debe a que, realmente, cada partici´on del histograma de la Figura 6 (a) consiste en la suma de todas las particiones sumadas del eje vypara un mismo valor de vxdel (b). Tras conocer la distribuci´on inicial, nos interesa una estimaci´on de la cantidad de colisiones por part´ıcula que requerir´a para alcanzar el equilibrio. Como comentamos anteriormente, los momentos de las velocidades de las part´ıculas nos ser´an de gran utilidad (Figura 7), ya que podemos observar claramente c´omo alcanzan el valor estacionario. Debemos recalcar que, como esper´abamos, alcanzan los valores te´oricos con gran precisi´on. Aunque en el primer 32 (a) (b) Figura 6: Distribuci´on inicial de velocidades plana para un sistema de 1000 part´ıculas y 200 trayectorias en representaci´on marginal (a) y completa (b). momento, Figura 7 (a), parezca no alcanzarlo, si nos fijamos en el eje de abscisas vemos que la diferencia entre estos valores es de la mil´esima, error que puede entrar en el rango de las imprecisiones de los m´etodos num´ericos. Efectivamente, como podemos observar en la Figura 7, los momentos se estabilizan alrededor de la colisi´on 5000, es decir, para s= 5 colisiones por part´ıcula. Sin embargo, debido a que la representaci´on en histograma llega a asemejarse a una gaussiana mucho antes, en la Figura 8 hemos presentado el sistema tras 400 y 1000 colisiones, es decir, s= 0,4 y s= 1 colisiones por part´ıcula. 33 (a) (b) (c) (d) (e) (f) Figura 7: Momentos de las velocidades (⟨v⟩hasta v6, respectivamente) frente a las colisiones por part´ıcula sacontecidas. En rojo, se marca el valor te´orico en el equilibrio. La tendencia general es la esperada: exceptuando al segundo momento, que debe ser invariante al ser proporcional a la energ´ıa cin´etica total, todos tienden al valor de equilibrio. 34 (a) (b) (c) (d) Figura 8: Distribuciones de velocidades para 1000 part´ıculas y 200 trayectorias, con distribuci´on inicial de velocidades plana. (a) y (b) corresponden a la distribuci´on para s= 0,4 y, (c) y (d) para s= 1. En la Figura 9 observamos la distribuci´on final de velocidades para la distribuci´on plana tras 10000 colisiones, es decir, 10 colisiones por part´ıcula. Sabemos, como comentamos anteriormente, que el sistema ha alcanzado el estado estacionario gracias a la estabilizaci´on de los momentos, por lo que esperamos que esta distribuci´on est´e modelada por la gaussiana f(v) = nm 2πkT exp−mv2 2kT ,(94) fijando m= 1kg ykT = 1J. Adem´as, debido a que trabajamos con histogramas divididos en un cierto n´umero de celdas, vamos a trabajar con la funci´on normalizada, F(v) = 1 √2πexp−v2 2≈0,4 exp−v2 2.(95) Para comprobar c´omo se adec´ua nuestro resultado al estado de equilibrio, hemos normalizado a 1 el histograma de la distribuci´on final (Figura 9 (a)). Efectivamente, los resultados se adec´uan de manera excelente a la predicci´on te´orica. 35 (a) (b) Figura 9: Distribuci´on final de velocidades para una distribuci´on inicial plana, un sistema de 1000 part´ıculas y 200 trayectorias en representaci´on marginal normalizada a 1 (a) y completa (b) tras 10 colisiones por part´ıcula. En rojo (Figura (a)), se muestra la tendencia te´orica esperada para el equilibrio. 5.2. Distribuci´on sim´etrica Continuaremos con la segunda de las distribuciones iniciales mencionadas, la sim´etrica. En este caso, hemos dispuesto los discos en la red de tal forma que las part´ıculas tienen velocidades pr´acticamente diagonales, variando en un peque˜no ´angulo que est´a relacionado con la anchura de la distribuci´on inicial de velocidades (Figura 10 (a)) alrededor de los dos extremos. Este ´angulo es aleatorio para cada disco, de tal forma que es posible que se den las colisiones. De lo contrario, ´unicamente tendr´ıamos un sistema cuyas part´ıculas se mueven diagonalmente con colisiones perfectamente frontales. (a) (b) Figura 10: Distribuci´on inicial de velocidades sim´etrica para un sistema de 1000 part´ıculas y 200 trayectorias en representaci´on marginal (a) y completa (b). 36 Se presentan los momentos de las velocidades ⟨vn⟩en funci´on del n´umero de colisiones por part´ıcula del sistema sen la Figura 11. A diferencia del caso anterior, esta distribuci´on no alcanza el equilibrio hasta pasada la s´eptima colisi´on por part´ıcula. (a) (b) (c) (d) (e) (f) Figura 11: Momentos de las velocidades (⟨v⟩hasta v6, respectivamente) frente a las colisiones acontecidas para una distribuci´on inicial de velocidades sim´etrica. Cabe remarcar la secci´on constante al comienzo de cada momento. Existe una diferencia fundamental con los momentos anteriores en el segmento horizontal por el que comienzan todas las gr´aficas. Este se debe, precisamente, a lo comentado previamente sobre las primeras colisiones del sistema, que son casi frontales. Un sistema compuesto por 37 part´ıculas que s´olo colisionan frontalmente es, en esencia, un conjunto de sistemas unidimensionales. Aunque no se haya tratado en este trabajo, puede demostrarse con facilidad que los momentos de las velocidades en sistemas unidimensionales son constantes, ya que la ley de cambio de velocidades se convierte en una trivialidad. Justamente, esa es la tendencia inicial que observamos. (a) (b) (c) (d) Figura 12: Distribuciones de velocidades para 1000 part´ıculas y 200 trayectorias, con distribuci´on inicial de velocidades sim´etrica. (a) y (b) corresponden a la distribuci´on para s= 0,6 y, (c) y (d) para s= 2. Seguidamente, en la Figura 12 podemos apreciar diferentes instantes en la evoluci´on del sistema. Es realmente interesante notar c´omo en este caso se evidencia la superioridad de la representaci´on completa frente a la marginal: cuando uno podr´ıa pensar que se est´a acercando al equilibrio tras comprobar la Figura 12 (c), pues se comienza a asemejar a una gaussiana, basta con pasar a la Figura 12 (d) para constatar la disparidad entre la distribuci´on actual y la de equilibrio. De nuevo, esto se debe a que la representaci´on marginal est´a relacionada con la proyecci´on de la completa en uno de los ejes. 38 (a) (b) Figura 13: Distribuci´on final de velocidades para una inicial sim´etrica, un sistema de 1000 part´ıculas y 200 trayectorias en representaci´on marginal normalizada a 1 (a) y completa (b) para 10 colisiones por part´ıcula. En rojo (Figura (a)), se muestra la tendencia te´orica esperada para el equilibrio. Finalmente, en la Figura 13 (a) vemos c´omo nuestra distribuci´on de velocidades para s= 10 colisiones por part´ıcula encaja a la perfecci´on con las predicciones te´oricas. Efectivamente, obtenemos una funci´on de distribuci´on id´entica al caso de la distribuci´on inicial plana. 5.3. Distribuci´on asim´etrica Finalmente, trataremos la ´ultima de las tres distribuciones, la asim´etrica. Es muy similar a la anterior, con las particularidades de que tanto el n´umero de part´ıculas como las velocidades difieren para cada extremo. En nuestro caso, la partici´on de part´ıculas ser´a de una en el extremo de alta velocidad a dos en extremo de baja velocidad, cada uno de ellos contando con cierta anchura en la distribuci´on de velocidades. En la Figura 15 representamos la evoluci´on de los seis primeros momentos de las velocidades de nuestro sistema en funci´on del n´umero de colisiones por part´ıcula. Al igual que en la distribuci´on sim´etrica, pero de forma atenuada, vemos una ligera tendencia en las primeras colisiones diferente a la general. Aunque en este caso resulte menos obvio que en el anterior, tambi´en existe cierta inclinaci´on al caso unidimensional, ya que las velocidades son diagonales. Sin embargo, al diferir en m´odulo, este tramo acaba mucho antes. 39 Referencias [1] J. J. Brey Abalo y J. De La Rubia Pacheco. Introducci´on a la mec´anica estad´ıstica. Ediciones del Castillo S.A., 1978. isbn: 84-219-0164-8. [2] P. R´esibois y M. De Leener. Classical Kinetic Theory of Fluids. John Wiley y Sons, 1977. isbn: 978-0471716945. [3] C. Cercignani. The Boltzmann Equation and Its Applications. Springer, 1988. isbn: 978- 0387966373. [4] M. P. Allen y D. J. Tildesley. Computer Simulation of Liquids (2nd edn). Oxford University Press, 2017. isbn: 978-0198803201. 46