clearvars; clc
%--- Cantidad de elementos y nodos del sistema ---%
e = 6;%input('Ingresar la cantidad de elementos finitos (debe ser múltiplo de 3): '); % Número de elementos a utilizar
n = e+1; % Cantidad de nodos (siempre igual a la cantidad de elementos + 1)
%--- Excitación ---%
Fext = 1000*9.81; % Fuerza aplicada, [N]
%--- Distribución de los nodos en cada cuerpo ---%
% Se asume que se desea una distribución igualitaria de nodos en los 3 cuerpos
c = e/3 + 1; % Cantidad de nodos en cada cuerpo
%--- Parámetros de los elementos de cada cuerpo ---%
LC1 = 0.01; % Largo del cuerpo 1, [m]
LC2 = 0.008; % Largo del cuerpo 2, [m]
LC3 = 0.006; % Largo del cuerpo 3, [m]
Le(1) = LC1/(e/3); % Largo de los elementos dentro del cuerpo 1
Le(2) = LC2/(e/3); % Largo de los elementos dentro del cuerpo 2
Le(3) = LC3/(e/3); % Largo de los elementos dentro del cuerpo 3
% Largo acumulado (utilizado para plotear al final) %
Lac = zeros(1,e);
Lac(1) = Le(1);
for i=2:e/3
Lac(i) = Lac(i-1) + Le(1);
end
for i=e/3+1:2*e/3
Lac(i) = Lac(i-1) + Le(2);
end
for i=2*e/3+1:e
Lac(i) = Lac(i-1) + Le(3);
end
AC1 = pi*0.004^2/4; % Ãrea del cuerpo 1, [m2]
AC2 = pi*0.002^2/4; % Ãrea del cuerpo 2, [m2]
AC3 = pi*0.002^2/4; % Ãrea del cuerpo 3, [m2]
EC1 = 21000*9.81*10^6; % Módulo elástico del cuerpo 1, [Pa]
EC2 = 8000*9.81*10^6; % Módulo elástico del cuerpo 2, [Pa]
EC3 = 1000*9.81*10^6; % Módulo elástico del cuerpo 3, [Pa]
%--- Rigidez de los elementos de cada cuerpo ---%
keC1 = EC1*AC1/Le(1);
keC2 = EC2*AC2/Le(2);
keC3 = EC3*AC3/Le(3);
%--- Matriz de rigidez (MMRR) ---% en SI
K = zeros(n);
kcan = [1 -1;-1 1]; % Matriz canónica de rigidez
keC1 = keC1*kcan; % Matriz "canónica" de rigidez de los elementos
keC2 = keC2*kcan;
keC3 = keC3*kcan;
% Acoplamiento %
% Adición de la rigidez del cuerpo 1 a la MMRR %
for i=1:e/3+1
for j=1:e/3+1
if i==j && (i==1 || i==e/3+1)
K(i,j) = K(i,j) + keC1(1,1);
elseif i==j && i>=2
K(i,j) = K(i,j) + 2*keC1(1,1);
elseif i~=j && (i == j+1 || j == i+1)
K(i,j) = K(i,j) + keC1(1,2);
end
end
end
% Adición de la rigidez del cuerpo 2 a la MMRR %
for i=e/3+1:2*e/3+1
for j=e/3+1:2*e/3+1
if i==j && (i==e/3+1 || i==2*e/3+1)
K(i,j) = K(i,j) + keC2(1,1);
elseif i==j && i>=e/3+2
K(i,j) = K(i,j) + 2*keC2(1,1);
elseif i~=j && (i == j+1 || j == i+1)
K(i,j) = K(i,j) + keC2(1,2);
end
end
end
% Adición de la rigidez del cuerpo 3 a la MMRR %
for i=2*e/3+1:n
for j=2*e/3+1:n
if i==j && (i==2*e/3+1 || i==n)
K(i,j) = K(i,j) + keC3(1,1);
elseif i==j && i>=2*e/3+2
K(i,j) = K(i,j) + 2*keC3(1,1);
elseif i~=j && (i == j+1 || j == i+1)
K(i,j) = K(i,j) + keC3(1,2);
end
end
end
%--- Desplazamientos y fuerzas externas ---%
syms q [1 n]
syms F [1 n]
% Desplazamientos %
q(1) = 0; q(n) = 0.05e-03; % Condiciones de borde
q = q.';
% Fuerzas externas %
% Sólo existen 3 fuerzas externas en el problema: El empotramiento en el nodo 1, la fuerza aplicada en el nodo e/3+1
% y la reacción que ocurrirá en el último nodo debido a que se elonga más del espacio que existe.
for i=2:n-1
F(i) = 0;
end
F(e/3+1) = Fext;
F = F.';
%--- Condensación estática: Matriz de Parámetros Condensados, MPC ---%
K = K([2:n],[2:n]); % Se elimina la primera fila y columna de la MMMRR
q = q(2:n);
%--- Resolución de la MPC ---%
for i=1:n-1
ecn(i) = K(i,:)*q == F(i+1);
end
[solq2, solq3, solq4, solq5, solq6, solF7] = solve(ecn, [q2 q3 q4 q5 q6 F7]);
S = [solq2, solq3, solq4, solq5, solq6];
plot(Lac(2:e),S)
