Carga crítica de una columna, Trabajo 9, (Grupo 5)

Trabajo realizado por estudiantes
Título	Carga crítica de una columna. Grupo 5-A
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2015-16
Autores	Noemi Palomino Maria Lucia Cueva
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1. Introducción

En este tema se estudia el fenómeno de inestabilidad de piezas simples solicitadas a compresión centrada. Se trata por lo tanto de la primera aproximación a los fenómenos de pandeo de pilares, en la que se restringe el enfoque el caso de cortante puro. Así mismo, en este tema se tratará de forma general, la carga crítica a la que se somete una columna se sección circular constante de una pieza simple e ideal. Aunque la información puede ser extrapolada a los pilares con sección variable, su caso se analizará más adelante. Para la exposición del pandeo por composición pura de piezas simples se comenzará con los fundamentos teóricos del problema.

2. La pieza ideal

Se denomina pieza ideal a una abstracción matemática de los pilares reales, basada en las siguientes hipótesis:

• El material tiene una rigidez constante caracterizada por un módulo de deformación E.

• El material es infinitamente elástico, sin que se alcance nunca un límite de elasticidad, o tensión de agotamiento.

• La pieza no tiene imperfecciones de ninguna clase, ajustándose de manera perfecta a su directriz teórica.

• No existen tensiones residuales dentro de la sección

Con estas simplificaciones se pueden plantear de manera sencilla las ecuaciones que gobiernan el problema de inestabilidad, y su resolución matemática

3. Carga crítica

Vamos a analizar la posible inestabilidad de la columna biapoyada de la figura siguiente. Es evidente que la solución trivial de la pieza consiste en permanecer recta, sin deformación. Sin embargo, la posible aparición de una flecha y(x), hace que el axil (P) provoque un momento flector que a su vez puede incrementar la flecha. Para que la deformación sea posible, debe encontrarse una solución de equilibrio sobre la pieza deformada, entre los esfuerzos solicitantes y la respuesta de la sección. Esta solución de equilibrio es un análisis de segundo orden. El planteamiento del equilibrio en cualquier sección de la pieza deformada es inmediato:

Momento solicitante:

EI: Módulo de elasticidad lineal (o módulo de Young) que depende de las propiedades elásticas del material, y que supondremos constante

I(x): Momento de inercia de la sección transversal respecto del centro

3.1 Estudio de la pieza de equilibrio

Para que la pieza esté en equilibrio, se debe resolver la ecuación diferencial de segundo orden. Para resolverla, suponemos que los desplazamientos son pequeños y aproximamos y(x), resolviendo la ecuación de la curva elástica que proviene del equilibrio de momentos:

Se trata de una ecuación diferencial homogénea que admite solución general:

3.2 Problema de contorno de la pieza de equilibrio. Concepto de carga crítica

Para obtener la solución se deben particularizar las condiciones de contorno:

Lo que da lugar a:

Este sistema tiene una solución trivial que es A=0, con lo que la pieza está en equilibrio en la situación sin deformar, basado en la inexistencia de momentos flectores. Sin embargo, también admite una solución cuando la función trigonométrica se anula:

En este caso, el parámetro A queda indeterminado, pudiendo adoptar cualquier valor, este es un problema de contorno homogéneo, por ello tendremos que determinar los valores de P que hacen que el problema de contorno tenga infinitas soluciones. Esto supone que la pieza admite infinitas posiciones de equilibrio en situación deformada, cuando se alcanzan unos determinados valores del axil aplicado (carga crítica). Hemos conseguido resolver el problema de autovalores asociado, dados por la ecuación siguiente:

A partir de esta ecuación podemos observar que la deformada de la pieza vendrá dada por los distintos valores de cargas críticas, y que estos dependen a su vez de n (n=1,2,3…)

En la siguiente figura se ha hecho una representación de los diferentes axiles críticos. Se muestran los casos de los tres axiles críticos con los valores de n=1,2,3. La deformada que representa en la figura se le llama modo de pandeo.

Con ello se explica que para efectos prácticos, la carga crítica relevante es la del primer modo, ya que una vez alcanzada la flexión, la pieza seguirá deformándose ante cualquier perturbación transversal. Se denomina carga crítica de Euler:

4. Estabilidad de P

Se dice que una columna es estable cuando la único solución posible es la trivial y(x)=0.

4.1 Ejemplo de estabilidad de una columna con sección circular

Sea

Para estudiar la estabilidad de la columna, necesitamos saber cuáles son los valores de P, para ello sustituimos en (1), los datos enunciados:

Podemos deducir que habrá estabilidad cuando la carga no tenga un valor superior a 0,0470 KN

4.1.1 Gráfica de estabilidad en función del radio

Para comprobar que la carga crítica que se ha calculado es correcta, vamos a realizar un programa que asocia a cada radio el valor más pequeño de la carga crítica para el cual la columna deja de ser estable.

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)=1.2*(pi^3*R(1,k).^4)/(4*(9^2));
 end
plot(R,P,'b-')
xlabel('Radio (m)');
ylabel('Carga critica (kN)')

La carga crítica, se obtiene particularizando con n=1, como hemos comentado anteriormente, ya que es la carga que provoca el pandeo del pilar. En el apartado anterior, solo hay un radio, sin embargo para este ejemplo se ha utilizado un intervalo de radios (de 1 a 5)

En la gráfica se observa que un radio mayor permite un mayor valor para la carga crítica.

4.2 Ejemplo de estabilidad de una columna con sección cuadrada

Otro ejemplo particular, es el estudio de una columna con una sección cuadrada. Este ejemplo se resuelve de la misma forma que el apartado 4.1 solo que el momento de inercia cambia:

Vamos a suponer que la base (B), tiene la misma longitud que el radio de la circunferencia, por lo que nuestro nuevo valor para la carga crítica será:

Podemos comprobar que en una sección cuadrada la carga crítica disminuye, esto se debe al momento de inercia

5.Estabilidad por diferencias finitas

Sea nuestra ecuación diferencial

Escribimos una aproximación por diferencias finitas del apartado 4.1:

El método de Diferencias Finitas es un método de carácter general que permite la resolución aproximada de ecuaciones diferencias en derivadas parciales definidas en recintos finitos.

En este caso vamos a calcular una aproximación de la carga crítica. Vamos a calcular el problema matricialmente de forma que podamos representarlo numéricamente

Aproximando:

De esta forma, la ecuación diferencial nos quedará de la forma:

Y pasándolo a forma matricial:

Por lo tanto, la matriz de coeficientes nos queda de la forma:

Se trata por lo tanto, un problema en el que tendremos que aproximar nuestro autovalor (que es la carga crítica) y encontrar una solución de éste:

Numéricamente, se emplea el comando “eig”, que nos permite calcular los valores propios de la matriz A:

%Calcular numéricamente por diferencias finitas 
%los valores de P(cr). Para ellos se realizará
%un problema de autovalores
 
%Parte 1. Condiciones de Contorno
xa=0;
xb=9;
ya=0;
yb=0;
%Parte 2. Paso y tamaño de paso 
N=100;
h=(xb-xa)/N;
x=xa:h:xb;
%Parte 3. 
%%Matriz K
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
%%Vector de longitudes (Radio de 1 a 5)
R=linspace(1,5,N+1);
 
%Parte 4. Bucle
for i=1:length(R)
A(i,:)=eig(1.2.*0.25*pi*R(i)^4.*K);
%La matriz A es EI/M*K
end
p=A(:,1);
 
 
%Dibujo del problema de contorno
plot(R,p,'r')
xlabel('Radio (m)')
ylabel('Carga crítica (kN)')

clear all, clc
 
%Caso 1. Calculado a partir del autovalor Pr
 
L=input('introduce una longitud: ');
r=0:0.1:5;
%función del pcr 
Pcr=((pi^3).*(r.^4).*1.2)./(4.*(L^2))
 
%Caso calculado por aproximación de diferencias finitas
%Parte 1. Condiciones de Contorno
xa=0;
xb=9;
ya=0;
yb=0;
%Parte 2. Paso y tamaño de paso 
N=100;
h=(xb-xa)/N;
x=xa:h:xb;
%Parte 3. 
%%Matriz K
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
%%Vector de longitudes (Radio de 1 a 5)
R=linspace(1,5,N+1);
 
%Parte 4. Bucle
for i=1:length(R)
A(i,:)=eig(1.2.*0.25*pi*R(i)^4.*K);
%La matriz A es EI/M*K
end
p=A(:,1);
 
 
%Dibujo de la comparación de los 2 programas
plot(R,p,'r',r,Pcr,'*')
xlabel('Radio (m)')
ylabel('Carga crítica (kN)')

En la figura 4 vemos que las gráficas coinciden. Esto es porque la aproximación numérica difieren a partir del tercer decimal, por lo que podemos considerar ese error despreciable.

Es decir, calculando la carga crítica con n=1 por autovalores o aproximándolo por diferencias finitas, la carga crítica difiere en 〖10〗^(-3) m, considerando así los valores calculados correctos.

5.1 Comparación de la Carga crítica en sección circular y rectangular

Uno de los estudios más interesantes es la comparación de la carga crítica en función de su sección

%Parte 1. Condiciones de Contorno
xa=0;
xb=9;
ya=0;
yb=0;
%Parte 2. Paso y tamaño de paso 
N=100;
h=(xb-xa)/N;
x=xa:h:xb;
%Parte 3. 
%%Matriz K
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
%%Vector de longitudes (Radio de 1 a 5)
R=linspace(1,5,N+1);

%Parte 4. Bucle
for i=1:length(R)
A(i,:)=eig(1.2.*(0.25)*pi*R(i)^4.*K);
%La matriz A es EI/M*K
end
p=A(:,1);

for i=1:length(R)
A(i,:)=eig(1.2.*(1/12)*pi*R(i)^4.*K);
%La matriz A es EI/M*K
end
m=A(:,1);
%Dibujo del problema de contorno
plot(R,p,'r',R,m,'b')
xlabel('Radio (m)')
ylabel('Carga crítica (kN)')

En la figura 5 Observamos que la sección rectangular admite menos carga que la circular

6. Pandeo en sección variable

En ocasiones, a la hora de dimensionar diferentes columnas nos encontramos con el caso de columnas de sección variable, y por tanto de inercia variable, en los que no se puede aplicar de la misma forma el cálculo del pandeo.

Un ejemplo de secciones variables pueden ser soportes diseñados con un solo perfil o con varios perfiles enlazados entre sí de manera continua, y con la sección con un área ligeramente variable.

El cálculo de estos soportes se realiza obteniendo un valor del radio, y determinado como se indica en el siguiente apartado

6.1 Pandeo con R(x) variable

A continuación, se va a analizar el pandeo de un cambio de sección en la columna planteando el caso de radio variable más sencillo: un polinomio

Vamos a comprobar que R(x) se puede escribir en función de un único parámetro b

6.1.1 Gráfica de R(x) para diferentes valores de b

Como R(x) es variable, vamos a tener distintos valores del radio en función de b como hemos indicado en el apartado 1.1. Ahora, se va a dibujar numéricamente los radios de las curvas de [0,L], estudiaremos la gráfica par un intervalo de b [-10,10] y para valores positivos de b[0,10]

%Dibujar R(x) para diferentes valores de b
%Parte 1. Condiciones de Contorno
xa=0;
xb=9;
ya=0;
yb=0;
%Parte 2. Paso y tamaño de paso 
N=100;
h=(xb-xa)/N;
x=xa:h:xb;
%Parte 3. 
%%Matriz K
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
%%Vector de longitudes (Radio de 1 a 5)
R=linspace(1,5,N+1);
b=linspace(-10,10,N+1);
L=9;
%Parte 4. Bucle
 
 for k=1:length(b)
      a=(((-4.*b(k)*((L/2)).^3)/3)-sqrt((4.*b(k)*(L/2).^3)/3).^2+((4.*(L/2)^5)/5).*((25.83/pi)+b(k).^2*L))/((4*(L/2)^5)/5)
  for j=1:length(x)
      r(k,j)=a.*((x(j)-L/2)^2)+b(k);
  end
 end
%Dibujo del problema de contorno
plot(b,r)

Como podemos comprobar en la figura 6 y 7 tendremos un radio distinto para distintos valores de b. Sin embargo, para cálculos posteriores solo se tendrá en cuenta los valores positivos. También se puede comprobar que la gráfica de radios siempre será positiva para cualquier valor de b.

6.2 Valores de la carga crítica para los diferentes valores de b

Otro de los estudios interesantes es comprobar para que valores de R(x) la carga Crítica será máxima o mínima. En el cálculo de piezas de hormigón resulta de importancia, ya que podremos saber cuál es la carga máxima que necesitaremos aplicar a una columna para evitar que produzca pandeo.

A continuación, se detallará el programa que se ha utilizado para su cálculo, hemos tomado valores

clear all clc
%Calcular numéricamente por diferencias finitas 
%los valores de P(cr). Para ellos se realizará
%un problema de autovalores
 
%Parte 1. Condiciones de Contorno
xa=0;
xb=9;
ya=0;
yb=0;
;
%Parte 2. Paso y tamaño de paso 

h=0.1;
N=(xb-xa)/h;
x=xa:h:xb;
b=linspace(1,10,N+1)
L=9;
%Parte 3. 
%%Matriz K
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));


 for k=1:length(b)
      a=(((-4.*b(k)*((L/2)).^3)/3)-sqrt((4.*b(k)*(L/2).^3)/3).^2+((4.*(L/2)^5)/5).*((25.83/pi)+b(k).^2*L))/((4*(L/2)^5)/5);
  for j=1:length(x)
      r(k,j)=a.*((x(j)-L/2)^2)+b(k);
  end

    for i=1:length(r)
    A(i,:)=eig(1.2.*0.25*pi.*r(i)^4.*K);
%La matriz A es EI/M*K
    end
 end
p=A(:,1);
figure
plot(b,r)
figure
plot(r,p)

A partir de este programa, hemos creado los distintos valores máximos y mínimos de las cargas críticas.

Se puede deducir que para radios pequeños la carga crítica será menor, y para radios mayores la carga crítica será mayor. Esto es porque tanto la carga crítica como la inercia dependen del radio, y al ir multiplicados, radios mayores darán cargas críticas mayores.

También se puede comprobar que en el apartado 4 de este trabajo, se llega la misma conclusión, comparando los momentos de inercia del cuadrado y del rectángulo, solo que en este caso hay que añadir que el radio es otra variable y no es constante, es variable.

'6.2.1 Valores de la carga crítica en función del parámetro b

Vamos a representar de las cargas críticas dependiendo del valor b , así como la gráfica de cargas críticas para n=1 para diferentes valores de b

clear all clc

 
%Parte 1. Condiciones de Contorno
xa=0;
xb=9;
ya=0;
yb=0;
;
%Parte 2. Paso y tamaño de paso 

h=0.1;
N=(xb-xa)/h;
x=xa:h:xb;
b=linspace(1,10,N+1);
L=9;
%Parte 3. 
%%Matriz K
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));


 for k=1:length(b)
      a=(((-4.*b(k)*((L/2)).^3)/3)-sqrt((4.*b(k)*(L/2).^3)/3).^2+((4.*(L/2)^5)/5).*((25.83/pi)+b(k).^2*L))/((4*(L/2)^5)/5);
  for j=1:length(x)
      r(k,j)=a.*((x(j)-L/2)^2)+b(k);
  end

    for i=1:length(r)
    A(i,:)=eig(1.2.*0.25*pi.*r(i)^4.*K);
%La matriz A es EI/M*K
    end
 end
p=A(:,1);
figure
plot(b,r)
figure
plot(b,A)
figure
plot(b,p)
figure
plot(r,p)

Observamos que en ambas gráficas se cumple que a mayor valor de b mayor carga crítica.

Al igual que se ha hecho en apartados anteriores hemos elaborado un gráfico que represente los radios ( que en este caso dependerán de x , R(x)) y las cargas críticas para n=1

Estudiando la figura 11 llegamos a la misma conclusión que anteriormente, un mayor radio implica una mayor carga crítica.

6.2.2 Valor de R para Prc máximo y mínimo

Hemos hallado los valores mínimos y máximos de r, R(x) y a Pcr son directamente proporcionales, por ello hemos utilizados los comando max y min para obtener estos valores.

xa=0;
xb=9;
ya=0;
yb=0;

%Parte 2. Paso y tamaño de paso 

h=0.1;
N=(xb-xa)/h;
x=xa:h:xb;
b=linspace(1,10,N+1);
L=9
%Parte 3. 
%%Matriz K
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
%%Vector de longitudes (Radio de 1 a 5)

 for k=1:length(b)
      a=(((-4.*b(k)*((L/2)).^3)/3)-sqrt((4.*b(k)*(L/2).^3)/3).^2+((8.*(L/2)^5)/5).*((25.83/pi)+b(k).^2*L))/((4*(L/2)^5)/5);
  for j=1:length(x)
      r(k,j)=a.*((x(j)-L/2)^2)+b(k);
  end 
 end

    for i=1:length(r)
    H(i,:)=eig(1.2.*0.25*pi.*r(i)^4.*K);

    end
     
p=  H(:,1);
  for i=1:length(r)
      for i=1:length(r)
      Rmin(i)=min(r(i,:));
      end
    for i=1:length(r)
      Rmax(i)=max(r(i,:));
     end 
    A(i,:)=eig(1.2.*0.25*pi.*Rmin(i)^4.*K);
    P(i,:)=eig(1.2.*0.25*pi.*Rmax(i)^4.*K);
    end
   t= A(:,1);  
    m= P(:,1);  


figure(1)
plot(min(Rmin),min(t),'*')
xlabel(' Rmin')
ylabel(' Prc minimo')
figure(2)
plot(max(Rmax),max(m),'*')
xlabel('Rmax')
ylabel(' Pcr maximo')


rmin=(min(Rmin))
prcmin=(min(t))
rmax=(max(Rmax))
pcrmax=(max(m))
G=[t,m]

La carga crítica máxima se alcanza en 2.0967e+17 con un radio máximo de 3.6760e+04

La carga crítica mínima se alcanza en 0 con un radio minimo de 1

Hemos representados las cargas mínimas y máximas en dos columnas ( su longitud es de 91x1) pero al repetirse los valores sólo hemos hecho captura a lo más relevante

Para comparar la ganancia ( pcr con radio variable entre pcr con radio constante) se ha creado un programa matlab

xa=0;
xb=9;
ya=0;
yb=0;
;
%Parte 2. Paso y tamaño de paso 

h=0.1;
N=(xb-xa)/h;
x=xa:h:xb;
b=linspace(1,10,N+1);
L=9;
%Parte 3. 
%%Matriz K
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)
B(i,:)=eig(1.2.*(0.25)*pi*R(i)^4.*K);
%La matriz A es EI/M*K
end
p=B(:,1);

 for k=1:length(b)
      a=(((-4.*b(k)*((L/2)).^3)/3)-sqrt((4.*b(k)*(L/2).^3)/3).^2+((4.*(L/2)^5)/5).*((25.83/pi)+b(k).^2*L))/((4*(L/2)^5)/5);
  for j=1:length(x)
      r(k,j)=a.*((x(j)-L/2)^2)+b(k);
  end

    for i=1:length(r)
    A(i,:)=eig(1.2.*0.25*pi.*r(i)^4.*K);
%La matriz A es EI/M*K
    end
 end
t=A(:,1);
G=100*(t./p)

El programa nos saca por pantalla el valor ( %) de la ganancia