Modos de vibración de una membrana circular

Consideremos una membrana circular de radio R, que está fija por su borde. La ecuación del movimiento ondulatorio en coordenadas polares se escribe.

2 Ψ t 2 = c 2 ( 2 Ψ r 2 + 1 r Ψ r + 1 r 2 2 Ψ θ 2 )

Simetría radial

Estudiamos primero, el caso más sencillo, las soluciones Ψ(r,t) independientes del ángulo θ.

2 Ψ t 2 = c 2 ( 2 Ψ r 2 + 1 r Ψ r )

Condiciones iniciales y de contorno

Las condiciones de contorno son

Ψ(R,0)=0t0

Las condiciones iniciales son:

t=0{ Ψ(r,0)= Ψ 0 (r) Ψ(r,t) t | t=0 = Ψ ˙ 0 (r)

Modos de vibración

Aplicamos el método de separación de variables

Ψ(r,t)=Φ(r)·T(t)

Sustituimos en la ecuación diferencial

1 T(t) · d 2 T d t 2 = c 2 Φ(r) ( d 2 Φ d r 2 + 1 r dΦ dr )

Como el miembro izquierdo depende solamente de t y el derecho solamente de r, entonces, igualamos ambos a una constante que denominaremos -ω2.

1 T(t) · d 2 T d t 2 = c 2 Φ(r) ( d 2 Φ d r 2 + 1 r dΦ dr )= ω 2 { d 2 Φ d r 2 + 1 r dΦ dr + ω c 2 2 Φ=0 d 2 T d t 2 + ω 2 T(t)=0

Definimos una nueva variable s=ω·r/c.

d 2 Φ d s 2 + 1 s dΦ ds +Φ=0

Esta es la ecuación de Bessel de orden cero. La solución de esta ecuación es

Φ=C· J 0 (s)+D· Y 0 (s)

Como Y0(s) se hace infinito cuando s se aproxima a cero, el coeficiente D deberá ser cero.

Φ(r)= J 0 ( ω c r )

La condición de contorno implica que para todo t, Ψ(R,t)=0 por lo que Φ(R)=0. Las raíces de la ecuación

J0(α)=0 con α=ω·R/c.

son α1=2.405, α2=5.520, α3=8.654,.... Los sucesivos valores de α determinan las frecuencias ω1, ω2 .... ωn, de los distintos modos de vibración.

Representamos gráficamente J0(α)

>> alfa=linspace(0,20,200);
>> plot(alfa,besselj(0,alfa),'b')
>> xlabel('\alpha')
>> ylabel('y')
>> title('Raíces de J_0')
>> grid on

Con la herramienta Data cursor, determinamos los valores aproximados de las raíces: 2.4, 5.5, 8.6, 11.8, 14.9...

Con la función fzero de MATLAB calculamos las raíces pasándole su valor aproximado.

>> f=@(x) besselj(0,x);
>> for x0=[2.4 5.5 8.6 11.8 14.9]
fzero(f,x0)
end

ans =    2.4048
ans =    5.5201
ans =    8.6537
ans =   11.7915
ans =   14.9309

La diferencia entre dos raíces consecutivas se aproxima a π.

Representamos los primeros modos normales de vibración

alpha=[2.4048,5.5201,8.6537,11.7915,14.9309];
R=1; %radio de la membrana circular

r=0:0.01:1;
ang=0:pi/50:2*pi;
[r,ang]=meshgrid(r,ang);
x=r.*cos(ang);
y=r.*sin(ang);
z=besselj(0,alpha(3)*r/R); %tercer modo
surf(x,y,z)
xlabel('x')
ylabel('y')
zlabel('z')
title('Vibración de una membrana circular')

La solución de la ecuación diferencial dependiente del tiempo es tambien conocida

d 2 T d t 2 + ω 2 ·T=0 T n (t)= A n cos( ω n t)+ B n sin( ω n t)

La solución de la ecuación diferencial correspondiente al modo n de vibración de la membrana circular es

Ψ n (r,t)= Φ n (r) T n (t)= J 0 ( α n R r )( A n cos( ω n t)+ B n sin( ω n t) )

Solución completa

La solución general que satisface las condiciones de contorno es la superposición

Ψ(r,t)= n=1 Ψ n (r,t) = n=1 J 0 ( α n R r )( A n cos( ω n t)+ B n sin( ω n t) )

Esta ecuación describe todos los posibles modos de vibración de la membrana circular. La vibración particular que experimenta la membrana está únicamente determinada por las condiciones iniciales, que a su vez determinan los valores de las constantes An y Bn.

Ψ(r,0)= Ψ 0 (r) Ψ(r)= n=1 J 0 ( α n R r ) A n A n = 2 R 2 J 1 2 ( α n ) 0 R r Ψ 0 (r) J 0 ( α n R r ) drn=1,2,3... Ψ(r,t) t | t=0 = Ψ ˙ 0 (r) Ψ ˙ (r)= n=1 J 0 ( α n R r ) ω n B n B n = 2 R 2 J 1 2 ( α n ) ω n 0 R r Ψ ˙ 0 (r) J 0 ( α n R r ) drn=1,2,3...

Modos de vibración de una membrana circular (solución general)

Consideremos una membrana circular de radio R, que está fija por su borde. La ecuación del movimiento ondulatorio en coordenadas polares se escribe.

2 Ψ t 2 = c 2 ( 2 Ψ r 2 + 1 r Ψ r + 1 r 2 2 Ψ θ 2 )

Aplicamos el método de separación de variables y sustituimos en la ecuación diferencial

Ψ(r,t)=Φ(r,θ)·T(t) 1 T(t) · d 2 T d t 2 = c 2 Φ(r,θ) ( 2 Φ r 2 + 1 r Φ r + 1 r 2 2 Φ θ 2 )= ω 2

Como el miembro izquierdo depende solamente de t y el derecho solamente de r y θ, entonces, igualamos ambos a una constante que denominaremos -ω2.

{ 2 Φ r 2 + 1 r Φ r + 1 r 2 2 Φ θ 2 + ω 2 c 2 Φ=0 d 2 T d t 2 + ω 2 T(t)=0

Aplicamos el método de separación de variables a la primera ecuación diferencial

Φ(r,θ)=R(r)·Θ(θ) Θ(θ)( r 2 d 2 R d r 2 +r dR dr )+R(r) d 2 Θ d θ 2 + ω 2 c 2 r 2 R(r)Θ(θ)=0 1 R(r) ( r 2 d 2 R d r 2 +r dR dr + ω 2 c 2 r 2 R(r) )= 1 Θ(θ) d 2 Θ d θ 2 = m 2

Obteniendo el sistema de ecuaciones diferenciales

{ d 2 R d r 2 + 1 r dR dr +( ω 2 c 2 m 2 r 2 )R=0 d 2 Θ d θ 2 + m 2 Θ=0

La solución de estas ecuaciones diferenciales es

{ R(r)= C r J m ( ω c r )+ D r Y m ( ω c r ) Θ(θ)= C θ cos(mθ)+ D θ sin(mθ)

Donde m es un número entero m=0,1,2,3...

Como Y(r) se hace infinito cuando r se aproxima a cero, el coeficiente Dr deberá ser cero, como en el caso anterior.

Representamos las primeras funciones Jm(x), para m=0,1,2,3

color=['b' 'g' 'r' 'k'];
alfa=linspace(0,20,200);
hold on
for m=0:3
    plot(alfa,besselj(m,alfa),color(m+1),'displayName',num2str(m))
end
hold off
xlabel('\alpha')
ylabel('y')
title('J_m')
legend('-DynamicLegend','location','Best')
grid on

Como la función radial R(r) se hace cero cuando r=R en el borde de la membrana circular. Calculamos las raíces de 1,2,3...n de la función Jm(α) que denominaremos αm,n, Los sucesivos valores de α·R/c determinan las frecuencias ωm1, ωm2 .... ωmn, de los distintos modos de vibración.

x=linspace(0,20,20);
for m=0:3;
    alfa=raices(@(x) besselj(m,x),x);
    disp([m,alfa])
end
    0         2.4048    5.5201    8.6537   11.7915   14.9309   18.0711
    1.0000    3.8317    7.0156   10.1735   13.3237   16.4706   19.6159
    2.0000    5.1356    8.4172   11.6198   14.7960   17.9598
    3.0000    6.3802    9.7610   13.0152   16.2235   19.4094

Nota: la función raices se describe en la página titulada Funciones de Bessel donde se calculan las raíces y se representan las funciones Jm(x), para m=0,1,2..,

J0(x) J1(x) J2(x)
2.4048 3.8317 5.1356
5.5201 7.0156 8.4172
8.6537 10.1735 11.6198
11.7915 13.3237 14.7960
14.9309 16.4706 17.9598
18.0711 19.6159 21.1170

El modo de vibración m,n se describe por la función

Ψ m,n (r,θ,t)= J m ( α m,n r R )·( C n cos(mθ)+ D n sin(mθ) )·( A n cos( ω m,n t)+ B n sin( ω m,n t) ) m=0,1,2,...n=1,2,3...

Los primeros modos de vibración de una membrana circular. Las líneas en color rojo señalan los nodos, los puntos de la membrana que no vibran para un modo particular. Más abajo, veremos la representación en 3D de los distintos modos (m,n) de vibración

Los nodos corresponden a las raíces de la ecuación

J m ( α m,n r R )=0

Consideremos el caso n=2. Supongamos que hemos calculado las raíces (véase la tabla más arriba) α01=5.5201 (segunda raíz de J0(x)), α11=7.0156 (segunda raíz de J1(x)), α21=8.4172 (segunda raíz de J2(x))

alfa=[5.5201,7.0156,8.4172];
for m=0:2
    f=@(x) besselj(m,alfa(m+1)*x);
    r=fzero(f, 0.5);
    disp(r)
end

Radios r/R de los nodos

    0.4356
    0.5462
    0.6101

Consideremos el caso n=3. Supongamos que hemos calculado las raíces (véase la tabla más arriba) α02=8.6537 (tercera raíz de J0(x)), α12=10.1735 (tercera raíz de J1(x)), α22=11.6198 (tercera raíz de J2(x))

alfa=[8.6537,10.1735,11.6198];
for m=0:2
    f=@(x) besselj(m,alfa(m+1)*x);
    r1=fzero(f, 0.2);
    r2=fzero(f, 0.8);
    disp([r1,r2])
end

Radios r/R de los nodos

    0.2779    0.6379
    0.3766    0.6896
    0.4420    0.7244

Representamos el modo de vibración de la membrana circular, para m=1 y para n=3, tercera raíz de J1(x)

Ψ m,n (r,θ,t)= J m ( α m,n r R )·cos(mθ)·cos( ω m,n t)

La función seno se podría poner en vez de coseno.

alpha=10.1735; % raiz n=3 de J1(x) 
R=1; %radio de la membrana circular
r=0:0.01:R;
ang=0:pi/50:2*pi;
[r,ang]=meshgrid(r,ang);
x=r.*cos(ang);
y=r.*sin(ang);
fichero = 'placa_c.gif';
P=1;
m=1; %modo de vibración
n=3;
for t=0:P/20:P
    z = besselj(m,alpha*r/R).*cos(m*ang)*cos(2*pi*t/P);
    surf(x,y,z)
    xlabel('x')
    ylabel('y')
    zlabel('z')
    title('Vibración de una membrana circular')
    axis([-1 1 -1 1 -1 1]);
    
    frame=getframe;
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    if t==0
        imwrite(imind,cm,fichero,'gif','DelayTime',0,'loopcount',inf);
    else
        imwrite(imind,cm,fichero,'gif','DelayTime',0,'writemode','append');
    end
end

Modos de vibración de un anillo

Supongamos una membrana en forma de anillo de radio exterior a y radio interior b

Partimos de la solución de las dos ecuaciones diferenciales, que hemos deducido en el apartado titulado 'Modos de vibración de una membrana circular (solución general)'

{ R(r)= C r J m ( ω c r )+ D r Y m ( ω c r ) Θ(θ)= C θ cos(mθ)+ D θ sin(mθ)

La función R(r) se tiene que hacer cero para r=a y para r=b. Obtenemos un sistema homogéneo de dos ecuaciones con dos incógnitas C y D.

{ C J m ( ω c a )+D Y m ( ω c a )=0 C J m ( ω c b )+D Y m ( ω c b )=0

Despejamos C en la primera ecuación y lo sustituimos en la segunda, obteniendo una ecuación trascendente en x=ω/c

J m ( ω c a ) Y m ( ω c b ) Y m ( ω c a ) J m ( ω c b )=0

a=4; %radios
b=2;
x=linspace(0,10,20);
for m=0:3;
    f=@(x) bessely(m,x*b).*besselj(m,x*a)-besselj(m,x*b).*bessely(m,x*a);
    k=raices(f,x);
    disp([m,k])
end
         0    1.5615    3.1367    4.7091    6.2807    7.8520    9.4231
    1.0000    1.5983    3.1562    4.7222    6.2906    7.8599    9.4297
    2.0000    1.7035    3.2139    4.7614    6.3202    7.8837    9.4496
    3.0000    1.8644    3.3080    4.8261    6.3692    7.9231    9.4825

La solución R(r) es

R mn (r)= Y m ( ω mn c b ) J m ( ω mn c r ) J m ( ω mn c b ) Y m ( ω c r )

Comprobamos que cumple R(a)=0 y R(b)=0.

La solución completa es similar a la membrana circular, excepto en la parte radial, que ya hemos determinado

Ψ m,n (r,θ,t)={ Y m ( ω mn c b ) J m ( ω mn c r ) J m ( ω mn c b ) Y m ( ω c r ) }· ( C n cos(mθ)+ D n sin(mθ))·( A n cos( ω m,n t)+ B n sin( ω m,n t)) m=0,1,2,...n=1,2,3...

Representamos el modo de vibración de la membrana, para m=1 y para n=3, la raíz de la ecuación transcendente es ω13/c=4.7222 (segunda fila, cuarta columna de la tabla)

Ψ m,n (r,θ,t)={ Y m ( ω mn c b ) J m ( ω mn c r ) J m ( ω mn c b ) Y m ( ω c r ) }·cos(mθ)·cos( ω m,n t)

alfa=4.722232462741136;  
a=4; %radios del anillo circular
b=2;
r=b:0.01:a;
ang=0:pi/50:2*pi;
[r,ang]=meshgrid(r,ang);
x=r.*cos(ang);
y=r.*sin(ang);
fichero = 'anillo_c.gif';
P=1;
m=1; %modo de vibración
n=3;
for t=0:P/20:P
    z = (bessely(m,alfa*b)*besselj(m,alfa*r)-besselj(m,alfa*b)*bessely(m,alfa*r))
.*cos(m*ang)*cos(2*pi*t/P);
    surf(x,y,z)
    xlabel('x')
    ylabel('y')
    zlabel('z')
    title('Vibración de anillo circular')
     
    frame=getframe;
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    if t==0
        imwrite(imind,cm,fichero,'gif','DelayTime',0,'loopcount',inf);
    else
        imwrite(imind,cm,fichero,'gif','DelayTime',0,'writemode','append');
    end
end