Demo entry 6356405

koen

   

Submitted by anonymous on Apr 18, 2017 at 12:42
Language: Matlab. Code size: 6.9 kB.

clc
clear
close all
 
%% PROBLEMA 1
%% Datos
 
E = 60e9; % Módulo de Young [Pa]
P = 9500; % Fuerza exterior P [N]
H1 = 450e-3; % [m]
H2 = 900e-3; % [m]
L = 1200e-3; % [m]
 
%% Geometría
 
N = 4; % Nodos
Nd = 2; % Dimensiones del problema
Nel = 4; % Número de elementos
Nig = 3; % Número de grados de libertad por nodo
Ngl = N * Nig; % Total de grados de libertad (Ngl = Nig * N)
Nnod = 2; % Número de nodos en una barra
 
 
% Matriz de coordenadas nodos
 
x = [0 L 0 0; (H2+H1) H2 H2 0]; 
 
% Matriz de conectividades
 
T1 = [1 1 3 3;...
      2 3 2 4]; % Matriz conectividades nodos
T2 = [1 1 7 7;...
      2 2 8 8;...
      3 3 9 9;...
      4 7 4 10;...
      5 8 5 11;...
      6 9 6 12]; % Matriz conectividades grados de libertad
  
  Ic = [0; 3e-6; 1.5e-6; 3e-6]; % Inercias de cada barra [m^4]zsw
  A = [12e-6; sqrt(Ic(2)*12); sqrt(Ic(3)*12); sqrt(Ic(4)*12)]; % Vector secciones barras
  
  %% Computación de matrices de rigidez
  
   Kel = zeros(6, 6, Nel);
 for e = 1:Nel % matriz de elementos
     [x1, y1, x2, y2, l] = locales(e, T1, x); % función locales
     s = (y2 - y1) / l; % Sinus
     c = (x2 - x1) / l; % Cosinus
     K1 = ((A(e) * E) / l) * [1 0 0 -1 0 0;...
                            0 0 0 0 0 0;...
                            0 0 0 0 0 0;...
                            -1 0 0 1 0 0;...
                            0 0 0 0 0 0;...
                            0 0 0 0 0 0]; 
     K2 = ((E * Ic(e)) / l^3) * [0 0 0 0 0 0;...
                            0 12 6*l 0 -12 6*l;...
                            0 6*l 4*l^2 0 -6*l 2*l^2;...
                            0 0 0 0 0 0;...
                            0 -12 -6*l 0 12 -6*l;...
                            0 6*l 2*l^2 0 -6*l 4*l^2];
     R = [c s 0 0 0 0;...
          -s c 0 0 0 0;...
          0 0 1 0 0 0;...
          0 0 0 c s 0;...
          0 0 0 -s c 0;...
          0 0 0 0 0 1];
      Kprima = K1 + K2; 
      K = R' * Kprima * R;
      for r = 1:(Nnod * Nig) % proceso de almacenamiento
         for s = 1:(Nnod * Nig)
             Kel(r, s, e) = K(r, s);
         end
         
     end
 end
 
 %% Matriz de rigidez global
 
KG = sparse(Ngl, Ngl);
for e = 1:Nel
    for i = 1:(Nnod * Nig)
        I = T2(i, e);
        for j = 1:(Nnod * Nig)
            J = T2(j, e);
            KG(I, J) = KG(I, J) + Kel(i, j, e); % Matriz de rigidez en globales
        end
    end
end
 
%% Sistema de ecuaciones globales
 
vL = [1 2 3 4 5 6 7 8 9]; % Grados libres 
vR = [10 11 12]; % Grados restringidos 
KLL = KG(vL, vL);
KLR = KG(vL, vR);
KRL = KG(vR, vL);
KRR = KG(vR, vR);
Fext = [0; 0; 0; 0; -P; 0; 0; 0; 0]; % Fuerzas exteriores
uR = [0; 0; 0]; % Desplazamientos restringidos 
 
% Resolución del sistema de ecuaciones
 
uL = KLL \ (Fext - KLR * uR); % Cálculo de desplazamientos libres
R = KRR * uR + KRL * uL; % Cálculo de reacciones
u(vL) = uL; 
u(vR) = uR; % Vector desplazamientos
 
%% Cálculo de los momentos y esfuerzos cortantes
 
% Barra 1
 
[x1, y1, x2, y2, l] = locales(1, T1, x);
s = norm((y2 - y1) / l); % Sinus
c = norm((x2 - x1) / l); % Cosinus
syms a b c d r
U = a*r^3 + b*r^2 + c*r + d; % Ecuación de desplazamientos
U1 = diff(U, r); % Ecuación para los giros
U0 = [subs(U, r, 0) == u(1)*s + u(2)*c; subs(U, r, l) == u(3)*s + u(5)*c];
U01 = [subs(U1, r, 0) == u(3); subs(U1, r, l) == u(6)]; 
Uf = [U0; U01]; % Vector de 4 ecuaciones con 4 incógnitas
Res = solve(Uf, a, b , c, d); % Parámetros ecuación cúbica
Me1 = [E * Ic(1) * 2*Res.b; E * Ic(1) * (6*Res.a*l + 2*Res.b)]; % Momento primera barra
Se1 = [-E * Ic(1) * 6*Res.a; -E * Ic(1) * 6*Res.a]; % Esfuerzo cortante primera barra
figure
l1 = linspace(0, l, 10);
Se1t = zeros(1, 10);
Se1t (1,:) = Se1(2);
subplot(1, 2, 1);
plot(l1, Se1t, 'b');
title 'Cortantes Barra 1';
xlabel 'x [m]';
ylabel 'Se [N]';
Me1t = zeros (10, 1);
for i = 1:10
Me1t(i)= Me1(2) + Se1(2)*l1(i);
end
subplot(1, 2, 2);
plot(l1, Me1t, 'r')
title 'Momentos Barra 1';
xlabel 'x [m]';
ylabel 'Me [N*m]';
 
% Barra 2
 
[x1, y1, x2, y2, l] = locales(2, T1, x);
syms a b c d r
U = a*r^3 + b*r^2 + c*r + d; % Ecuación de desplazamientos
U1 = diff(U, r); % Ecuación para los giros
U0 = [subs(U, r, 0) == u(1); subs(U, r, l) == u(6)];
U01 = [subs(U1, r, 0) == u(3); subs(U1, r, l) == u(9)]; 
Uf = [U0; U01]; % Vector de 4 ecuaciones con 4 incógnitas
Res = solve(Uf, a, b , c, d); % Parámetros ecuación cúbica
Me2 = [E * Ic(2) * 2*Res.b; E * Ic(2) * (6*Res.a*l + 2*Res.b)]; % Momento seguna barra
Se2 = [-E * Ic(2) * 6*Res.a; -E * Ic(2) * 6*Res.a];  % Esfuerzo cortante seguna barra
figure
l2 = linspace(0, l, 10);
Se2t = zeros(1, 10);
Se2t (1,:) = Se2(2);
subplot(1, 2, 1);
plot(l2, Se2t, 'b');
title 'Cortantes Barra 2';
xlabel 'x [m]';
ylabel 'Se [N]';
Me2t = zeros(10, 1);
for i = 1:10
Me2t(i) = Me2(2) + Se2(2)*l2(i);
end
subplot(1, 2, 2);
plot(l2, Me2t, 'r');
title 'Momentos Barra 2';
xlabel 'x [m]';
ylabel 'Me [N*m]';
 
% Barra 3
 
[x1, y1, x2, y2, l] = locales(3, T1, x);
syms a b c d r
U = a*r^3 + b*r^2 + c*r + d; % Ecuación de desplazamientos
U1 = diff(U, r); % Ecuación para los giros
U0 = [subs(U, r, 0) == u(8); subs(U, r, l) == u(5)];
U01 = [subs(U1, r, 0) == u(9); subs(U1, r, l) == u(6)]; 
Uf = [U0; U01]; % Vector de 4 ecuaciones con 4 incógnitas
Res = solve(Uf, a, b , c, d); % Parámetros ecuación cúbica
Me3 = [E * Ic(3) * 2*Res.b; E * Ic(3) * (6*Res.a*l + 2*Res.b)];  % Momento tercera barra
Se3 = [-E * Ic(3) * 6*Res.a -E * Ic(3) * 6*Res.a]; % Esfuerzo cortante tercera barra
figure
l3 = linspace(0, l, 10);
Se3t = zeros(1, 10);
Se3t (1,:) = Se3(2);
subplot(1, 2, 1);
plot(l3, Se3t, 'b');
title 'Cortantes Barra 3';
xlabel 'x [m]';
ylabel 'Se [N]';
Me3t = zeros(10, 1);
for i = 1:10
Me3t(i) = Me3(2) + Se3(2)*l3(i);
end
subplot(1, 2, 2);
plot(l3, Me3t, 'r')
title 'Momentos Barra 3';
xlabel 'x [m]';
ylabel 'Me [N*m]';
 
% Barra 4
 
[x1, y1, x2, y2, l] = locales(4, T1, x);
syms a b c d r
U = a*r^3 + b*r^2 + c*r + d; % Ecuación de desplazamientos
U1 = diff(U, r); % Ecuación para los giros
U0 = [subs(U, r, 0) == u(7); subs(U, r, l) == u(11)];
U01 = [subs(U1, r, 0) == u(9); subs(U1, r, l) == u(12)]; 
Uf = [U0; U01]; % Vector de 4 ecuaciones con 4 incógnitas
Res = solve(Uf, a, b , c, d); % Parámetros ecuación cúbica
Me4 = [E * Ic(4) * 2*Res.b; R(3)]; % Momento cuarta barra
Se4 = [-E * Ic(4) * 6*Res.a; R(1)]; % Esfuerzo cortante cuarta barra
figure
l4 = linspace(0, l, 10);
Se4t = zeros(1, 10);
Se4t(1,:) = Se4(2);
subplot(1, 2, 1);
plot(l4, Se4t, 'b');
title 'Cortantes Barra 4';
xlabel 'x [m]';
ylabel 'Se [N]';
Me4t = zeros(10, 1);
for i = 1:10
Me4t(i) = Me4(2) + Se4(2)*l4(i);
end
subplot(1, 2, 2);
plot(l4, Me4t, 'r')
title 'Momentos Barra 4';
xlabel 'x [m]';
ylabel 'Me [N*m]';
 
 
 

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).