Estudio de la carga crítica de una columna

Trabajo realizado por estudiantes
Título	Estudio de la carga crítica de una columna. Grupo 14-B
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2013-14
Autores	Manuel Campayo Fernández, Alejandro Arturo Serrano Leo, David Suárez Cuesta, Arturo Rodríguez Gárate, Jorge Martín Sebastián, Joaquín Sánchez Molina
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Se nos plantea analizar cual es la carga crítica de una columna, para lo cual realizamos el siguiente estudio. Cogemos una columna homogénea y de sección transversal circular. La colocamos de forma vertical y aplicamos sobre ella una carga P. El eje de simetría de la columna, que se extiende desde x=0 (extremo superior de la viga) hasta x=L (extremo inferior), sufrirá con la carga P desplazamientos horizontales (eje Y) respecto a la posición de equilibrio. Para medir estos desplazamientos utilizaremos una función y(x) denominada deflexión, la cual (si los desplazamientos son pequeños) es posible aproximar mediante la resolución de la ecuación de la curva elástica que proviene del equilibrio de momentos y que mostramos a continuación:: [math] y''(x) = \frac{M(x)}{EI(x)} \left\{\begin{matrix} E=Módulo De Elasticidad Lineal. Lo Supondremos Constante. \\ I(x)=Momento De Inercia De La Sección Trasversal Restecto Al Centro \\ M(x)=Momento Flector \end{matrix}\right. [/math]

En este caso el momento depende de la deflexión, por lo que [math] M(x)=-PI(x) [/math]

1 Problema de contorno para y(x)

El problema físico planteado no se resuelve mediante una ecuación diferencial en derivadas parciales, sino con un problema de contorno, que en nuestro caso, aplicándole todas las condiciones antes explicadas queda de la siguiente manera:: [math] (PA) = \left\{\begin{matrix} E*I(x)*y''(x)=M(x) \\ y(0)=y(L)=0 \end{matrix}\right\}. [/math]

La condición y(0)=y(L)=0 viene determinada por el hecho de que la viga no sufre ningún tipo de desplazamiento horizontal en sus extremos.

Sustituyendo en lo anterior la condición de [math] M(x)=-PI(x) [/math] , el problema de contorno que buscamos queda:: [math] (PA) = \left\{\begin{matrix} E*I(x)*y''(x)+P*y(x)=0 \\ y(0)=y(L)=0 \end{matrix}\right\}. [/math]

2 Valores de P que hacen estable a la columna

Una vez tenemos el problema de contorno planteado, queremos obtener para qué valores de P la columna es estable. Para ello, suponemos que E=1, la sección de la columna circular es de radio constante R=1 m, su densidad es ρ=1[math]Kg/m^3[/math] y L=10 m.

Antes de nada, se debe aclarar que una columna es estable si la única solución posible de la elástica es la trivial [math] y(x)=0 [/math].

El problema (PA) es homogéneo, por tanto siempre verificará la solución trivial. Ahora bien, este problema es estable sólo si verifica esa y ninguna otra, por tanto, cuando haya infinitas soluciones, el problema será inestable. Esto se traduce en un pandeo de la columna y en el posterior derrumbamiento de la estructura. Por ende, sólo tendremos que determinar los valores de P que hacen que el problema de contorno (PA) tenga infinitas soluciones.

Sustituyendo los valores de L, E, R y densidad que nos dan en el enunciado, el problema de autovalores a determinar queda:: [math] (PA) = \left\{\begin{matrix} \frac{1}{4} π*y''(x)+P*y(x)=0 \\ y(0)=y(10)=0 \end{matrix}\right\}. [/math]

Resolviendo el problema de autovalores obtenemos los valores de la carga P que hacen que (PA) tenga infinitas soluciones, provocando la inestabilidad de la barra y su posterior pandeo. La solución del problema de autovalores es la siguiente:

Los autovalores son : [math] P_k = \frac{k^2*π^3}{400} [/math] y las autofunciones asociadas son \[\begin{array}{crl}

\\

Φ_k = Asin\; \left(\frac{kπx}{10}\right)\\

\\

\end{array}\]

con k=1,2,...,N y A número real.

Por tanto, los valores de P que hacen la columna estable son todos los P menores o iguales a 0 y aquellos que no verifican:: [math] P_k = \frac{k^2*π^3}{400} [/math] con k=1,2,...,N.

Así la columna sería estable.

3 Función de las cargas críticas mínimas para las cuales la columna deja de ser estable

En el apartado anterior, estudiábamos el problema para R=1. Sin embargo, en este caso existen varios valores de R a estudiar, concretamente los pertenecientes al intervalo [1,5]. Si resolvemos el problema de autovalores (PA) en función del radio (aportado por la expresión del momento de inercia del círculo), obtenemos la expresión de los autovalores en función del radio de la sección del pilar.

Los autovalores son : [math] P_k = \frac{k^2*π^3*R^4}{400} [/math] y las autofunciones asociadas son \[\begin{array}{crl}

\\

Φ_k = Asin\; \left(\frac{kπx}{10}\right)\\

\\

\end{array}\]

con k=1,2,...,N y A número real.

Introduciendo esa expresión en Matlab (código 1) y evaluándola en el intervalo [1,5], obtenemos los valores de P_crit. Para ello, en vez de tomar únicamente cinco valores de R, hemos tomado un salto más pequeño, en concreto h=0.1. Por tanto, el código utilizado es:

%Representamos los valores de la carga critica en funcion del radio de la sección.
clear; clf; clc;
r1=1;
r2=5;
h=0.01;
N=(r2-r1)/h;
R=r1:h:r2;
P=zeros(1,N+1)
for k=1:(N+1)
P(1,k)=(pi^3*R(1,k).^4)/400
end
plot(R,P,'b-')
xlabel('Radio (m)');
ylabel('Carga critica (kN)')

Introduciendo la función en Matlab, la gráfica que se obtiene es la siguiente:

Al haber tomado un valor tan pequeño del módulo de elasticidad, E=1, y resolviendo, los primeros valores de la gráfica son muy pequeños, pero todos distintos de cero.

Nota: la carga crítica (P_crit) se obtiene particularizando la expresión de los autovalores para k=1 y físicamente se corresponde con la carga que provoca el pandeo del pilar. El resto de los autovalores (k=2,3, …) tienen valores más elevados y su significado físico es más complejo, el cual tampoco es interesante determinar puesto que la sobrecarga que hay que introducir para que se produzca es mucho más elevada que en el caso del primer autovalor. Estudiamos el problema para el caso k=1.

3.1 Cálculo numérico

Ahora se tratará de resolver por diferencias finitas el problema de autovalores cuya solución analítica se calculó en los apartados anteriores. Se trata de un problema sencillo, con condiciones del tipo Dirichlet en la frontera y homogéneo.: [math] (PA) = \left\{\begin{matrix} \frac{1}{4} π*R^4*y''(x)+P(x)=0 \\ y(0)=y(10)=0 \end{matrix}\right\}. [/math]

En este caso, en nuestro problema desconocemos R (pertenece a un intervalo, no podemos introducirlo directamente), y desconocemos también P. Sabiendo que:: [math] y'' = \frac{y(x_n-1)-2y(x_n)+y(x_n+1)}{h^2}+σ(h^2) [/math] Por tanto la matriz \begin{equation*}

K

=

\begin{bmatrix}

2 & -1 & 0 & 0 & ... & 0 \\

-1 & 2 & -1 & 0 & ... & 0 \\

... & ... & ... & ... & 2 & -1

\end{bmatrix}

\end{equation*}

Para resolverlo numéricamente, empleamos el comando “eig”, al que introducimos la matriz cuadrada K. La respuesta del programa es una matriz A que contiene las cargas críticas correspondientes a cada R y a cada valor de k. La columna que nos interesa a nosotros es la primera, dado que es donde aparecen los valores de la carga crítica correspondiente a k=1.

El código de Matlab que hemos utilizado para esto es:

% Calculamos numéricamente los valores de P_cr. El primer autovalor es el que se corresponde con dicha carga (esto está almacenado en la primera columna de la matriz A)
clear, clf;
xa=0;
xb=10;
ya=0;
yb=0;
N=100;
h=(xb-xa)/N;
x=xa:h:xb;
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
R=linspace(1,5,N+1)
for i=1:length(R)
A(i,:)=eig(0.25*pi*R(i)^4.*K)
end p=A(:,1)
%Ahora ya podemos resolver el problema de contorno.
plot(x,p,'-')

Dándonos de este modo la gráfica:

4 Estudio de una columna con radio variable

4.2 Gráfica de P_crit para diferentes valores de R_o

A continuación vamos a ver cómo varía el valor de la Pcrit cuando cambiamos el valor del radio de la viga. Para ello, hemos ampliado el código de Matlab utilizado en el apartado anterior. Así conseguimos saber cómo varía la P_crit cuando la R₀ varía desde 0,5 hasta 1,5. Además, en apartados posteriores, querremos saber para qué valores de R₀ la P_crit es mínima y máxima, por lo que lo reflejamos también en el código. El código ampliado es el siguiente:

%Resolvemos ahora el problema suponiendo el abombamiento de la viga.
clear, clf;
 
xa=0;xb=10;
ya=0;yb=0;
N=10; h=(xb-xa)/N;
x=xa:h:xb;
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
%Definimos los parámetros a b y c
R0=linspace(0.5,1.5,N+1)
%Es un bucle doble anidado.
for j=1:length(R0)  
    a=(-333.33*R0(j)+sqrt(-22222.22*(R0(j))^2+133333.33))/6666.67
    b=-10*a
    c=R0(j)
    for k=1:length(x)
        R(j,k)=a*x(k)^2+b*x(k)+c
    end
end
for i=1:length(R)
    Rsm(i)=min(R(i,:))
end
%El radio de las columnas ahora será variable. Calcularemos los autovalores en el punto de la sección en el que el radio sea menor, así obtendremos el Pcr en toda la columna
 
for i=1:length(R)
    A(i,:)=eig(0.25*pi*Rsm(i).^4.*K) %Los autovalores en la sección más estrecha de cada R_0
end
p=A(:,1) %De aquí sacamos que la columna con menor P es la número 2 (con R_0=0.6m) y la de mayor la número 11 (con R_0=1.5m)
figure (1) %Dibujamos la columna número 2
plot(x,R(2,:),'-') %En este caso R_0 es 1, por lo que la sección es como la resuelta en los apartados anteriores (circular, sin cambio de radio)
grid on
title('Columna con menor carga crítica (admisible)')
xlabel('x (m)'); ylabel('Radio de la columna (m)')
figure (2)
plot(x,R(11,:),'-')
title('Columna con mayor carga crítica (admisible)')
xlabel('x (m)'); ylabel('Radio de la columna (m)')
grid on
figure (3)
plot(R0,p,'x')
title('Evolución de la carga crítica con el valor de R_0')
xlabel('R_0'); ylabel('P_cr')
grid on
figure (4) %Dibujamos la columna número 1
plot(x,R(1,:),'-') %Dibujamos la columna con R_0=0,5 para demostrar que su volumen es mayor que el de R_0=0,6
grid on
title('Columna con R_0=0,5 (admisible)')
xlabel('x (m)'); ylabel('Radio de la columna (m)')

La gráfica en la cual se ven los valores de P_crit desde R₀=0,5 hasta R₀=1,5 es la siguiente:

4.3 Valores de R₀ para los que P_crit es máximo y mínimo

Los valores exactos de P_crit son los siguientes:

Como podemos ver, el valor mínimo de P_crit no es el del menor radio (R₀=0,5) como podría decirnos la lógica, sino el segundo menor (R₀=0,6). Esto sucede porque si vemos la gráfica de la columna con radio R₀=0,5 llegamos a una conclusión:

Como podemos ver, el gráfico nos dice que hay una zona en la que el radio de la columna sería negativo. Aunque esto sea físicamente imposible, matemáticamente esto provocaría que el volumen del sólido de revolución de esta parábola alrededor de la línea de “Radio de la columna”=0 sería mayor que el del sólido de revolución de la parábola correspondiente a R₀=0,6. Dicha parábola es la siguiente:

Fijándonos en los valores exactos de la P_crit vemos que el máximo valor de la P_crit se da para R₀=1,5. La parábola correspondiente es la siguiente: