Ensamblaje de la matriz de rigidez global en elementos finitos (barras, MATLAB)

Cálculos Preliminares y Definiciones

A continuación, se definen las variables iniciales necesarias para el proceso de ensamblaje de la matriz de rigidez global.

• s=size(CON,1); % Número de elementos
• n=length(COOR); % Número de nodos
• N=2*n; % Número total de grados de libertad
• p=length(inddD); % Número de grados de libertad con desplazamiento conocido (prescrito)
• m=N-p; % Número de grados de libertad con desplazamiento desconocido (incógnita)

Se establece el formato de visualización para números:

format short e

Se define la matriz de rigidez de una barra en su sistema de referencia local, sin el factor EA/L. Esta matriz es fundamental para el cálculo dentro del bucle de elementos.

Kebarraprov=[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0];

Se inicializa un vector columna para almacenar las longitudes de cada barra:

LONG=zeros(s,1);

Ensamblaje de la Matriz de Rigidez Global (K)

La matriz de rigidez global K se inicializa con ceros y se procede a su ensamblaje elemento por elemento.

K=zeros(N);

El siguiente bucle itera sobre cada elemento (barra) para calcular su contribución a la matriz de rigidez global.

for e=1:s

Cálculo de Coordenadas y Longitud de la Barra

Se obtienen las coordenadas x1, y1, x2, y2 de los nodos 1 y 2 de la barra actual. (Sección pendiente de implementación)

I=CON(e,1);
J=CON(e,2);
x1=COOR(I,1);
y1=COOR(I,2);
x2=COOR(J,1);
y2=COOR(J,2);
vector1=[x1,y1];
vector2=[x2,y2];

Se calcula la longitud de la barra y se almacena en el vector LONG. (Sección pendiente de implementación)

LONG(e,1)=norm(vector1-vector2);

Cálculo de la Matriz de Rigidez Local y Matriz de Transformación

Se calcula la matriz de rigidez de la barra en su referencia local, kebarra. (Sección pendiente de implementación)

kebarra=(BARRAS(e)/LONG(e,1))*Kebarraprov;

Se calcula la matriz de transformación de coordenadas (matriz de rotación), Re. (Sección pendiente de implementación)

A=(x2-x1)/LONG(e,1);
B=(y2-y1)/LONG(e,1);
Re=[A -B; B A];

Cálculo de la Matriz de Rigidez del Elemento en Coordenadas Globales (ke)

La matriz de rigidez del elemento en coordenadas globales, ke, se obtiene mediante la transformación de la matriz de rigidez local. (Sección pendiente de implementación)

ke=[Re zeros(2,2); zeros(2,2) Re]*kebarra*[Re' zeros(2,2); zeros(2,2) Re'];

Ensamblaje de Submatrices de Rigidez en K

Se procede al ensamblaje de las submatrices de ke (ke11, ke12, ke21, ke22) en la matriz de rigidez global K. (Secciones pendientes de implementación)

Ensamblaje de ke11

ke11=ke([1,2],[1,2]);
I=CON(e,1);
J=CON(e,1);
K([2*I-1,2*I],[2*J-1,2*J])=K([2*I-1,2*I],[2*J-1,2*J])+ke11;

Ensamblaje de ke12 y ke21

ke12=ke([1,2],[3,4]);
I=CON(e,1);
J=CON(e,2);
K([2*I-1,2*I],[2*J-1,2*J])=K([2*I-1,2*I],[2*J-1,2*J])+ke12;
K([2*J-1,2*J],[2*I-1,2*I])=K([2*J-1,2*J],[2*I-1,2*I])+ke12;

Ensamblaje de ke22

ke22=ke([3,4],[3,4]);
I=CON(e,2);
J=CON(e,2);
K([2*I-1,2*I],[2*J-1,2*J])=K([2*I-1,2*I],[2*J-1,2*J])+ke22;

end

Nota: Para usuarios con mayor experiencia, es posible ensamblar la matriz ke de una sola vez, sin necesidad de separar y ensamblar ke11, ke12, ke21 y ke22 individualmente.

Resolución del Sistema de Ecuaciones Global

Se procede a la resolución del sistema de ecuaciones Kaux*d=Faux para determinar los desplazamientos d y las fuerzas F. (Sección pendiente de implementación)

Kaux=K;
Faux=zeros(N,1);
inddI=[3 4 5 6 7 8 9 10 11];
G=1e18;
Faux(inddI)=FD;
Faux(inddD)=G*dD;

Se aplica una penalización para las condiciones de contorno de desplazamiento conocido.

for i=1:p
Kaux(inddD(i),inddD(i))=G;
end

Finalmente, se resuelven los desplazamientos:

d=Kaux\Faux;