1 |
- function [U,R] = statique
% Calcul, assemblage et resolution d'un probleme statique KU = F
% Les conditions aux limites en deplacement sont prises en compte par la
% methode du terme unite sur la diagonale.
%================================================================
% appel : [U,R] = statique
% en sortie valeurs nodale du champ inconnu (deplacements)
% et des flux inconnus (reactions)
%
% H.Oudin
% version MATLAB 7.4.0
%================================================================
global nddln nddlt nelt nnode ncld
global Connec Typel Ncl Vcl F
t0=cputime; %----- initialisation du temps de calcul
K=zeros(nddlt); %----- initialisation de la matrice rigidite
Fg=zeros(nddlt,1); %----- initialisation du vecteur force generalisee
Kr=zeros(ncld,nddlt);R=zeros(ncld,1);
for iel=1:nelt %----- boucle sur les elements
[Ke,Fe] = feval(Typel(iel,:),iel);
loce=[]; %---- localisation des ddl de l'element
for i=1:nnode
if Connec(iel,i) > 0
loce=[loce,(Connec(iel,i)-1)*nddln+[1:nddln]];
end
end;
K(loce,loce)=K(loce,loce) + Ke; %----- assemblage
Fg(loce)=Fg(loce) + Fe;
end
F = F + Fg; %----- vecteur du chargement
ir=0; %----- prise en compte des conditions aux limites
for i=1:nddlt % technique du terme unite sur la diagonale
if ( Ncl(i) == 1 )
ir=ir+1; R(ir,1) = -F(i);
F = F - K(:,i)*Vcl(i);
Kr(ir,:) = K(i,:);
K(i,:)=zeros(1,nddlt); K(:,i)=zeros(nddlt,1); K(i,i)=1;
F(i) = Vcl(i);
end
end
%disp(K);disp(F);
U = K\ F; %----- resolution
%disp(Sol);
R(:,1) = R(:,1) + Kr*U; %----- calcul des reactions (ddl imposes)
ir=ncld; %----- on etend a l'ensemble des ddl
for i=nddlt:-1:1
if ( Ncl(i) == 0 ) R(i,1) = 0;
else R(i,1)= R(ir,1);ir=ir-1;
end
end
|