Cable de una estructura civil. (Grupo 16B)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Cable de una estructura civil. Grupo 16 |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2013-14 |
| Autores | Javier Díez Olaya 121 Javier Lozano Aragoneses 248 Enrique Martínez Mur 271 Begoña Bigeriego Alvarez 637 |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
1 Modelización de problema
El problema a estudio es una modelización de la ecuación de cuerda vibrante, que se presenta de la siguiente forma:
[math]
\rho u_{tt}-Zu_{xx}=f(x,y); x∈[0,L]; t\gt0;
[/math]
[math]\rho=\rho(x,t)[/math] = densidad lineal de la cuerda.
[math]Z=Z(x,t)[/math] = tracción o tensión.
[math]c=\sqrt{\frac{Z}{\rho}}[/math] = celeridad.
Está sometido a unas condiciones de contorno, CC, que gracias a nuestras condiciones del ejercicio sabemos que son unas condiciones Dirichlet puesto que tenemos nuestros extremos a cota conocida::
[math]
u(0,t)=g_{1}(t)\\
u(L,t)=g_{2}(t)
[/math]
[math]g_{1}(t)[/math] y [math]g_{2}(t)[/math] = Funciones de contorno de los extremos de la cuerda.
A su vez la cuerda vibrante funciona bajo unas condiciones iniciales, CI, que nos indicaran el perfil inicial de la cuerda así como la velocidad inicial que tienen sus puntos::
[math]
u(x,0)=A(x)\\
u_{t}(x,0)=B(t)\\
[/math]
[math]A(x)[/math] = Función que describe el perfil inicial de la cuerda situado a cota inicial [math]x=x_{0}[/math].
[math]B(t)[/math] = Velocidad a la que están sometidos todos los puntos de la cuerda en el instante inicial del movimiento[math]t=t_{0}[/math].
En nuestro caso, en el que consideramos un cable de una estructura civil de longitud [math]L = 10m[/math] sujeto por ambos extremos cuya sección es mucho menor con respecto de la longitud del mismo y en que sus vibraciones se pueden modelizar mediante la ecuación de ondas tendremos el siguiente probelma:
Vamos a modelizar una ecuación de ondas en la que consideraremos una cuerda vibrante sujeta por sus extremos::
[math]
\begin{cases}
u_{tt}-u_{xx}=0; x∈[0,10]; t\gt0\\
\begin{cases}
u(0,t)=0\\
u(10,t)=0\\
\end{cases}\\
\begin{cases}
u(x,0)=0\\
u_{t}(x,0)=0\\
\end{cases}\\
\end{cases}\\
[/math]
Interpretando en términos de cuerda vibrante tendremos:
[math]EDP[/math]: Cuerda homogénea de densidad lineal [math]\rho=1[/math] y tracción [math]Z=1[/math], por tanto tenemos también una celeridad de [math]c=1[/math], no está sometido a fuerzas por unidad de longitud en sus puntos interiores ya que [math]f(x,t)=0[/math] y la cuerda ocupa un intervalo de [math][0,10]m[/math].
[math]CC[/math]: Ambos extremos al estar empotrados están a cota [math]x=0[/math].
[math]CI[/math]: El perfil inicial de la cuerda al ser una sección pequeña con respecto de su longitud es de nula, [math]u(x,0)=0[/math] así como la velocidad inicial a la que están sometidos todos los puntos de la cuerda, [math]u_{t}(x,0)=0[/math].
2 Desplazamiento vertical del cable
2.1 Método del Trapecio
Sujetamos el cable desde el centro y lo desplazamos 2 m en la dirección perpendicular. Al soltarlo, este empieza a vibrar. Aproximar [math]u(x; t)[/math] por el método de diferencias finitas con [math]Δx = 0.1[/math], y usar el método del trapecio tomando [math]Δx = Δt[/math]. Dibujar la solución en tiempo [math]t ɛ [0,40][/math]
% Datos del problema
L=10;
T=40;
% Discretización espacial
dx=0.1;
N=L/dx;
% Vector de puntos espaciales
x=0:dx:L;
% Vector de espacio en los nodos interiores
xint=dx:dx:L-dx;
% Discretización temporal
dt=dx;
% Vector de tiempos
t=0:dt:T;
% Aproximación de -u_xx
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=(1/dx^2)*K;
F=zeros(N-1,1);
% Datos iniciales
u0=(2*xint)/5.*(xint<=5)+(2*(10-xint)/5).*(xint>5);
v0=zeros(N-1,1);
% Aproximación en tiempo
% Matriz M
M=[zeros(N-1), eye(N-1); -K, zeros(N-1)];
% Dato inicial
W0=[u0,v0']';
%Método del trapecio
WW=W0;
U=zeros(length(t),length(x));
% Definimos la matriz sol con la u para pintarla
sol=zeros(length(t),2*N);
sol(1,:)=[0,W0',0];
U(1,:)=[0,u0,0];
% Iteraciones W^j->W^j+1
for j=1:length(t)-1
WW=(eye(2*N-2)-(dt/2)*M)\((eye(2*N-2)+(dt/2)*M)*WW);
sol(j+1,:)=[0,WW',0];
U(j+1,:)=[0,WW(1:N-1)',0];
end
% Dibujamos la solución
figure(1)
[xx,tt]=meshgrid(x,t);
surf(xx,tt,U)
2.2 Método de Euler Explícito
Con las mismas condiciones que mediante el método del trapecio, realizamos la modelización a de nuestro sistema mediante Euler Explicito:
%Metodo Euler explicito
WW=W0;
U=zeros(length(t),length(x));
U(1,:)=[0, [W0(1:N-1)]',0];
for n=1:length(t)-1
WW=WW+dt*M*WW;
U(n+1,:)=[0, [WW(1:N-1)]',0];
end
[xx,tt]=meshgrid(x,t);
% Dibujamos la solución
figure(2)
surf(xx,tt,U)
2.3 Método de Euler Modificado
De la misma forma se realizará para el Euler modificado:
%Metodo de Euler Modificado
WW=W0;
U=zeros(length(t),length(x));
U(1,:)=[0, [W0(1:N-1)]',0];
for n=1:length(t)-1
k1=M*WW;
WW=(eye(2*N-2)+dt*M/2)*WW+(dt/2)*k1+(dt^2)*M*k1/2;
U(n+1,:)=[0, [WW(1:N-1)]',0];
end
[xx,tt]=meshgrid(x,t);
% Dibujamos la solución
figure(3)
surf(xx,tt,U)
2.4 Comparación de Métodos
La comparación de los Métodos del Trapecio, Euler explicito y Euler modificado dejan ver.................................................
2.5 Método de Fourier
% Datos del problema
L=10;
T=40;
k=1;
% Discretización espacial
dx=0.1;
N=L/dx;
% Vector de puntos espaciales
x=0:dx:L;
% Discretización temporal
dt=dx;
% Vector de tiempos
t=0:dt:T;
sol=zeros(length(t),N+1);
% Posición inicial
u0=(2*x)/5.*(x<=5)+(2*(10-x)/5).*(x>5);
v0=zeros(N-1,1);
% Fourier
for i=1:k
p=sin((i*pi/L)*x);
c=trapz(x,(u0.*p))/trapz(x,p.*p);
fp=c.*cos(i*pi*t/L);
sol=sol+[fp]'*p;
end
[xx,tt]=meshgrid(x,t);
% Dibujamos la solución
figure(4)
surf(xx,tt,sol)
3 Energía del cable
El método de las energías , también llamado de los multiplicadores ; se utiliza en este tipo de funciones para asegurar la unicidad de solución. Ésto , se demuestra partiendo de que el problema homogéneo asociado sólo admita la solución trivial, es decir, la nula; lo cual demostraremos , primero numéricamente y analíticamente a posteriori. La fórmula de la energía del cable:: [math]E(t)=\int^{L}_{0} (u^{2}_{t}(x,t)+(u^{2}_{x}(x,t)) dx[/math] Obtenemos las gráficas en las que nos apoyamos para determinar el valor de la energía, si ésta se conserva en el tiempo,..
% Calculamos la energía
ut=[WW(N-1:length(WW));0];
ux=zeros(length(t),length(x));
ux(1,:)=(2/5).*(x<=5)+(-2/5).*(x>5);
for q=2:length(t)-1
for J=1:length(x)-1
ux(q,J)=(U(q,J+1)-U(q,J))/dx;
ux(q+1,J)=(U(q,J)-U(q+1,J))/dx;
end
end
E=trapz(x,ut.^2)+trapz(x,ux.^2);
figure(2)
plot(E,t)Resolución analítica:
Tenemos que la energía total es igual a la suma de las energías potencial y cinética.
[math]E(t)=E_{c}(t)+E_{p}(t)[/math]
[math]E(t)=\int^{10}_{0} (u^{2}_{t}(x,t)+(u^{2}_{x}(x,t)) dx[/math]
[math]E'(t)=\int^{10}_{0} (2u_{t}(x,t)u_{tt}(x,t)+2(u_{x}(x,t)u_{xt}(x,t)) dx[/math]
Resolución de la integración por partes:
[math]
\begin{cases}
u=u_{x} ; du=u_{xx}dx\\
dv=u_{xt}dx ; v=u_{t}\\
\end{cases}\\
[/math]
[math]E'(t)=2\int^{10}_{0} u_{t}(x,t)u_{tt}(x,t)+2[u_{x}(x,t)u_{t}(x,t)|^{10}_{0}-\int^{10}_{0}u_{t}(x,t)u_{xx}(x,t) dx][/math]
[math]E'(t)=2\int^{10}_{0} u_{t}(x,t)u_{tt}(x,t)+2[u_{x}(10,t)u_{t}(10,t)-u_{x}(0,t)u_{t}(0,t)-\int^{10}_{0}u_{t}(x,t)u_{xx}(x,t) dx][/math]
Del enunciado del problema se tiene::
[math]u_{t}=\lim_{x \to{0}}{\frac{u(0,t+h)-u(0,t)}{h}}=0[/math]
:
[math]u_{t}=\lim_{x \to{0}}{\frac{u(10,t+h)-u(10,t)}{h}}=0[/math]
Resultando:
[math]E'(t)=2[\int^{10}_{0} u_{t}(x,t)u_{tt}(x,t)+u_{t}(x,t)u_{xx}(x,t) dx]=2\int^{10}_{0} u_{t}(x,t)\underbrace{[u_{tt}(x,t)-u_{xx}(x,t)]}_{0}dx=0[/math]
Lo que nos indica que [math]E(t)=cte [/math]
[math]E(t)=E(0)=\int^{10}_{0} \underbrace{[u_{t}^2(x,t)+u_{x}^2(x,t)]}_{0}dx[/math]
Entonces [math]u(x,t)=cte, cualquiera que sea t\gt=0 [/math]
Y como : [math]u(x,0)=0 [/math] \Longrightarrow [math]u(x,t)=0 para todo x en [0,10] y t\gt=0[/math]
3.1 Cable sumergido en medio viscoso
Al introducir el cable en un medio viscoso, se introduce una variable nueva hasta ahora desconocida que modifica nuestra ley de movimiento::
[math] u_{tt}-u_{xx}+au_t=0[/math]
[math]a[/math] = constante que depende del amortiguamiento que produce el medio.
Nosotros trabajaremos entonces para [math]a=0,1,4,10,100[/math] obteniendo una comparativa de todas en la siguiente gráfica:
3.2 Cable sujeto a vibraciones periódicas en un extremo
Ahora suponemos que nuestro cable esta sujeto en su extremo derecho a una estructura que sufre vibraciones periódicas con frecuencia [math]F_0[/math] bajo la función [math]f(t)=sin(2*pi*F_0*t)[/math]. Para lo que consideraremos[math]F_0=\frac{1}{L}+0.01[/math], que también el cable parte del reposo. Consideraremos un tiempo [math]t∈[0,60][/math].
En este caso varía la CI del extremo derecho, pasando el problema a ser:
[math] \begin{cases} u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,60]\\ u(0,t)=0\\ u(10,t)=sin(2*pi*F_0*t)\\ u(x,0)=0\\ u_{t}(x,0)=0\\ \end{cases}\\ [/math]
% Posición inicial
U=[0*xint];
V=[0*xint];
sol(1,:)=[0,U,0];
f=zeros(1,N-1);
for n=1:length(t)-1
Z=U+dt*V;
W=[V']-dt*K*[U'];
U=Z;
V=[W'];
sol(n+1,:)=[0,U,sin(2*(pi)*F0*t(n))];
end
% Dibujamos la solución
figure(3)
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
% Calculamos la energía
ut=[0,V];
ux=[U];
% Utilizamos el método de Euler implícito
for n=1:length(x)-1
ux=(eye(N-1)+dx*K)\(ux+dx.*f);
UX(n+1,:)=[0, ux'];
end
E=trapz(x,ut.^2)+trapz(x,UX.^2);
figure(4)
plot(E,t)
3.3 Cable sujeto a un aparato que envía una respuesta a la vibración que recibe
Nuestro extremo derecho del cable está ahora sujeto a un aparato que envía una respuesta a la vibración que recibe [math]u_x(L,t)=bu(L,t)[/math]. Bajo las mismas CI, condiciones iniciales, que se aplicaron al método del trapecio, calculamos el comportamiento de la energía para [math]b=2[/math] y [math]b=-2[/math]
