AYUDA EN CODIGO DE MATLAB PARA UN ESTRUCTURA DE PUNETE EN 3D

6 次查看（过去 30 天）
Lizandro
Lizandro 2024-11-14，17:30
回答： Torsten 2024-11-14，18:17
Quiero que me funcione para porder hallar los desplazamiento y reacciones, pero me sale el siguiente error. Warning: Matrix is singular to working
precision.
> In TrabajoF (line 80)
ADJUNTO EL CÓDIGO MATLAB:
% Armadura espacial
% Módulo de elasticidad:
E = 210 * 1e6; % en N/m² (Pa)
% Área:
A = 400 / 1e6; % en m²
% Coordenadas de los nodos:
XYZ = [16 5 0; 12 5 0;8 5 0; 4 5 0; 0 5 0; 12 5 4; 4 5 4; 0 0 0; 4 0 0;...
8 0 0; 12 0 0; 16 0 0; 12 0 4; 4 0 4];
x = XYZ(:,1);
y = XYZ(:,2);
z = XYZ(:,3);
% Topología:
IJ = [1 2;2 3;3 4;4 5;1 6;2 6;3 6;3 7;4 7;5 7;6 7;12 11;11 10;10 9;...
9 8;12 13;11 13;10 13;10 14;9 14;8 14;13 14;6 13;14 7;5 8;1 12];
% Dibujo de la estructura:
figure
for e = 1:26
Q = [XYZ(IJ(e,1),1) XYZ(IJ(e,1),2) XYZ(IJ(e,1),3); ...
XYZ(IJ(e,2),1) XYZ(IJ(e,2),2) XYZ(IJ(e,2),3)];
plot3(Q(:,1), Q(:,2), Q(:,3), "-b")
hold on
end
xlabel("x")
ylabel("y")
zlabel("z")
axis equal
% Longitudes de las barras
l = zeros(26,1);
for e = 1:26
l(e) = sqrt((x(IJ(e,2))-x(IJ(e,1)))^2 + (y(IJ(e,2))-y(IJ(e,1)))^2 + (z(IJ(e,2))-z(IJ(e,1)))^2);
end
% Matriz de rigidez global
K = zeros(42,42);
for e = 1:26
eta = (x(IJ(e,2)) - x(IJ(e,1))) / l(e);
mu = (y(IJ(e,2)) - y(IJ(e,1))) / l(e);
nu = (z(IJ(e,2)) - z(IJ(e,1))) / l(e);
% Matriz de rigidez local del elemento
K_e = (E * A / l(e)) * [eta^2 eta*mu eta*nu -eta^2 -eta*mu -eta*nu; ...
eta*mu mu^2 mu*nu -eta*mu -mu^2 -mu*nu; ...
eta*nu mu*nu nu^2 -eta*nu -mu*nu -nu^2; ...
-eta^2 -eta*mu -eta*nu eta^2 eta*mu eta*nu; ...
-eta*mu -mu^2 -mu*nu eta*mu mu^2 mu*nu; ...
-eta*nu -mu*nu -nu^2 eta*nu mu*nu nu^2];
% Grados de libertad globales del elemento
g_e = [3*IJ(e,1)-2 3*IJ(e,1)-1 3*IJ(e,1) 3*IJ(e,2)-2 3*IJ(e,2)-1 3*IJ(e,2)];
% Expande K_e en la matriz de rigidez global K
DeltaK_e = zeros(42,42);
DeltaK_e(g_e, g_e) = K_e;
K = K + DeltaK_e;
end
% Definición de desplazamientos y reacciones
% Grados de libertad restringidos y libres
a = (5:42)'; % Grados de libertad libres
b = (1:4)'; % Grados de libertad restringidos
% Dividir la matriz de rigidez
K_aa = K(a, a);
K_ab = K(a, b);
K_ba = K(b, a);
K_bb = K(b, b);
% Vector de fuerzas
P = zeros(42,1);
P(3) = -2;
P(10) =-2% Fuerza en el nodo 1, dirección z
% Calcular desplazamientos en los grados de libertad libres
U_a = K_aa \ (P(a) - K_ab * zeros(length(b), 1));
% Reacciones en los grados de libertad restringidos
R_b = K_ba * U_a + K_bb * zeros(length(b), 1);
% Mostrar resultados
disp('Desplazamientos en los grados de libertad libres:')
disp(U_a)
disp('Reacciones en los grados de libertad restringidos:')
disp(R_b)

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2024-11-14，18:17
If you insert the lines
% Calcular desplazamientos en los grados de libertad libres
rank(K_aa)
size(K_aa)
before the command
U_a = K_aa \ (P(a) - K_ab * zeros(length(b), 1));
you will see that K_aa is 38x38, but that its rank is only 25 - thus deeply rank-deficient.

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by