Ecuación de vigas: Modelo de Euler-Bernoulli (13A)

De MateWiki
Revisión del 21:59 13 may 2015 de Isabel Roselló (Discusión | contribuciones) (→‎Flexión de una viga apoyada en sus extremos dependiendo del tiempo. Problema dinámico.)

Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Ecuación de vigas: Modelo de Euler-Bernoulli (13A)
Asignatura Ecuaciones Diferenciales
Curso Curso 2014-15
Autores María Aguilera, Paula Martínez, Miguel Sánchez, Laura García, Isabel Roselló, Sarah Boufounas
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Introducción

El presente trabajo trata de estudiar, mediante un modelo matématico, el comportamiento de una viga de 10 metros de longitud y sección rectangular, al someterla a la acción de diferentes esfuerzos y variar las dimensiones de su sección transversal, así como se estudiará la respuesta de esta misma viga al cambiar sus condiciones de apoyo.

Sección transversal de la viga
centro

Supondremos que la viga ocupa el intervalo x∈[0,L] y denotamos el desplazamiento vertical de su eje por y(x). Para desplazamientos pequeños, y(x) satisface la ecuación de la curva elástica que proviene del equilibrio de momentos:

[math]y'''=\frac{M(x)}{E I(x)}[/math]

donde E es el módulo de elasticidad lineal (o módulo de Young) que depende de las propiedades elásticas del material, y que supondremos constante, mientras que I(x) es el momento de inercia de la sección transversal respecto al centro. La función M(x) se conoce como momento flector y representa el momento de las fuerzas aplicadas sobre la viga.


2 Viga biapoyada sometida a la acción de momentos flectores

En este apartado se estudiará una viga biapoyada de la que en cada caso conoceremos distintos datos.

centro


2.1 [math] \ E=5x10^4 ; a=0.6 ; b=0.3 ; M(x)=L/2−|x−L/2|\ [/math]

Sabiendo que [math] \ E=5x10^4 , a=0.6 , b=0.3 , M(x)=L/2−|x−L/2|\ [/math], plantaremos un problema de contorno, con objeto de conocer la deformada de la viga al serle aplicado este momento flector.Una vez calculada la deformada estudiaremos sus puntos y encontraremos el de mayor deflexión (mayor y). Para realizar estos cálculos empleamos el siguiente código de matlab. 

% datos
L=10;               % longitud de la viga
E=5e4;              % módulo de Young
a=0.6;              % alto sección 
b=0.3;              % ancho sección 
I=(1/12)*b*a^3;     % inercia
 
%discretización 
x0=0;
xN=L;
N=50;
h=(xN-x0)/N;
x=x0:h:xN;
xx=x(2:N);
% f(x)
y0=0;yL=0;
M=L/2-abs(xx-L/2);
f=(M/(E*I))';       
f(1)=f(1)-y0/(h^2);
f(N-1)=f(N-1)-yL/(h^2);
% matriz K
K=(1/h^2)*(-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));

%solución
y=K\f;
 
%valores del contorno
y=[y0;y;yL];            
 
%deflexión máxima
fmax=min(y)
 
% dibujamos
plot(x,y)


A continuación se muestra la gráfica que representa la ley de deflexiones de la viga

Ley de deflexiones de la viga

El valor de la deflexión máxima es -0.1544 que como podemos observar en la gráfica, se produce en el centro de vano.

2.2 [math] \ E=5x10^4 ; a∈[0,1;0,9] ; b=1−a ; V_[viga]=cte \ [/math]

En esta ocasión, nuestra viga tendrá canto y ancho variable a lo largo de su longitud. Sabiendo que el volumen de esta permanecerá constante con respecto al caso anterior [math] \ (a=0.6 ; b=0.3 ; M(x)= L/2−|x−L/2|)\ [/math],procederemos a calcular la ley de deformadas de esta nueva viga. Para realizar estos cálculos volveremos a utilizar matlab. 

% datos 
L=10;        % longitud viga
E=5e4;       % módulo de Young
y0=0;
yL=0;

% discretización        
x0=0;
xN=L;
N=50;
h=(xN-x0)/N;
x=x0:h:xN;
xx=x(2:N);
M=L/2-abs(xx-L/2);
% matriz K
K=(1/h^2)*(-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1));
n=0;
fmax=zeros(1,9);
for a=0.1:0.1:0.9      % alto sección
    n=n+1;
    b=1-a;              % ancho sección
    I=(1/12)*b*a^3;     % inercia
    % f(x)
    f=(M/(E*I))';       
    f(1)=f(1)-y0/(h^2);
    f(N-1)=f(N-1)-yL/(h^2);
 
    %solución
    y=K\f;
    %valores del contorno
    y=[y0;y;yL];           
    fmax(n)=min(y);
 
 
    % dibujamos
    figure(1)
    hold on
    plot(x,y)
end
 fmax(n)=min(y)
hold off
figure(2)
plot(0.1:0.1:0.9,fmax,'-g')


De esta forma obtenemos que la menor deflexión máxima, con un valor de 0.0973, se produce cuando la sección tiene un canto de 0.7 y un ancho de 0.3

A continuación se muestran las gráficas que representan las deflexiones sufridas por cada sección de la viga y el gráfico que representa las máximas deflexiones para cada canto en esta misma viga.

Ley de deflexiones en cada sección
centro

2.3 Viga con sección variable respecto de x

En este apartado supondremos que la viga tiene seccion cuadrada con el lado a(x) variable:

[math] a(x)=cos(c*(x-L/2)) +d [/math]

Con c y d como variables pero que siempre cumplan que el volumen total de la viga sea igual al de la viga de sección cuadrada y lado a=0,5 . Por tanto esta condición nos permite sacar la relación entre c y d.

[math] L*a^2=\int_{0}^{L} [cos(c(x-L/2))+d]^2 dx[/math]

Operando llegamos a la siguiente relación:

[math] d=\frac{\frac{-4}{c}*sen(\frac{c*L}{2}) ± \sqrt{\frac{16}{c^2}*sen^2(c*L/2)-\frac{2*L}{c}*sen(c*L)-4*L^2(1/2-a^2)}}{2*L} [/math]

Tomando solo el valor positivo y hallando el intervalo de c haciendo que el discriminante de la raiz sea mayor que cero, planteamos el siguiente código. 

% Partición espacial
L=10;               % longitud de la viga
x0=0;xN=L;
N=50;dx=(xN-x0)/N;
x=x0:dx:xN;
xi=(x0+dx):dx:(xN-dx);

% Matriz K
KK=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
K=(1/dx^2)*KK;

% Datos 
E=5E4;              
y0=0;yL=0;
M=L/2-abs(xi-L/2);

n=0;
c0=-0.4354;cf=0.4354;
dc=0.02;
for c=c0:dc:cf         
    n=n+1;
    d=-4/c*sin(c*L/2)+sqrt((16*sin(5*c)^2)/(c*c)-(20*sin(10*c))/c-100);
    a=cos(c*(xi-L/2))+d; 
    I=(1/12).*a.^4;     
    
    % f(x)
    f=(M./(E*I))';       
    f(1)=f(1)-y0/(dx^2);
    f(N-1)=f(N-1)-yL/(dx^2);
    
    %solución
    y=K\f;
    
    y=[y0;y;yL];            
    fle_max(n)=-max(abs(y));
     
    figure(1)
    hold on
    plot(x,y)
    hold on
    end
hold off
figure(2)
plot(c0:dc:cf ,fle_max,'-ok')


izquierda
centro
La primera gráfica representa los distintos perfiles de las vigas para cada parametro de c y d.
La segunda las flechas según el valor de c.
Gracias a la ecuación obtenemos el valor óptimo y pésimo. La viga con menor deflexión (en rojo en la figura) se produce cuando [math] c=-0,0046 [/math] y [math] d=-9,9982 [/math] dando una flecha máxima en el centro de vano de [math] -1,5264·10^-6 [/math] m. La viga con mayor deflexión (en azul en la figura) se produce cuando [math] c=-0,3745 [/math] y [math] d=-4,3407 [/math] y una flecha máxima en el centro de vano de [math] -6,85·10^-5 [/math] m.
centro

3 Viga biempotrada sometida a la acción de una carga

centro

En el presente apartado se realizará un estudio homólogo al del apartado anterior con la diferencia de que en este caso en lugar de estar los extremos apoyados, estarán encastrados. La sección transversal tendrá el mismo valor para el canto y el ancho, es decir, nos encontramos ante una sección cuadrada con a=b=0.5. Con respecto a la acción que actúa sobre la viga nos encontramos ante otra diferencia con los casos anteriores: en esta ocasión en lugar de actuar momentos flectores actúa una carga de valor [math] \ w(x)= L/2−|x−L/2| \ [/math]. Para resolver este problema deberemos resolver la siguiente ecuación:

[math]y''''=\frac{-w(x)}{E I(x)}[/math]

Esta ecuación resulta ser de orden cuatro por lo que necesitaremos cuatro condiciones de contorno, que sacaremos de las condiciones de sustentación de la viga (empotramientos en los dos extremos). Las condiciones de contorno que resultan son las siguientes: \[\left\{\begin{matrix} y(0)=0 \ & \\ y'(0)=0\ \end{matrix}\right.\] \[\left\{\begin{matrix} y(L)=0\ & \\ y'(L)=0\ \end{matrix}\right.\] Para resolver el problema de contorno que tenemos planteado será necesario que planteemos el método de diferencias finitas, nos ayudará a calcular la solución mediante una aproximación. Utilizaremos las siguientes aproximaciones:

[math] y''''(x_j)=\frac{y(x_j-2)-4y(x_j-1)+6y(x_j)-4y(x_j+1)+y(x_j+2)}{h^4}+\theta(h^2)[/math]
[math] y'(x_0)=\frac{y(x_1)-y(x_-1)}{2h}+\theta(h^2)[/math]
[math] y'(x_N)=\frac{y(x_N+1)-y(x_N-1)}{2h}+\theta(h^2)[/math]

Para resolver el problema anterior y calcular el valor de mayor deflexión y el punto en el que se alcanza desarrollamos el siguiente programa de matlab: 

% datos generales
L=10;               % longitud de la viga
E=5e4;              % módulo de Young
a=0.5;              % altura sección 
b=1-a;              % ancho sección 
I=(1/12)*b*a^3;     % inercia
 

%discretización 
x0=0;
xN=L;
N=10;
h=(xN-x0)/N;
x=x0:h:xN;
xx=x(2:N); 


% f(x)
w=L/2-abs(xx-L/2);
f=-(w/(E*I))';       % vector columna
 
 
% matriz K
K=(1/h^4)*(6*diag(ones(1,N-1))-4*diag(ones(1,N-2),1)-4*diag(ones(1,N-2),-1)+diag(ones(1,N-3),2)+diag(ones(1,N-3),-2));
K(1,1)=7;
K(N-1,N-1)=7;

%solución
y=K\f;
 
y=[0;y;0];            % añadimos los valores del contorno
fmax=min(y)
 
W=L/2-abs(x-L/2);

% dibujamos
plot(x,W,x,y)


De esta manera obtenemos que la deflexión máxima se produce en la mitad del vano, es decir, en L/2=5, y con un valor de 0.3806

Representación de la carga (en rojo) y de la deflexión de la viga (en azul)

4 Viga cuya deformada dependede del tiempo

5 Flexión de una viga apoyada en sus extremos dependiendo del tiempo. Problema dinámico.

Consideramos ahora el problema dinámico. En este caso la deflexión de la viga depende también del tiempo y(x,t) quedando la ecuación:

[math]ρy_tt + EIy_xxxx = -w(x,t) [/math]

Siendo "p" la densidad de la viga que mantendremos constante y de valor la unidad. De esta forma nuestra ecuación final será:

[math] y_tt + EIy_xxxx = -w(x,t) [/math]

Al tratarse de un problema dependiente de 2 variables necesitamos dos condiciones iniciales [math] y(x,0)=y_0(x) \ ; \ y_t(x,0)=y_(x) [/math]

De esta forma planteamos el sistema de ecuaciones completo y su aproximación de diferencias finitas considerando el método del trapecio para resolver el problema en tiempo.

\[\left\{\begin{matrix}\ y=\ ky + b(t)\ , & \\ y(0)=y_0 \ , \\ y'(0)=z_0\ & \end{matrix}\right.\]

Tras realizar las derivaciones respecto al tiempo y respecto a x y despejando los valores conocidos nuestro sistema queda definido por:

[math] vec{y''}=k*vec{y} +vec{b(1)} [/math]

donde si tenemos en cuenta las condiciones iniciales resultaria: \[\left\{\begin{matrix}\vec{y}=k*\vec{y} +\vec{b(1)} \ , & \\\vec{y(0)}=\vec{y_0} \ , \\ \vec{y'(0)}=\vec{z_0}\ & \end{matrix}\right.\] Realizando el cambio de variable:

[math] \vec{y'}=\vec{z} [/math]
[math] \vec{z'}=k*\vec{y} +\vec{b} [/math]

Nuestro sistema final sobre el cual aplicaremos el método del trapecio es:

\[\left\{\begin{matrix} \vec{W'}= K*\vec{W} + \vec{B} \ , & \\ \vec{W(0)}=\vec{W_0} \ & \end{matrix}\right.\]

El código que a continuación mostramos nos permite calcular el sistema anteriomente mostrado siendo:

[math] y(0,t)=sin(\frac{16pix}{L}) [/math]
[math] y_t(0,t)=0 [/math]

Además nos permitirá obtener el comportamiento de y(x,t) en el punto x=0.7 en un intervalo de t [0,5]: 

%%% (ro)ytt+EIyxxxx=-w(x)
%%% ytt+(EI/ro)yxxxx=-w(x)/ro=f(x)
%%% y(0,t)=y'(0,t)=y(L,t)=y'(L,t)=0
%%% y(x,0)=(1/3)*sin(8*pi*x/L); yt(x,0)=0

clear all
close all

% datos generales
L=10;               % longitud viga
E=5E4;              % módulo de Young
ro=1;               % densidad
a=0.5;              % altura sección rectangular
b=1-a;              % anchura sección rectangular
I=(1/12)*b*a^3;     % momento de inercia

% partición espacial
x0=0;xN=L;
N=100;dx=(xN-x0)/N;
x=x0:dx:xN;
xi=(x0+dx):dx:(xN-dx);

% condiciones iniciales
y0=((1/3)*sin(8*pi*xi/L))';      % y(x,0) // column vector
z0=zeros(1,length(xi))';        % yt(x,0) // column vector
W0=[y0;z0];

% f(x)=-w(x)/ro
w=1*(L/2-abs(xi-L/2));
f=-(w/ro)';         % vector columna

% matriz K
KK=6*diag(ones(1,N-1))-4*diag(ones(1,N-2),1)-4*diag(ones(1,N-2),-1)+diag(ones(1,N-3),2)+diag(ones(1,N-3),-2);
KK(1,1)=7;KK(N-1,N-1)=7;
K=-(E*I/ro)*(1/dx^4)*KK;  

% matriz Q=[O,I;K,O]
I=eye(N-1);O=zeros(N-1);
Q=[O,I;K,O];

% partición temporal
tM=5;M=100;dt=tM/M;t=0:dt:tM;

% theta=0(Euler),0.5(Crank-Nicholson),1(Euler implícito)
theta=1/2;  % Crank-Nicholson

% matriz P=inv(L)*R
I2=eye(2*(N-1)); 
L=I2-dt*theta*Q;
R=I2+dt*(1-theta)*Q;
% P=L\R;

% término B
B=[zeros(length(f),1);f];
% B=[zeros(length(f),1);zeros(length(f),1)];

% programa
W=zeros(2*(N-1),M+1);   
W(:,1)=W0;

for n=1:M
    W(:,n+1)=L\(R*W(:,n)+dt*B);  
end

Y=W(1:N-1,1:M+1);
Y0=zeros(1,M+1);
YN=zeros(1,M+1);
Y=[Y0;Y;YN];           

figure(100)     % gráfico y(x,t)
[X,T]=meshgrid(x,t);
figure(100)
surf(X,T,Y') 

figure(200)     % curva y(0.5,t)
xc=0.5;
plot(t,Y(round(xc/dx+1),:))
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_de_vigas:_Modelo_de_Euler-Bernoulli_(13A)&oldid=29827»