Ecuación de vigas: Modelo de Euler-Bernoulli (14A)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Ecuación de vigas: Modelo de Euler-Bernoulli (14A) |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2015-16 |
| Autores | Adrián Piñeiro Novas, Antón Fernandez, Francisco Javier Vela Cobos, Ricardo Sanz Marin, Alexandro Cortijo Garcia |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Introducción
- 2 Viga biapoyada sometida a la acción de momentos flectores
- 3 Diseño de viga con menos deflexión
- 4 Viga de sección variable con deflexión mínima.
- 5 Viga biempotrada sometida a la acción de una carga
- 6 Problema dinámico en el que la flexión de una viga apoyada en sus extremos depende del tiempo.
1 Introducción
Estudiaremos mediante un modelo Euler-Bernouilli, 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.
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
Supongamos que E = 5 x 10^4 y que la viga se encuentra apoyada. En este caso, calcularemos las condiciones de contorno para y(x) aplicando el momento flector dado en el enunciado M(x) = L/2 - |x - L/2| Supondremos también que la sección transversal de la viga es un rectángulo de lados a = 0,6 m. y b = 0,3 m. Planteado el problema de contorno que verifica y(x), resolvemos por el método de diferencias finitas. Por último, determinaremos cuál es el punto de mayor deflexión de la viga, que será la menor ‘y’ dado que la deflexión es hacia abajo.
% 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
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.
3 Diseño de viga con menos deflexión
Tomamos diferentes valores de a en el intervalo [0.1;0.9] y de b = 1 - a de manera que el volumen de la viga permanezca constante. Calculamos la deflexión máxima para cada uno de los valores de a y deducimos qué diseño de viga tiene menor deflexión. El momento flector (M(x)) y el volumen de la viga permanecerán constantes con respecto al caso anterior (a=0.6 ; b=0.3).
% 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.
4 Viga de sección variable con deflexión mínima.
En este apartado supondremos que la viga tiene sección cuadrada con el lado a(x) variable: a(x)=cos(c*(x-L/2)) +d 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. Elegimos los valores de c y d que nos permitan tener la deflexión mínima. Para calcular el volumen de la viga integramos la sección a(x)^2 a lo largo de toda la longitud, lo que nos da la relación:
L·a^2=∫_0^L▒〖〖[cos(c(x-L/2))+d]〗^2 dx〗 Operamos para sacar la relación entre c y d. Para facilitar el calculo separamos la integral en tres sumandos: L·a^2=∫_0^L▒〖[1]dx+∫_0^L▒〖[2]dx+ 〗〗 ∫_0^L▒〖[3]dx 〗
∫_0^L▒〖[1]dx=∫_0^L▒〖〖cos〗^2 (〗 c·(x-L/2)) dx=L/2+sen(2c·(x-L/2))/4c]■(L@0)〗 = L/2+(sen(cL))/2c
∫_0^L▒〖[2]dx=2d∫_0^L▒〖cos(〗 c·(x-L/2)) dx=(2d·sen(c(x-L/2)))/c]■(L@0)〗 = (4d·sen(cL/2))/c
∫_0^L▒〖[3]dx=∫_0^L▒d^2 dx〗=Ld^2 L·a^2=L/2+(sen(cL))/2c+(4d·sen(cL/2))/c+Ld^2
Despejamos d, tomando el valor positivo: d=((-4)/c sen(cL/2)+√(16/c^2 〖sen〗^2 (cL/2)-4L^2 (1/2-a^2 )+2L/c sen(cL)))/2L
Calculamos los valores de c y d para los cuales el volumen permanece constante.
% 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')
- 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.
5 Viga biempotrada sometida a la acción de una carga
emos) 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 w(x) = L/2-|x-L/2| Para resolver este problema deberemos resolver la siguiente ecuación: y’’’’ = (-w(x))/(E.I) 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: y(0)=0, y'(0)=0, y(L)=0, y'(L)=0. 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: y'(x_j)=(y(x_(j-2) )-4.y(x_(j-1) )+6.y(x_j )-4.y(x_(j+1) )+y(x_(j+2) ))/h^4 +О(h^2 ) y’(x_0)=(y(x_1 )-y(x_(-1) ))/(2.h)+О(h^2) y’(x_N)=(y(x_(N+1) )-y(x_(N-1) ))/(2.h)+О(h^2) Definiendo los nodos artificiales x_(-1) 〖y x〗_(N+1) de forma adecuada.
% 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
6 Problema dinámico en el que la flexión de una viga apoyada en sus extremos depende del tiempo.
Aproximación de yxxxx=y(4(xi) por el método de diferencias finitas Desarrollando por Taylor y(x) en xi+2 y xi-2 se obtiene: xi- xi+2 =2h y(xi+2)=y(xi) + 2hy’(xi) + 〖(2h)〗^2/2!y’’(xi) + 〖(2h)〗^3/3!y’’’(xi) + 〖(2h)〗^4/4!y(4(xi) + O(h4) y(xi-2)=y(xi) - 2hy’(xi) + 〖(2h)〗^2/2!y’’(xi) - 〖(2h)〗^3/3!y’’’(xi) + 〖(2h)〗^4/4!y(4(xi) + O(h4) Sumando 1) y 2) y sustituyendo h^2=(y(x_(i+1) )-〖2y(x〗_i)+y(x_(i-1)))/(y(x_i)) Obtenemos: y^((4) (x_i)=(y(x_(i+2) )-4y(x_(i+1) )+6y(x_i )-4y(x_(i-1) )+y(x_(i-2)))/h^4 Es decir, se aproxima en cada nodo: ρy_tt+EIy_xxxx= -w(x,t)
% Datos
L=10; % longitud viga
E=5E4; % módulo de Young
a=0.5;
b=1-a;
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(16*pi*xi/L))'; % y(0,t)
z0=zeros(1,length(xi))'; % yt(0,t)
W0=[y0;z0];
% f(x)=-w(x)/ro
w=1*(L/2-abs(xi-L/2));
f=-(w)';
% 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)*(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=1/2; % theta=0.5(trapecio)
% Matriz P=inv(L)*R
I2=eye(2*(N-1));
L=I2-dt*theta*Q;
R=I2+dt*(1-theta)*Q;
% Término B
B=[zeros(length(f),1);f];
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];
% Grafico y(x,t)
figure(1)
[X,T]=meshgrid(x,t);
figure(1)
surf(X,T,Y')
figure(2)
xc=0.7;
plot(t,Y(round(xc/dx+1),:))
- La primera gráfica representa el comportamiento de y(x,t) entre los intervalos dados. La segunda representa el comportamiento frente al tiempo para x=0,7
