[FEM] Barra empotrada en Matlab

Ok, esto es algo que se me planteó en la Universidad, en el ramo de Elementos Finitos.

Si bien es cierto no había que hacerlo en MATLAB, pero gracias al ramo de Elementos Finitos Generalizados he empezado a tomarle el gusto a este software.

El problema:

Sea una barra empotrada en uno de sus extremos, mientras que en el otro, se aplica una fuerza P=300 [N]. Las dimensiones según la Figura 1 son:

  • P=300 [N]
  • l=1 [m]
  • h=b=0,02 [m]
  • E=2,07E11 [N/m]

Solución

Bien, aquí es algo complejo. Me basaré en la siguiente literatura:

  • The Finite Element Method a Practical Course – G.R. Liu & S. S. Quek [ON LINE]

Código
function barra(ne,largo)
nxe=2;          % 2 Nodos por elemento
nodos=ne+1;     % Total de nodos
le=largo/ne;    % Largo de cada elemento
dof=2;          % Grados de libertad por nodo
doft=dof*nodos; % Grados de libertad del sistema
% %
% Defino las funciones de forma
syms xi;
N1=1/4*(2-3*xi+xi^3);
N2=le/8*(1-xi-xi^2+xi^3);
N3=1/4*(2+3*xi-xi^3);
N4=le/8*(-1-xi+xi^2+xi^3);
N=[N1 N2 N3 N4];
NT=[N1;N2;N3;N4];
% % Defino la matriz de desplazamiento-deformacion
B=-diff(N,xi,2);
BT=-diff(NT,xi,2);
% %
EI=2760;
K_global=zeros(doft);
Q_global=zeros(doft,1);
% Para el elemento e
for e=1:ne
    % % Defino las matrices de conectividad
    Ae=zeros(dof*nxe,doft);
    i=1;
    while i<=dof*nxe
        Ae(i,dof*e-2+i)=1;
        i=i+1;
    end% Fin Matrices de Conectividad % % Defino la matriz de rigidez para el elemento e
    Ke=EI/(le/2)^3*int(BT*B,xi,-1,1);
    % % Comienzo el ensamble de la matriz global
    K_global=K_global+Ae'*Ke*Ae;
end % Finalizo el calculo para cada elemento % %Imposicion de condiciones de borde
[K,Q]=boundary_conditions(K_global,Q_global);
% %Resolución del sistema
solucion=inv(K)*Q;
% Comienzo el ploteo de la solución
clf('reset');
[sp,rad,xsp]=sol_printed(ne,largo,nxe,solucion);
subplot(2,1,1);
plot(xsp, sp, '-r');
title('Deformacion de la barra, en metros.');
axis([0 largo min(sp)*2 min(sp)*-1]);
hold on;
plot(max(xsp), min(sp), 'or');
text(max(xsp)*.95, min(sp), '\rightarrow');
leyenda=legend('Deformación de la barra','Deflexión Máxima',3);
set(leyenda,'Interpreter','none');
hold on;
subplot(2,1,2);
plot(xsp, rad, '*r');
title('Variacion de la inclinacion de la barra en radianes');
hold off;
end

Descarga de barra .m

function [K,Q]=boundary_conditions(K_global,Q_global)
%Carga de 300 N
Q_global(length(Q_global)-1)=Q_global(length(Q_global)-1)-300;
%Empotramiento en el origen
Q_global(1)=0; Q_global(2)=0;
K_global(:,1)=0; K_global(1,:)=0;
K_global(:,2)=0; K_global(2,:)=0;
K_global(1,1)=1; K_global(2,2)=1;
%
K=K_global;
Q=Q_global;
end

Descarga de boundary_conditions.m

function [sp,rad,xsp]=sol_printed(ne,largo,nxe,solucion)
nodos=ne+1;
sp=zeros(nodos*nxe/2,1);
rad=zeros(nodos*nxe/2,1);
for i=1:nodos
    sp(i)=solucion(i*nxe-1);
end
for j=1:nodos
    rad(j)=solucion(j*nxe);
end
xsp=global_coord(ne,largo);
end

Descarga de sol_printed.m

function [coordglobales]=global_coord(ne,l)
% Devuelve las coordenadas globales para una barra lineal
% ne=numero de elementos en la barra
% l= largo de la barra
% nodos=cantidad total de nodos en la barra
% le=largo de cada elemento
% coordglobales(i)=coordenada global del nodo i en la barra
%
nodos=ne+1;
le=l/ne;
coordglobales=zeros(nodos,1);
for i=1:nodos
    coordglobales(i)=le*(i-1);
end
end

Descarga de global_coord.m

Resultado

ne=10 ; largo=1

Cosas por mejorar

Este método, como se puede apreciar, realiza el cálculo de las matrices de rigidez via integrales. Si bien es cierto que la teoría indica esto, resolverlas cuando hay muchos elementos puede ser lento. Así que una mejora a la rutina, sería aplicar el método de Gauss para no tener que integrar. Este método da buenos resultados, los cuales espero hacer eventualmente.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: