Diferencia entre revisiones de «Cable de una estructura civil. (Grupo 16B)»

Revisión del 23:06 19 may 2014

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:

Modelo de cable sujeto por ambos extremos expuesto a vibraciones.

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)

Método del Trapecio .

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:

% Ecuación de ondas con Euler explícito

% u_tt-u_xx=0, x en (0,L)
% u(0,t)=0
% u(L,t)=0
% u(x,0)=u0(x)
% u_t(x,0)=v0(x)

clear all

% 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=0*xint;


U(1,:)=[0,u0,0];
W=[u0]';
V=[v0]';

for n=1:length(t)-1
W=W+dt*V+F;
V=V-dt*K*W+F;
U(n+1,:)=[0,W',0];
end

[xx,tt]=meshgrid(x,t);

% Dibujamos la solución
figure(2)
surf(xx,tt,U)

Método del Euler Explicito .

2.3 Método de Euler Modificado

De la misma forma se realizará para el Euler modificado:

% Ecuación de ondas con Euler modificado

% u_tt-u_xx=0, x en (0,L)
% u(0,t)=0
% u(L,t)=0
% u(x,0)=u0(x)
% u_t(x,0)=v0(x)

clear all

% 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=0*xint;

U(1,:)=[0,u0,0];
W=[u0]';
V=[v0]';


for n=1:length(t)-1
a=K*W;
k1=[V a]';
k2=[V a]';
W=W+(dt/2)*(k1(1,:)'+k2(1,:)')+F;
a=K*W;
k1=[V a]';
k2=[V a]';
V=V-(dt/2)*(k1(2,:)'+k2(2,:)')+F;
U(n+1,:)=[0, W',0];
end

[xx,tt]=meshgrid(x,t);

% Dibujamos la solución
figure(3)
surf(xx,tt,U)

Método del Euler Modificado .

2.4 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)

Método de Fourier con 1 término.

Método de Fourier con 3 términos.

Método de Fourier con 5 términos.

Método de Fourier con 10 términos.

Método de Fourier con 20 términos.

2.5 Comparación de Métodos

En la comparativa de los diferentes métodos para la obtención de solución por diferencias finitas, a partir de las gráficas, podremos valorar:

• El método de Fourier lo consideraremos poco estable, ya que sólo va a ser fiable cuando consideremos un cuantioso número de términos. Vemos claramente en las gráficas como a partir de 10 términos , es cuando la gráfica se asemeja a los otros métodos.

• En cuanto a los otros 3 métodos, el menos preciso es el método de Euler; que crea un término a partir del valor del siguiente ( [math]U_{n}[/math] a partir de [math]U_{n+1}[/math]).

• Por la poca estabilidad de este método se introdujo el Euler modificado , también llamado Runge Kutta 2; que como su nombre indica va a tener orden 2 y va a ser fiable.

• Por último , el método del trapecio, es un método implícito de orden 2, que va a ser el que más se aproxime a la solución que buscamos. ( mejor estabilidad que los otros dos métodos).

Conclusión : El mejor método será el de Fourier siempre que podamos tomar un gran número de términos, sino, el del trapecio será el de mayor estabilidad.

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,..

% Aproximar la ecuacion de ondas
% u_tt-u_xx=0, x en (0,L)
% u(0,t)=0
% u(L,t)=0
% u(x,0)=u0(x)
% u_t(x,0)=v0(x)=0

clear all

% 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;
% Diferencias finitas
% 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);

% Discretización temporal
dt=dx;
% Vector de tiempos
t=0:dt:T;

% 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


% 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(5)
plot(E,t)

\underline{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]E(t)=cte[/math] cualquiera que sea [math]t\geq 0 [/math]
Y como :
[math]E(0)=0 [/math] => [math]E(t)=0[/math] para todo [math] x [/math] en [0,10] y [math] t\geq 0[/math]

Por tanto, la energía se mantiene constante siempre, aunque el paso de discretización cambie.

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]