Diferencia entre revisiones de «Carga crítica en una columna (G 7A)»

% 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 all; clf
a=0;L=9;
ya=0;yL=0;
h=0.1;N=(L-a)/h;
x=a:h:L;
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);
l=sqrt(pi.*R.^2);       %lado del cuadrado (resultado de igualar masas)
for i=1:length(l)
    A(i,:)=eig((1/12)*1.2*l(i)^4*L.*K);     %momento inercia prisma rectangular
end
hold on
Pcr=A(:,1);
plot(l,Pcr,'-')
xlabel('Radio (m)');
ylabel('Carga critica (kN)')
hold off

clear all  clf
 
xa=0;xb=9;
ya=0;yb=0;
L=xb;
p1=1.2;
N=20; 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
a0=0.13;
af=0.5;
a=linspace(a0,af,N+1);
%Es un bucle doble anidado.
for j=1:length(a)  
    b=((-4/a(j))*sin(a(j)*(L/2))+sqrt((16/(a(j)^2))*(sin(a(j)*L/2))^2)-4*L*(((1/(2*a(j)))*sin(a(j)*L))-31/(pi*p1)+L/2));
    for k=1:length(x)
        R(j,k)=cos(a(j)*(x(k)-L/2))+b;
    end
end
for i=1:length(R)
    Rsm(i)=min(R(i,:));
end
for i=1:length(R)
    A(i,:)=eig(0.25*pi*1.2*Rsm(i).^4.*K); %Los autovalores
end
p=A(:,1); %De aquí sacamos que la columna con menor P es la número 1 (a0) y la de mayor la número N+1 (af)
figure (1) %Dibujamos la columna mas estrecha
plot(x,R(1,:),'-') %En este caso a es cte, por lo que la sección es como la resuelta en el 2 circular y 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(N+1,:),'-')
title('Columna cerca de que mayor carga crítica tiene')
xlabel('x (m)'); ylabel('Radio de la columna (m)')
grid on
figure (3)
plot(a,p,'x')
title('Evolución de la carga crítica con el valor de a')
xlabel('a'); ylabel('P_cr')
grid on
R2=0.955859;   % radio de la columna constante si su masa es 31 kg

for k=1:length(p)
    g(k)=(p(k)-(pi^3*p1*R2^4)/(4*L^2))*100/(pi^3*p1*R2^4)/(4*L^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 (4) %Dibujamos el cambio en porcentaje de Pcrit para cada valor de a respecto a la columna de seccion circular constante 
plot(a,g,'.-')
grid on
title('Ganancia en porcentaje de Pcrit respecto a la columna circular de seccion constante')
xlabel('R0'); ylabel('Ganancia en porcentaje de Pcrit')
title('Columna con R=0.9558 (31 kg r=cte)')
xlabel('a'); ylabel('Radio de la columna (m)')
figure(5)   %como cambia el radio en funcion de los valores de a y x
[xx,aa]=meshgrid(x,a);
surf(x,a,R)

clear all clc
%apartado 8
a=0;b=9;
h=0.001;
N=(b-a)/h;
x=a:h:b;
E=1;
I=1/4*1.2*pi*0.8^4;
P=0.8846678;
%introducimos la función -x/EI 
f=-x./(E*I);
%indica numero de sumando que tomamos de la serie de fourier
G=[3,5,10,40];
Q=0;
hold on
for i=1:4
    for k=1:G(i)
   
        %introducimos los autovalores
        muk=(k*pi/b)^2;
        %autofunciones
        pk=sin(k*pi*x/b);
        %calculmo los coeficientes de fourier
        ak=trapz(x,f.*pk)/trapz(x,pk.^2);
        %calculo de la aproximación de fourier
        Q=Q+(ak/(-muk+P/E*I))*pk;
    end
    plot(Q,x)
    ylabel('x (m)'); xlabel('y(x)')
end
hold off