Esempio di soluzione di una struttura reticolare

Questo script illustra la soluzione matriciale di una struttura reticolare, nel caso generale D-dimensionale. Il codice è valido per D >= 1; viene presentato un esempio per D == 2

Definizione della struttura

Costanti

D : numero delle dimensioni dello spazio fisico
nnodel : nodi per elemento

dimensioni della struttura:

N : numero di nodi
M : numero di aste

geometria e topologia

XY : matrice delle coordinate nodali, di dimensione (D,N) letta per colonne contiene le coordinate degli N nodi
INC : matrice delle incidenze nodali, di dimensione (nnodel, M) letta per colonne contiene il numero del primo e del secondo nodo
kel : vettore delle rigidezze elementari, di dimensione M

% problema piano 2D
D = 2;
% due nodi per elemento
nnodel = 2;

% struttura con quattro nodi
N = 4;

% struttura con cinque aste
M = 5;

% coordinate nodali (mm)
XY = [ 0.0, 10.0,  0.0, 10.0;
       0.0,  0.0, 10.0, 10.0 ];

% incidenze
INC = [1, 2, 3, 2, 4;
       2, 3, 1, 4, 3];

% rigidezze (N/mm)
kel = [10, 10, 10, 10, 10];

% verifico la coerenza delle dimensioni delle strutture dati
assert(all(size(XY) == [D, N]), 'matrice XY: errore dim');
assert(all(size(INC) == [nnodel, M]), 'matrice INC: errore dim');
assert(length(kel) == M, 'matrice kel: errore dim');

% verifico che in INC ci siano tutti e soli gli interi 1:N
assert(all(unique(INC) == (1:N)'), 'errore matrice INC')

Vincoli e carichi

vincoli

vinc: matrice dei vincoli, di dimensione (numero nodi vincolati, 3). Ogni riga rappresenta un vincolo semplice (carrello) e contiene tre valori: numero del nodo, indice del grado di libertà vincolato (x → 1, y → 2, …) , spostamento imposto (zero per vincoli perfetti)

carichi

forc: matrice delle forze nodali, di dimensione (numero di nodi caricati, 3). Ogni riga contiene una componente di forza assegnata secondo uno schema simile ai vincoli: numero del nodo, indice del grado di libertà, componente di forza.

% spostamenti imposti (mm)
vinc = [ 1, 1, 0.0;
         1, 2, 0.0;
         2, 2, 0.0 ];

% forze nodali (N)
forc = [ 4, 1, 100.
         4, 2, 100. ];

% controlli sulle dimensioni delle matrici
assert(size(vinc, 2) == 3);
assert(size(forc, 2) == 3);

% controlli sugli elementi delle matrici
assert(all(ismember(vinc(:, 1),1:N)))
assert(all(ismember(forc(:, 1),1:N)))
assert(all(ismember(vinc(:, 2),1:D)))
assert(all(ismember(forc(:, 2),1:D)))

output dati di input

disp '************'
disp 'Dati di input'
disp ' '

disp 'coordinate nodali:'
disp(XY)

disp 'incidenze e rigidezze'
disp(vertcat(INC, kel))

disp 'vincoli'
disp(vinc)

disp('forze nodali')
disp(forc)

************
Dati di input
 
coordinate nodali:
     0    10     0    10
     0     0    10    10

incidenze e rigidezze
     1     2     3     2     4
     2     3     1     4     3
    10    10    10    10    10

vincoli
     1     1     0
     1     2     0
     2     2     0

forze nodali
     4     1   100
     4     2   100

Matrici ausiliarie

% gradi di libertà per nodo
gdln = D;
% gradi di libertà totali
gdl = gdln*N;
% gradi di libertà per elemento
gdle = gdln*nnodel;

% mappa da (componente, nodo) -> gdl
maptogdl = reshape(1:gdl,D,N);

Assemblaggio della matrice di rigidezza globale

ciclo sugli elementi per assemblare la matrice di rigidezza globale

% preparo matrice di zeri delle dimensioni opportune
K = zeros(gdl,gdl);
eyegdle = eye(gdle, gdle);
% ciclo sugli elementi
for i=1:M
    % matrice di rigidezza dell'elemento
    ki = kbar(XY(:,INC(:,i)),kel(i));
    % creo la matrice di assemblaggio
    L = zeros(gdle,gdl);
    jj = maptogdl(:,INC(:,i));
    jj = jj(:);
    L(:, jj) = eyegdle;
    % assemblo la matrice di rigidezza
    K = K + L'*ki*L;
end

% matrice di rigidezza assemblata
% disp(K)

% verifiche finali matrice assemblata
% tol è una tolleranza relativa, definita rispetto all'epsilon macchina
tol = 1000 * eps;
% vettore che conterrà i moti rigidi del sistema
urig = zeros(gdl, D*(D+1)/2);
% calcolo traslazioni rigide
for i = 1:D
    % spostamento unitario in direzione i
    urig(maptogdl(i,:),i) = 1;
end
% calcolo rotazioni rigide
soD = rots(D);
for i = 1:size(soD, 3)
    rot = soD(:,:,i)*XY;
    urig(:,D+i) = rot(:);
end
assert(size(urig,2) == D+size(soD,3))
assert(rank(urig) == min(size(urig)))

% (K * urig) è sufficientemente piccolo?
assert(norm(K*urig,inf) < tol*norm(urig,inf)*norm(K,inf))
% è giusto il rango della matrice?
assert(rank(K) <= gdl - size(urig, 2))

Imposizione dei vincoli e dei carichi nodali

i vincoli vengono imposti mediante decomposizione a blocchi della matrice di rigidezza

% calcolo gdl vincolati
gdv = (vinc(:,1)-1)*gdln + vinc(:,2);
% uv sono gli spostamenti imposti
u_V = vinc(:, 3);

% verifico che non ci siano ripetizioni
assert(length(gdv) == length(unique(gdv)));

% calcolo gdl incogniti
gdi = setdiff(1:gdl, gdv);

% assemblo il vettore delle forze nodali
f = zeros(gdl, 1);
for i = 1:size(forc,1)
    j = maptogdl(forc(i, 2), forc(i, 1));
    f(j) = f(j) + forc(i, 3);
end

% procedo alla decomposizione a blocchi della matrice

K_II = K(gdi,gdi);
K_IV = K(gdi,gdv);
K_VV = K(gdv,gdv);
f_I = f(gdi);
f_V = f(gdv);

Calcolo di spostamenti e reazioni vincolari

si risolve il sistema ridotto e si calcolano le reazioni vincolari

% calcolo spostamenti nodali e reazioni vincolari
u_I = K_II\(f_I-K_IV*u_V);
r_V = K_IV'*u_I + K_VV*u_V - f_V;

% ricostruisco il vettore completo degli spostamenti
u = zeros(gdl,1);
u(gdi) = u_I;
u(gdv) = u_V;
r = zeros(gdl,1);
r(gdv) = r_V;

output dei risultati

disp '************'
disp 'Risultati'
disp ' '

disp 'spostamenti:'
disp(reshape(u,D,N))

disp 'forze nodali:'
disp(reshape(f,D,N))

disp 'reazioni vincolari:'
disp(reshape(r,D,N))

************
Risultati
 
spostamenti:
         0   10.0000   40.0000   50.0000
         0         0   10.0000   10.0000

forze nodali:
     0     0     0   100
     0     0     0   100

reazioni vincolari:
 -100.0000         0         0         0
 -100.0000   -0.0000         0         0

verifiche finali

verifico che le equazioni di equilibrio siano soddisfatte

assert(norm(f + r - K*u,inf) < tol*norm(f,inf))