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

Revisión actual del 23:58 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

Se va a resolver a continuación el problema para k=1. Las gráficas para el resto de valores de k requeridos, se obtendrán cambiando ese término en el programa. Los valores de k requeridos son, además de k=1, k=3.5,k=10 y k=20.

% 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]

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

Al ejecutar el programa expuesto inmediatamente encima, se produce un error en el comando "trapz" de la línea 73, el cual no hemos sabido subsanar, por lo que la gráfica resultado de la aplicación numérica del método de la energía no se la podemos ofrecer. En cualquier caso, a pesar de este inconveniente, podemos indicar que la energía es cero, tal como se demuestra a continuación de forma analítica

$\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} \qquad du=u_{xx}dx\\ dv=u_{xt}dx \qquad 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.

Se puede observar de la ecuación actual que la constante de amortiguamiento que produce el medio ira reduciendo el movimiento hasta finalmente su detención.
Nosotros trabajaremos entonces para [math]a=0,1,4,10,100[/math] .En el caso de que se obtuviera una gráfica comparativa de las diferentes constantes de amortiguamiento se podría observar con claridad que a medida que la constante va aumentando, el movimiento se va reduciendo hasta finalmente detenerse.: [math] \begin{cases} u_{tt}-u_{xx}+au_{t}=0; x∈[0,10]; t∈[0,60]\\ u(0,t)=0\\ u(10,t)=0\\ u(x,0)=0\\ u_{t}(x,0)=0\\ \end{cases}\\ [/math]

3.2 Cable sujeto a vibraciones periódicas en un extremo

En este caso suponemos que el extremo derecho del cable esta sujeto a una estructura que sufre vibraciones periódicas con frecuencia [math]F_0[/math]. El cable parte del reposo y estudiamos su comportamiento para [math]t∈[0,60][/math] El sistema que tenemos que estudiar en este caso será: [math] \begin{cases} u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,60]\\ \begin{cases} u(0,t)=0\\ u(10,t)=sin(2*pi*F_0*t)\\ \end{cases}\\ \begin{cases} u(x,0)=0\\ u_{t}(x,0)=0\\ \end{cases}\\ \end{cases}\\ [/math] Tomaremos los siguientes valores para [math]F_0[/math]:
[math]F_0=\frac{1}{L}+0.01[/math]
[math]F_0=\frac{1}{L}-0.01[/math]
[math]F_0=\frac{1}{L}[/math]
[math]L=10[/math]

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, y con las condiciones de frontera modificadas calculamos el comportamiento de la energía para [math]b=2[/math] y [math]b=-2[/math]:

[math] \begin{cases} u_{tt}-u_{xx}=0; x∈[0,10]; t∈[0,40]\\ u(0,t)=0\\ u_x(10,t)=bu(x,t)\\ u(x,0)=\begin{cases} \frac {2x}{5} & \text{ si } 0\ltx\lt5 \\ \frac {2(10-x)}{5} & \text{ si 5\ltx\lt10 } \end{cases}\\ u_{t}(x,0)=0\\ \end{cases}\\ [/math]