CARGA CRÍTICA DE UNA COLUMNA
| Trabajo realizado por estudiantes | |
|---|---|
| Título | CARGA CRÍTICA DE UNA COLUMNA (Grupo-1) |
| Asignatura | Ecuaciones Diferenciales |
| Curso | |
| Autores |
Javier Marrero Patrón Tejanni El Bannoudi Guanxiong Chen Fernando Díaz-Roncero González |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 INTRODUCCIÓN
En el análisis de la Mecánica Estructural es fundamental el estudio del equilibrio de un sólido. El presente trabajo está orientado a describir cualitativamente y cuantitativamente la inestabilidad elástica. Particularizando en el caso de una pieza homogénea de sección transversal circular sometida a compresión centrada. De esta forma se plantea un acercamiento aproximado a los fenómenos de pandeo (flexión lateral) que producen fallos frágiles en elementos esbeltos. Los aspectos anteriormente mencionados tienen múltiples aplicaciones en el Cálculo de Estructuras, debido a que la aparición de dicha flexión lateral, su rápido crecimiento y la pérdida total de estabilidad del elemento produce el colapso.
2 MODELIZACIÓN MATEMÁTICA
Para poder hacer un modelo matemático del comportamiento de la columna, se tiene que considerar un pieza ideal, es decir, la rigidez a flexión (EI) tiene que ser constante a lo largo de la columna, y también es indispensable considerar el comportamiento elástico de esta, entre otra características.
Teniendo en cuanta las propiedades del material (E), de la sección (I), y curvatura que sufre la columna al ser sometido a una acción, se define la ecuación de la curva elástica, que nos permite analizar la manera con la que se deforma la pieza.
A continuación se ha estudiado un columna sometido a una carga P constante, de longitud L y suponiendo el comportamiento de la columna como una pieza ideal. En la imagen se puede observar el comportamiento general de la columna al esta sometida a la carga critica.
Para estudiar el comportamiento de la columna, se planteará un problema de contorno Partiendo de la ecuación de la curva elástica aplicada a la columna:
\begin{matrix} y(x)=\frac{M(x)}{E I(x)} \end{matrix}
Siendo:
· [math]E[/math]: Módulo de elasticidad ó de Young.
· [math]I(x)[/math]: Momento de inercia de una sección trasversal de la columna respecto al eje "Z" (se define como eje "Z" el vector perpendicula al plano "OXY").
· [math]M(x)[/math]: Momento flector.
En el caso de la columna, el momento flector depende de la deflexión, de manera que: [math]M(x)=-Py(x)[/math].
Por tanto, sustituyendo la fórmula del momento flector en la ecuación inicial: \[\begin{array}{crl}
\\
y(x)=\frac{M(x)}{E I(x)} \Longrightarrow y(x)=\frac{-Py(x)}{E I(x)} \Longrightarrow y(x)+\frac{Py(x)}{E I(x)}=0 \\
\\
\end{array}\]
Es necesario determinar las condiciones de contorno para particularizar este problema: [math] \left\{\begin{matrix} y(0)=0 \\ y(L)=0\\ \end{matrix}\right. [/math] ya que en los extremo de la viga: [math]x=0[/math] y [math]x=L[/math] tiene deflexión cero.
Por tanto, el problema de contorno planteado para resolver el ejercicio sería:
[math] \left\{\begin{matrix} y''(x)+\frac{Py(x)}{E I(x)}=0 \\ y(0)=0 \\ y(L)=0\\ \end{matrix}\right. [/math]
Siendo:
· [math]y[/math]: El desplazamiento de la columna según el eje "y".
· [math]y''(x)[/math]: Es la curvatura que adopta la columna al pandear.
3 ESTUDIO DE LA ESTABILIDAD DE LA COLUMNA
A continuación se va a estudiar la estabilidad de la columna. Se considera que la columna es estable, si la única solución del problema de contorno es la solución trivial ([math]y(x)=0[/math]), esto se va a resolver tanto de forma numérico (haciendo uso de Matlab), y de forma analítica:
Para calcular la carga que hace pandear la columna debemos resolver la problema de contorno planteado anteriormente:
(Nota: para facilitar los cálculos, consideramos que: [math]y(x)=y[/math], [math]I(x)=I[/math]).
Por tanto, el problema quedaría
\[\begin{array}{crl}
\\
E Iy+Py=0 \\
\\
\end{array}\]
Resolución:
Al ser la ecuación diferencial homogénea y el coeficiente de y es independientes de "x", se puede hallar directamente el polinomio característico:
\[\begin{array}{crl}
\\
E Iy+Py=0 \Longrightarrow E Im^2+P=0 \Longrightarrow m^2=\frac{-P}{E I}=0 \Longrightarrow m=\; \pm\;\sqrt{\frac{P}{E I}} i \\
\\
\end{array}\] La solución general de la ecuación diferencial homogénea será \[\begin{array}{crl}
\\
y\;= y(x)\; =\; C_{1}\; cos\; \left(\sqrt{\frac{P}{E I}} x\right) + \; C_{2}\; sin\; \left(\sqrt{\frac{P}{E I}} x\right) = 0\\
\\
\end{array}\]
Aplicando las condiciones de contorno iniciales ya mencionadas, se particulariza la ecuación [math]y(0)=0[/math] y [math]y(L)=0[/math]
Para [math]y(0)=0[/math]
\[\begin{array}{crl}
\\
y(0)\; =\; C_{1}\; cos\; (0) + \; C_{2}\; sin\; (0) = 0 \Longrightarrow y(0)\; =\; C_{1}\; 1 + \; C_{2}\; 0 = 0\\
\\
\end{array}\]
Sabiendo que la solución de [math]y(0)[/math] tiene que ser 0, se llega a la conclusión de que, para que se cumpla la ecuación, tiene que ocurrir que [math]C_{1}=0[/math] necesariamente.
Por otro lado, para [math]y(L)=0[/math] \[\begin{array}{crl}
\\
y(L)\; =\; 0\; cos\; \left(\sqrt{\frac{PL^2}{E I}}\right) + \; C_{2}\; sin\; \left(\sqrt{\frac{PL^2}{E I}}\right) = 0 \Longrightarrow y(L)\; =\; 0\ + \; C_{2}\; sin\; \left(\sqrt{\frac{PL^2}{E I}}\right) = 0\\
\\
\end{array}\]
Como se puede comprobar, [math]y(L)=0[/math], implica que [math]\; C_{2}\; sin\; \left(\sqrt{\frac{PL^2}{E I}}\right) = 0[/math]. Si [math]C_{2}=0[/math], se tiene [math]y=0[/math], pero, si [math]C_{2}≠0[/math], entonces [math]sin\; \left(\sqrt{\frac{P}{E I}} x\right) = 0[/math]. La última condición indica que, para que la igualdad planteada se cumpla, el argumento de la función seno ha de ser un múltiplo entero de p. \[\begin{array}{crl}
\\
\sqrt{\frac{P}{E I}} L= nπ \Longrightarrow {\frac{P}{E I}}= {\frac{n^2π^2}{L^2}} \\
\\
\end{array}\] Por lo tanto, para todo real [math]C_{2}[/math] distinto de cero, es una solución del problema para cada n. Puesto que la ecuación diferencial es homogénea, y ya hemos supuesto que [math]sin\; (\frac{nπx}{L})= 0[/math], no necesitamos escribir [math]C_{2}[/math] si así lo deseamos, por lo que la ecuación nos quedaría \[\begin{array}{crl}
\\
y = sin\; \left(\frac{nπx}{L}\right) \\
\\
\end{array}\]
Volviendo a la ecuación anteriormente obtenida, se hallará el valor de [math]P_{cr}[/math] donde la columna se desvía. Esto sólo ocurre cuando la carga toma determinados valores de dicha ecuación (que coinciden con los autovalores de la ecuación). Esas fuerzas se llaman cargas críticas, que corresponden a la ecuación
\[\begin{array}{crl}
\\
\sqrt{\frac{P_{cr}}{E I}} L= nπ \Longrightarrow {\frac{P_{cr}}{E I}}= {\frac{n^2π^2}{L^2}} \Longrightarrow P_{cr}\; = {\frac{n^2π^2EI}{L^2}} \\
\\
\end{array}\]
De estos autovalores sólo no interesa el mas pequeños, ya que es el primer valor que provoca el pandeo de la columna. Como se puede observar en la formula de la carga critica, el menor valor se obieteine para [math]n=1[/math], ya que para [math]n=0[/math] la columna es estable.
A continuación se va exponer un ejemplo para la mejor comprensión del problema.
Se considera un columna de longitud [math]L=8m[/math], de una sección circular de radio constante [math]R=0.3m[/math], de modulo [math]E=27000MPa[/math] y cuya densidad es [math]ρ=2400 Kg/m^3[/math].
Para calcular la carga crítica necesitamos el momento de inercia de una sección circular respecto al eje "Z": [math]I = {\frac{πR^4}{4}}[/math].
Por tanto, la nueva ecuación de [math]P_{cr}[/math] quedaría expresada así \[\begin{array}{crl}
\\
P_{cr}\; = {\frac{n^2π^3ER^4}{4L^2}} \\
\\
\end{array}\]
Sustituyendo los datos del ejemplo en la formula se obtiene la siguiente carga crítica:
\[\begin{array}{crl}
\\
P_{cr}\; = {\frac{1^2π^30.3^4·2.7·10^6}{4·8^2}=26489kN} \\
\\
\end{array}\]
Para hacer un estudio mas detallado del comportamiento de la carga crítica, en función de la sección transversal de la columna, se va a variar el radio entre 0.1-0.5m, para obtener la representación gráfica de lo dicho, se va a hacer uso de Matlab.
%Características General de la columna
E=27000*10^3; %[KN/m^2]
L=8; %Longitud de la Columna[m]
R=0.1:0.1:5; %Variación del Radio
I=(pi*R.^4)./4; %Inercia de la Sección Tranversal(Circular)
%Se considera el menor valor posible de la carga crítica n=1
Pc=(pi^2*E*I)./L^2;
plot(R,Pc) %Representación gráfica de la carga en función del radio
xlabel('RADIO(m)')
ylabel('CARGA CRÍTICA(kN)')
Como se puede apreciar en el gráfico, a un mayor radio conlleva un aumento de la inercia por lo tanto mayor carga crítica.
Para contratar el método de resoluciones realizado anteriormente, se va a resolver de nuevo, pero esta vez numéricamente.
Para ello se va a hacer uso de una aproximación por diferencias finitas del problema, que consiste en plantear el problema matricialmente y calcular una aproximación de la carga crítica usado un tamaño de paso 0.1.
El problema de contorno planteado anteriormente es: [math] \left\{\begin{matrix} EIy''(x)+Py(x)=0 \\ y(0)=0 \\ y(L)=0\\ \end{matrix}\right. [/math]
Aproximamos: y"≈KY , donde K tiene dimensión de matriz.
Sustituyendo y" en el problema de contorno obtenemos: E·I·K·Y+PY=0
Despejando PY y definiendo A=EIK, se obtiene: [math](A-λI)Y=0[/math]
Se ha reducido el problema de contorno a un problema de autovalores, para la resoluciones de dichos autovalores en Matlab se ha usado el comando "eig", como se puede observar en el código de Matlab.
a=0; %Extremo Izquierdo (a) [Límite Inferior] del Intervalo
b=8; %Extremo Derecho (b) [Límite Superior] del Intervalo
ca=0; %Condición de Contorno en x=a [Condición Inicial y(a)]
cb=0; %Condición de Contorno en x=b [Condición Final y(b)]
%Discretización Espacial
h=0.1; %Tamaño de paso
%Número de Subintervalos
N=round((b-a)/h);
%Matriz K - Valores de la aproximación de la segunda derivada
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
R=linspace(0.1,5,N-1); %Variación del Radio
%Determinación de Autovalores (A=E*I*K)
for i=1:length(R)
A(i,:)=eig(27000*10^3*(pi*R(i).^4)./4.*K);
end
%Determinación de la Carga Crítica
Pc=A(:,1);%Solución Numérica
%Características Generales
%Material=Hormigón
E=27000*10^3; %[KN/m^2]
L=8; %Longitud de la Columna[m]
R1=linspace(1,5,N-1); %Variación del Radio
I1=(pi*R1.^4)./4; %Inercia de la Sección Tranversal(Circular)
%Se considera el menor valor posible de la carga crítica n=1
Pc1=(pi^2*E*I1)./L^2;%Solución Exacta
hold on
plot(R,Pc,'r')
plot(R1,Pc1,'b')
hold off
xlabel('RADIO(m)')
ylabel('CARGA CRÍTICA(KN)')
En el imagen anterior se ha dibujado la evolución de la carga crítica calculada analíticamente y numéricamente, como se puede observar el error del cálculo numérico es mínimo, siendo mas exacto en cálculo analítico.
A continuación se va a volver a realizar el procedimiento numérico realizado anteriormente con la única diferencia que ahora la sección transversal de la columna será una sección cuadrada de masa constante e igual a la masa de la columna con sección circular, para sacar el área de esta sección cuadrada vamos a realizar lo siguiente :
Como el material es el mismo se pueden igualar sus densidades siendo la formula de la densidad:
\begin{matrix} ρ=Masa/Volumen \end{matrix}
Al igualar las densidades se puede apreciar que sus áreas son iguales, al tener los dos la misma longitud (L=8m) y la misma masa.
\begin{matrix} π·R^2=b^2 \end{matrix} siendo: [math] b=R·sqrt{π} [/math]
La carga carga crítica para esta nueva sección quedaría: \[\begin{array}{crl}
\\
P_{cr}\; = {\frac{n^2π^2·E·b^4}{12·L^2}} \\
\\
\end{array}\]
Sustituyendo los datos de la columna: \[\begin{array}{crl}
\\
P_{cr}\; = {\frac{1^2π^2·2.7·10^6·0.3·sqrt{π}^4}{12·8^2}=27738kN} \\
\\
\end{array}\]
A continuación se va a exponer el código empleado en matlab para la resolución numérica y su comparación con la analítica para el caso de la sección cuadrada
a=0; %('Extremo Izquierdo (a) [Límite Inferior] del Intervalo: ');% 1
b=8; %('Extremo Derecho (b) [Límite Superior] del Intervalo: ');% 2
ca=0;%('Condición de Contorno en x=a [Condición Inicial y(a)] ');% 1
cb=0;%('Condición de Contorno en x=b [Condición Final y(b)] ');% 3
%Discretización Espacial
h=0.1; %Tamaño de paso
%Número de Subintervalos
N=round((b-a)/h);
%Matriz K - Valores de la aproximación de la segunda derivada
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
R=linspace(0.1,5,N-1); %Variación del Radio
b=sqrt(pi)*R; %Lado de la sección cuadrada
%Determinación de Autovalores (A=E*I*K)
for i=1:length(b)
A(i,:)=eig(27000*10^3*(b(i).^4)./12.*K);
end
%Determinación de la Carga crítica
Pc=A(:,1);
%Características Generales
%Material=Hormigón
E=27000*10^3; %[KN/m^2]
L=8; %Longitud de la Columna[m]
R1=linspace(1,5,N-1); %Variación del Radio
b1=sqrt(pi)*R1; %Lado de la sección cuadrada
I1=(b1.^4)./12; %Inercia de la Sección Tranversal(Cuadrada)
%Se considera el menor valor posible de la carga crítica n=1
Pc1=(pi^2*E*I1)./L^2;
hold on
plot(R,Pc,'r')
plot(R1,Pc1,'b')
hold off
xlabel('RADIO(m)')
ylabel('CARGA CRÍTICA(KN)')
end
Una vez calculada la carga crítica para la sección cuadrada se la va a comparar con la de la sección circular, para ello necesitamos crear un nuevo código donde aparezcan las dos secciones, para poder representarlas juntas, el código es el siguiente:
E=27000*10^3; %[KN/m^2]
L=8; %Longitud de la Columna[m]
R=1:0.1:2; %Variación del Radio
I=(pi*R.^4)./4; %Inercia de la Sección Tranversal(Circular)
%Se considera el menor valor posible de la carga crítica n=1
Pc=(pi^2*E*I)./L^2;
R1=linspace(1,2,N-1); %Variación del Radio
b1=sqrt(pi)*R1; %Lado de la sección cuadrada
I1=(b1.^4)./12; %Inercia de la Sección Tranversal(Cuadrada)
%Se considera el menor valor posible de la carga crítica n=1
Pc1=(pi^2*E*I1)./L^2;
hold on
plot(R,Pc,'r')
plot(R1,Pc1,'b')
xlabel('RADIO(m)')
ylabel('CARGA CRÍTICA(KN)')
legend('CIRCULAR','CUADRADA')
Ejecutando el programa nos sale la siguiente comparativa:
Como se puede apreciar en el gráfico anterior la columna con sección cuadrada será mejor, debido a que posee mayor inercia y por tanto mayor rigidez a la flexión lateral(pandeo). La diferencia de cargas críticas entre una y otra aumenta a mayor radio precisamente por el efecto de la inercia, que es único parámetro que no es constante.
Es decir, el momento de inercia de sección circula es: [math]I=6.3617·1^-3m^4[/math], por lo tanto la carga critica es:[math]P_{cr}=26489kN[/math], mientra que el momento de inercia de la sección cuadrada es [math]I=6.662·1^-3m^4[/math], y dando una carga critica de:[math]P_{cr}=27738kN[/math].
4 ANÁLISIS DE UNA COLUMNA DE SECCIÓN VARIABLE
Ahora vamos a suponer que la sección circular de la columna depende de "x", según la siguiente fórmula: [math]R(x)=a·|x-L/2|+b[/math] También se considera la masa constante, [math]m=1000Kg[/math]. A continuación vamos calcular la sección de la columna siguiendo la fórmula anterior, dando diferentes valores al parámetro b. Para ello es necesario poner el radio sólo en función de b. Como la masa es constante para todas las columnas que podemos formar variando b y su longitud no varia, entonces se puede poner la fórmula del radio en función de un sólo parámetro (en este caso b).
En el desarrollo anterior se ha considerado solo la raíz positiva, ya que para raíces negativas el radio da valores negativos.
Para poder dibujar R(x) para diferentes valores de b se a recurrido a Matlab cuyo código es el siguiente:
a=0; %('Extremo Izquierdo (a) [Límite Inferior] del Intervalo: ');% 1
b=8; %('Extremo Derecho (b) [Límite Superior] del Intervalo: ');% 2
%Discretización Espacial
h=0.5; %Tamaño de paso
%Número de Subintervalos
N=round((b-a)/h);
x=linspace(a,b,N-1);
b=linspace(0,0.22,N-1); %Variación de b
L=8; %Longitud de la Columna[m]
rho=2400; %Densidad [kg/m^3]
a=(-b.*L^2/2+sqrt((b.*L^2/2).^2-L^3/3.*(b.^2.*L-1000/(rho*pi))))/(L^3/6);
hold on
for j=1:length(b)
R=a(j)*abs(x-L/2)+b(j);
plot(x,R)
pause
end
hold off
xlabel('x')
ylabel('RADIO R(x)')
Este programa da la siguiente representación gráfica:
En la animación se muestra como varía la geometría de la columna cuando se considera un cambio de la sección. Cada gráfica representa la forma que obtiene la columna en función del parámetro "b". En el eje de las abcisas se representa los valores (en metros) de la posición y en el eje de las ordenadas, el valor del radio (en metros).
Para finalizar el profundizaje del estudio del pandeo de columnas, se va a representar en una gráfica el valor de P_{cr} para diferentes valores de b, el código para esta representación es el siguinte:
a=0; %('Extremo Izquierdo (a) [Límite Inferior] del Intervalo: ');% 1
b=8; %('Extremo Derecho (b) [Límite Superior] del Intervalo: ');% 2
%Discretización Espacial
h=0.1; %Tamaño de paso
%Número de Subintervalos
N=round((b-a)/h);
x=linspace(a,b,N-1);
b=linspace(0.1,0.22,N-1); %Variación de b
L=8; %Longitud de la Columna[m]
rho=2400; %Densidad [kg/m^3]
%Matriz K - Valores de la aproximación de la segunda derivada
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
a=(-b.*L^2/2+sqrt((b.*L^2/2).^2-L^3/3.*(b.^2.*L-1000/(rho*pi))))/(L^3/6);
hold on
for j=1:length(b)
R=a(j)*abs(x-L/2)+b(j);
%Determinación de Autovalores (A=E*I*K)
for i=1:length(R)
A(i,:)=eig(27000*10^3*(pi*R(i).^4)./4.*K);
end
plot(x,A(:,1))
pause
end
hold off
xlabel('x')
ylabel('CARGA CRÍTICA(KN)')
Al ejecutar este programa se puede observar
En la animación se muestra la variación de la carga crítica para todos los posibles puntos de la columna en sentido longitudinal (secciones transversales), considerando que cada gráfica corresponde a un valor específico de la constante"b". En el eje de las abcisas se representa los valores (en metros) de la posición y en el eje de las ordenadas, la carga crítica en kN.
