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\;respecto\;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, por existir ahí sendos apoyos fijos.

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 antes hemos expuesto y utilizando el momento de inercia de un círculo, ya que la columna es de sección transversal circular de radio R=1 m constante:: [math] I = \frac{π·R^4}{4} \Longrightarrow I = \frac{π}{4} [/math] 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. Para ello actuamos del siguiente modo (como herramienta de cálculo sustituimos [math] y(x)=y [/math]): \[\begin{array}{crl}

\\
\frac{1}{4} π·y+P*y=0 \Longrightarrow  \frac{1}{4} π·r^2+P·r=0 \Longrightarrow  r=\; \pm\;\sqrt{\frac{-4P}{π}}i \\
\\

\end{array}\] En este momento contemplamos tres casos:: [math] \left\{\begin{matrix} P\gt0 \\ P=0 \\ P\lt0 \end{matrix}\right. [/math] No existen autovalores que verifiquen P=0 y P<0, por lo que la columna es estable para toda carga P≤0. Pasamos, de esa manera, a estudiar el caso P>0:

La solución general de la ecuación diferencial homogénea será \[\begin{array}{crl}

\\

y\;= y(x)\; =\; A\; sin\; \left(2\sqrt{\frac{P}{π}} x\right) + \; B\; cos\; \left(2\sqrt{\frac{P}{π}} x\right) = 0\\

\\

\end{array}\]

Obligamos a que se satisfagan las condiciones iniciales de nuestro problema de contorno y obtenemos los siguientes resultados: \[\begin{array}{crl}

\\
y(0)=0 \Longrightarrow  A·0+B·1=0 \Longrightarrow  B=0 \\
\\

\end{array}\] \[\begin{array}{crl}

\\
y(10)=0 \Longrightarrow   0\; =\; A\; sin\; \left(20\sqrt{\frac{P}{π}} \right) \Longrightarrow  A\; sin\; \left(20\sqrt{\frac{P}{π}} \right) =0 \\
\\

\end{array}\] Así, despejando esa P en esa última igualdad, obtenemos que 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 m. 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) del mismo modo que lo hicimos en el anterior apartado, 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 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 función 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(R,p,'-')
xlabel('Radio (m)');
ylabel('Carga critica (kN)')

Dándonos de este modo la gráfica:

4 Estudio de una columna con radio variable a lo largo de la misma

Ahora, en vez de tomar un radio constante para toda la columna, el radio de cada sección irá variando en función de su distancia a la base. De esta forma, tendremos la expresión R(x), que será un polinomio de 2º grado simétrico respecto al centro de la columna. Con todo esto, tendremos $R(x)=ax^2+bx+c$, y las condiciones complementarias es que en los dos extremos, el radio de la sección es R₀. Además sabemos que la masa total de la columna es de 31 kg, que con ρ=1[math]Kg/m^3[/math], así que su volumen es de 31 [math]m^3[/math].

De esta forma tenemos que:

R(0)=R₀

R(10)=R₀

[math]\displaystyle\int_{0}^{10} πR(x)^2\,dx=31 \Longrightarrow \displaystyle\int_{0}^{10} π(ax^2+bx+c)^2\,dx=31[/math]

Calculando a, b y c tenemos que:: [math]c= R\ltsub\gt0\lt/sub\gt[/math], [math]b=-10a[/math] y [math] a = \frac{-333,3*R\ltsub\gt0\lt/sub\gt\;\pm\\;\sqrt{-22.222,22*(R\ltsub\gt0\lt/sub\gt)^2+133.333,33}}{6.666,67} [/math]

Sabemos que para el valor de R₀=1, R(x) es una recta, por lo que el término a es 0. Para ello, al resolver el término en a debemos tomar el signo positivo en el más menos. Una vez fijado a, tenemos b por la expresión anterior. Con estos datos ya tenemos R(x) en función de R₀ únicamente.

4.1 Gráfica de P_crit para diferentes valores de R₀

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)')
for k=1:length(p)
    g(k)=(p(k)-p(round(length(p)/2)))*100/p(round(length(p)/2)) %Calculamos así la ganancia (o pérdida) en porcentaje de Pcrit para cada valor de R0 respecto a la columna de seccion circular constante (R0=1)
end
figure (5) %Dibujamos el cambio en porcentaje de Pcrit para cada valor de R0 respecto a la columna de seccion circular constante (R0=1)
plot(R0,g,'.-')
grid on
title('Ganancia en porcentaje de Pcrit respecto a la columna circular de seccion constante')
xlabel('R0'); ylabel('Ganancia en orcentaje de Pcrit')
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:

siendo los valores exactos de P_crit 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 R₀=0,5 llegamos a una conclusión:

Radio de la sección a lo largo de la columna

4.2 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 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:

Carga mínima en función de la distancia a la base

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:

Carga máxima en función de la distancia a la base

4.3 Comparación de la ganancia en P_crit (en porcentaje) respecto a la columna de sección circular de radio constante

Por último, vamos a comparar las ganancias en Pcrit, en porcentaje, respecto al diseño en el que la columna tiene sección circular de radio constante, manteniendo la masa de la columna en 31 kg. Esto se da cuando la R0 es igual a 1. La fórmula utilizada para calcular esto es la siguiente: Pcrit(R0)-Pcrit(R0=1)*100 Pcrit(R0=1) Los valores que nos salen son los siguientes:

5 Método de Fourier para radio de la sección transversal constante

Para terminar nuestro estudio nos disponemos a resolver el problema de contorno con [math] M(x)=-PI(x)-x [/math] por el Método de Fourier, tomando de nuevo el radio de la sección transversal constante. Para hacer un estudio más detallado de este problema que nos planteamos, vamos a resolverlo tomando 3, 5, 10 y 40 términos de la serie.

Recordamos, en este punto, que el problema de contorno que teníamos en los primeros apartados de este trabajo era:: [math] (PA) = \left\{\begin{matrix} E·I(x)·y''(x)=M(x) \\ y(0)=y(L)=0 \end{matrix}\right. \Longrightarrow (PA) = \left\{\begin{matrix} \frac{1}{4} π·y''(x)+P·y(x)=0 \\ y(0)=y(10)=0 \end{matrix}\right. [/math]