Relaciones de dispersión

Cadena monoatómica lineal

Como vemos en la figura tenemos N partículas de masa m unidas a N+1 muelles iguales de constante k, cuyos extremos están fijos. La separación de equilibrio entre las partículas es a.

Supongamos que en un instante dado t, la partícula 1 se desplaza x1, la partícula 2 se desplaza x2, ... la partícula i se desplaza xi, etc.

La ecuación del movimiento para la partícula i será entonces

m d 2 x i d t 2 =k( x i x i1 )+k( x i+1 x i ) m d 2 x i d t 2 =k x i1 2k x i +k x i+1

i=1,2,3...N, x0=0, xN+1=0

Supongamos, que el sistema vibra en un modo de frecuencia ω. Cada partícula describirá un M.A.S. de la misma frecuencia ω , pero cuya amplitud Ai vamos a calcular.

xi=Ai·cos(ω t)

Introduciendo esta expresión en la ecuación diferencial que describe el movimiento de cada partícula, obtenemos, la relación de recurrencia.

A i+1 + A i1 = A i ( 2 m k ω 2 ) A 0 =0 A N+1 =0 A i+1 2s A i + A i1 =0,s=1 m 2k ω 2

Vamos a calcular las frecuencias ωn por dos procedimientos

Sistema homogéneo de ecuaciones

Creamos un script para calcular ωn para cualquier valor de N

syms s;
N=5; %número de partículas
matriz=-2*s*eye(N);
for i=2:N-1
    matriz(i,i+1)=1;
    matriz(i,i-1)=1;
end
matriz(1,2)=1;
matriz(N,N-1)=1;
eq=det(matriz);
sol=solve(eq);
w=sqrt(2*(1-sol));
double(w)'
w =
             3^(1/2)
             2^(1/2)
                   1
 (3^(1/2) + 2)^(1/2)
 (2 - 3^(1/2))^(1/2)
 
ans =    1.7321    1.4142    1.0000    1.9319    0.5176

Polinomios de Chebyshev

Tomando A0=0 y A1=1, las relación de recurrencia nos proporciona los coeficientes A2, A3...

A 1 =1 A 2 =2s A 1 A 0 =2s A 3 =2s A 2 A 1 =4 s 2 1 A 4 =2s A 3 A 2 =8 s 3 4s A 5 =2s A 4 A 3 =16 s 4 12 s 2 +1 ..... A N+1 =2s A N A N1 U N =2s U N1 U N2

Si vamos al final de la página los polinomios de Chebyshev de segunda especie veremos que A1 se corresponde con el polinomio U0, A2 con el polinomio U1... AN+1=0 (extremo fijo) con UN(s). Las raíces del polinomio UN(s) nos dan las frecuencias ωn de los modos de vibración que vamos a calcular a continuación

Creamos un script para obtener los coeficientes del polinomio UN(s). Mediante la función roots de MATLAB calculamos sus raíces y a continuación, las frecuencias ωn. Supondremos que k/m=1

N=5; %número de partículas
A=zeros(N,1);
A0=[0];
A1=[1];
for i=2:N+1
    A=2*[A1,0]-[0,A0];
    A0=[0,A1];
    A1=A;
end
sol=roots(A);
w=sqrt(2*(1-sol));
disp(w')
w =    0.5176    1.0000    1.4142    1.7321    1.9319

Alternativamente, calculamos las raíces del polinomio de Chebyshev UN(s) a través de su definición,

U N (s)= sin( (N+1)θ ) sinθ =0 θ n = nπ N+1 n=1,2,3...N s n =cos( nπ N+1 ) 1 m 2k ω n 2 = s n ω n =2 k m sin( nπ N+1 )

Las frecuencias ωn para N=5 partículas

N=5;
n=1:N;
w=2*sin(n*pi/(2*(N+1)))
w =    0.5176    1.0000    1.4142    1.7321    1.9319
Amplitudes

En un sistema de N partículas de masa m unidas a N+1 muelles elásticos de constante k, los extremos fijos tienen amplitudes A0=0 y AN+1=0

Para cada frecuencia ωn (n=1,2,3...N) calculamos las amplitudes: A1, A2,A3...AN o bien, los valores de los polinomios de Chebyshev U0(s)=1, U1(s), U2(s),...UN-1(s) para cada valor de sn.

A i = U i1 ( s n )= sin( i· θ n ) sin θ n = sin( i· nπ N+1 ) sin( nπ N+1 ) i=1,2,3..N

Una vez calculadas las amplitudes Ai, se multiplican por un factor para que la suma de sus cuadrados sea la unidad

i=1 N A i 2 =1

Finalmente, queda multiplicar la amplitud de cada partícula i, Aicos(ωnt) para observar el modo n de vibración de las N partículas unidas a N+1 muelles elásticos, tal como se ve en el programa interactivo de la página titulada Modos normales de vibración (I)

Suponiendo que k/m=1, representamos los dos primeros modos de vibración para un sistema de N=5 partículas. En el título de cada gráfica se proporcionan las frecuencias ω1 y ω2. Cambiando el índice n de 1 ... 5 se representan otros modos de vibración

Comprobamos en la página titulada Modos normales de vibración (I) que los valores de las frecuencias coinciden con las obtenidos por otros procedimientos (matricial)

N=5; %número de partículas
wn=2*sin((1:N)*pi/(2*N+2)); %frecuencias

subplot(2,1,1)
n=1; %primer modo de vibración
i=1:N;  %partículas
A=sin(i*n*pi/(N+1))/sin(n*pi/(N+1)); %U_0, U_1, ... U_N-1
suma=sum(A.^2);
A=A/sqrt(suma);
plot(0:N+1,[0,A,0],'-ro','markersize',4,'markeredgecolor','r','markerfacecolor','r')
xlim([0,N+1])
grid on
xlabel('n')
ylabel('A')
title(num2str(wn(n))) %frecuencia

subplot(2,1,2)
n=2; %segundo modo de vibración
i=1:N; %partículas
A=sin(i*n*pi/(N+1))/sin(n*pi/(N+1)); %valor del polinomio U_i correspondiente a w_n
suma=sum(A.^2);
A=A/sqrt(suma);
plot(0:N+1,[0,A,0],'-ro','markersize',4,'markeredgecolor','r','markerfacecolor','r')
xlim([0,N+1])
grid on
xlabel('n')
ylabel('A')
title(num2str(wn(n))) %frecuencia

Se representan los dos primeros modos de vibración para un sistema de N=20 partículas

Amortiguamiento

A la ecuación diferencial que describe el movimiento de la partícula i, estudiada al principio de esta página, añadimos un término que corresponde a una fuerza de rozamiento Fr=λdxi/dt proporcional a la velocidad de la partícula i, (véase las oscilaciones amortiguadas)

m d 2 x i d t 2 =λ d x i dt k( x i x i1 )+k( x i+1 x i ) d 2 x i d t 2 +2γ d x i dt + ω 0 2 ( 2 x i x i1 x i+1 )=0

donde ω 0 2 = k m ,2γ= λ m . Probamos una solución de la forma

x i (t)= A i exp(γt)cos(ωt)

Calculamos la derivada primera y segunda de xi

d x i dt =exp(γt) A i ( γcos(ωt)+ωsin(ωt) ) d 2 x i d t 2 =exp(γt) A i { ( γ 2 ω 2 )cos(ωt)+2γωsin(ωt) }

El resultado es la relación de recurrencia que hemos estudiado al principo de esta página, pero con otra expresión para la variable s

( ω 2 + γ 2 ) A i + ω 0 2 ( 2 A i A i1 A i+1 )=0 A i+1 2s A i + A i1 =0,s= 2 ω 0 2 ω 2 γ 2 2 ω 0 2

con i=1,2,3...N, siendo N el número de partículas de masa m unidas a N+1 muelles elásticos de constante k. Los extremos x0=0 y xN+1=0 están fijos por lo que A0=0 y AN+1=0 que equivale a la ecuación UN(s)=0, tal como se ha explicado en el primer apartado.

Las frecuencias ωn están relacionadas con las raíces sn del polinomio UN(s)=0 de Chebyshev.

s n =cos( nπ N+1 )n=1,2,...N ω n 2 =2 ω 0 2 ( 1cos( nπ N+1 ) ) γ 2 ω n 2 =4 ω 0 2 sin 2 ( nπ 2(N+1) ) γ 2

Calculamos las amplitudes Ai de la partícula i, tal como se hizo en el apartado, titulado Amplitudes. La posición xi de la partícula i en función del tiempo t, suponiendo oscilaciones amortiguadas, es

x i (t)= A 0 sin( i· nπ N+1 ) sin( nπ N+1 ) exp( γt )cos( ω n t )

Siempre que la constante de amortiguamiento γ no sea muy grande, de modo que ωn sea un número real, estaríamos entonces en el caso de oscilaciones amortiguadas, si fuera complejo, daría lugar a oscilaciones sobreamortiguadas

Actividades

Se ha fijado

Se introduce

Se pulsa el botón titulado Nuevo

Observamos el primer modo normal de vibración, en la parte superior izquierda se proporciona la frecuencia angular.

Se pulsa el botón titulado >>, observamos los siguientes modos normales

Se pulsa el botón titulado << para observar el modo normal de vibración anterior.

En la parte inferior, se muestra el desplazamiento de las partículas


Relación de dispersión, N=∞

Cadena monoatómica

Para N→∞, probamos una solución de la forma

Ai=A·sin(K·ia)

donde K es el número de onda K=. La relación de recurrencia se escribe

A i+1 2s A i + A i1 =0 Asin( K(i+1)a )2( 1 m 2k ω 2 )Asin( Kia )+Asin( K(i1)a )=0

Simplificando

2cos(Ka)=( 2 m k ω 2 ) ω 2 =4 k m sin 2 ( Ka 2 )

Esta ecuación que relaciona la frecuencia angular ω con el número de onda K, se denomina relación de dispersión.

En la figura, se representa el cuadrado de la frecuencia angular ω en función de Ka (donde K es el número de onda) en el intervalo (-π, +π).

fplot('4*sin(x/2)^2',[-pi,pi])
set(gca,'XTick',(-1:1)*pi)
set(gca,'XTickLabel',{'-\pi','0','\pi'})
xlabel('Ka')
ylabel('\omega^2')
grid on
title('Relación de dispersión')

Cadena diatómica

Consideremos una cadena lineal de moléculas diatómicas separados una distancia a en la situación de equilibrio, tal como se muestra en la figura.

Consideremos el movimiento de dos átomos contiguos.

Ecuaciones del movimiento

M d 2 y i d t 2 =k( y i x i1 )+k( x i+1 y i ) m d 2 x i+1 d t 2 =k( x i+1 y i )+k( y i+2 x i+1 )

Probamos una solución de la forma

yi=Aicos(ωt), con Ai=Asin(Kia)

xi+1=Ai+1cos(ωt), con Ai+1=Bsin(K(i+1)a)

donde K se denomina número de onda

Introducimos estas soluciones en las dos ecuaciones diferenciales y obtenemos un sistema homogéneo de dos ecuaciones con dos incógnitas A y B.

A(-ω2M+2k)sin(Kia)-Bk(sin(K(i+1)a)+sin(K(i-1)a))=0
-Ak
(sin(K(i+2)a)+sin(Kia))+B(-ω2m+2k) sin(K(i+1)a)=0

Teniendo en cuenta que

sin(K(i+1)a)+sin(K(i-1)a)=2sin(Kia)cos(Ka)
sin(K(i+2)a)+sin(Kia)= 2sin(K (i+1)a)cos(Ka)

y eliminamos A y B en el sistema homogéneo de dos ecuaciones

(-ω2M+2k) (-ω2m+2k)sen(Kia)sen(K(i+1)a)-4k2sen(Kia)cos2(Ka)sen(K(i+1)a))=0

Simplificando el factor común sin(Kia)sin(K(i+1)a)

ω 4 2k m+M mM +4 k 2 mM sin 2 (Ka)=0 ω 2 =k( 1 M + 1 m )±k ( ( 1 M + 1 m ) 2 4 sin 2 (Ka) mM )

Esta ecuación que relaciona la frecuencia ω con el número de onda K, se denomina relación de dispersión.

En la figura, se representa se representa el cuadrado de la la frecuencia angular ω en función Ka (donde K es del número de onda) en el intervalo (-π/2, +π/2). La curva superior (signo +) se denomina rama óptica, y la inferior (signo -) se denomina rama acústica.

M=0.75;
m=1;
k=1;
f1=@(x) k*(1/m+1/M)+k*sqrt((1/m+1/M)^2-4*sin(x)^2/(m*M));
f2=@(x) k*(1/m+1/M)-k*sqrt((1/m+1/M)^2-4*sin(x)^2/(m*M));
hold on
fplot(f1,[-pi/2,pi/2])
fplot(f2,[-pi/2,pi/2])
hold off
set(gca,'XTick',(-1:1)*pi/2)
set(gca,'XTickLabel',{'-\pi/2','0','\pi/2'})
xlabel('Ka')
ylabel('\omega^2')
grid on
title('Relación de dispersión')

Referencias

Runk R. B. Stul J. L. Anderson G. L. A laboratory analog for lattice dynamics. Am. J. Phys. (31) 1963, pp. 915-921

John H. Van Drie, Jr., Peter W. Murphy. Chebyshev polynomials in the loaded-string problem. American Journal of Physics Vol. 43, No. 4, April 1975, pp. 361-362

N Gauthier. Exact expressions for the amplitudes and damped normal frequencies of taut vibrating linear string of couped masses. Eur. J. Phys. 29 (2008) N21-N29