Zur Startseite

Technische Mechanik mit Computer

MATLAB-Script “RITZ-Verfahren mit Polynomansatzfunktionen”

Das nachfolgend gelistete MATLAB-Script realisiert das Näherungsverfahren von RITZ für gerade Biegeträger, die durch Einzelkräfte und -momente und Linienlasten belastet und zusätzlich zu starren Lagern auch in Federn und Drehfedern gelagert sein können. Für das elastische Potential

wird mit Ansatzfunktionen vi , die jede für sich die gemetrischen Randbedingungen erfüllen müssen, ein Näherungsansatz

mit den Freiwerten ai verwendet. Die ai werden aus den Minimalbedingungen

bestimmt, die auf ein lineares Gleichungssystem (m Gleichungen mit m Unbekannten) führen. Bei 3 Ansatzfunktionen kann es in der Form

aufgeschrieben werden, wobei die Elemente der Koeffizientenmatrix nach

und die Elemente des Vektors der rechten Seite nach

zu berechnen sind.

Diese Formeln gelten natürlich für beliebige Anzahl von Ansatzfunktionen.

Diesen Algorithmus (einschließlich Lösung des Gleichungssystems und Auswertung des Näherungsansatzes mit den berechneten Ansatzparametern) realisiert das nachfolgend gelistete Script. Dabei werden die Integrale numerisch nach der Simpson-Regel (mit beliebig fein einstellbarer Anzahl der Integrationsabschnitte) gelöst.

Das Script ist für maximal 5 Ansatzfunktionen eingerichtet. Wenn - wie im Script vorgesehen - ausschließlich Polynome verwendet werden, können diese mit den MATLAB-Funktionen polyder abgeleitet und Funktionswerte mit polyval berechnet werden. Bei Verwendung anderer Ansatzfunktionen müssen diese Berechnung entsprechend geändert werden.

Die hellblauen Zeilen, die an ein aktuelles Problem angepasst werden müssen, sind in diesem Beispiel-Script ausgefüllt für die Aufgabe, die hier beschrieben wird.

% Ritzsches Verfahren für Biegetraeger mit Polynom-Ansatzfunktionen,
% Musterprogramm: Es koennen bis zu m=5 Ansatzfunktionen definiert werden,
% integriert wird nach der Simpson-Regel mit n Abschnitten (n kann beliebig
% - aber geradzahlig - gewaehlt werden).
%
% Wenn man - wie unten empfohlen - ein Feld mit den Werten fuer die Biegesteifigkeiten
% aller Stuetzstellen und ein weiteres Feld mit den Linienlastintensitaeten an den
% Stuetzstellen aufbaut, dann muss man nur die vier mit den Kommentarzeilen
% “ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... ... BIS HIER”
% eingerahmten Bereiche an das aktuelle Problem anpassen.
%
% In dieser Datei sind in diesen Bereichen die Werte fuer die Aufgabe 33-9 eingetragen.

clear all

% 11111111111 ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... 111111111111111111111111111
% Parameter:
tl = 1000 ;
E  = 210000 ;
I0 = 200000 ;
q1 = 1 ;
F  = 500 ;
cT = 100000000 ;

% Ansatzfunktionen und deren 1. und 2. Ableitungen:
m  = 5 ;                       % Anzahl der Ansatzfuntionen (wenn m < 5 gesetzt ...
P1 = [1/tl   0] ;              % ... wird, werden nur die ersten ...
P2 = [1/tl^2 0 0] ;            % ... m Ansatzfunktionen verwendet)
P3 = [1/tl^3 0 0 0] ;
P4 = [1/tl^4 0 0 0 0] ;
P5 = [1/tl^5 0 0 0 0 0] ;
% 11111111111 ... BIS HIER 111111111111111111111111111111111111111111111111111111111111

P1D = polyder (P1);
P2D = polyder (P2);
P3D = polyder (P3);
P4D = polyder (P4);
P5D = polyder (P5);

P1DD = polyder (P1D) ;
P2DD = polyder (P2D) ;
P3DD = polyder (P3D) ;
P4DD = polyder (P4D) ;
P5DD = polyder (P5D) ;

% ------------------------------------------------------------------------------------
% Vorbereitung der numerischen Integration und der Ergebnisauswertung:
n  = 100 ;               % Anzahl der Abschnitte für numerische Integration
dz = tl / n ;            % Breite eines Integrationsintervalls
nS = n + 1 ;             % Anzahl der Stützstellen
zS = 0 : dz : tl ;       % Koordinaten der Stützpunkte

EIS = zeros (nS , 1) ;
qS = zeros (nS , 1) ;

% 22222222222222 ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... 222222222222222222222222
for i = 1:nS
   EIS(i) = E * I0 ;                  % Biegesteifigkeit und ...
   qS(i)  = q1 * (1 - zS(i)/tl) ;    % ... Linienlast
end
% 22222222222222 ... BIS HIER 222222222222222222222222222222222222222222222222222222222

% Funktionswerte und 1. und 2. Ableitung der Ansatzfunktionen an den Stützdtellen
vS (:,1) = polyval (P1  , zS)' ;
vS (:,2) = polyval (P2  , zS)' ;
vS (:,3) = polyval (P3  , zS)' ;
vS (:,4) = polyval (P4  , zS)' ;
vS (:,5) = polyval (P5  , zS)' ;
vdS (:,1) = polyval (P1D , zS)' ;
vdS (:,2) = polyval (P2D , zS)' ;
vdS (:,3) = polyval (P3D , zS)' ;
vdS (:,4) = polyval (P4D , zS)' ;
vdS (:,5) = polyval (P5D , zS)' ;
vddS(:,1) = polyval (P1DD , zS)' ;
vddS(:,2) = polyval (P2DD , zS)' ;
vddS(:,3) = polyval (P3DD , zS)' ;
vddS(:,4) = polyval (P4DD , zS)' ;
vddS(:,5) = polyval (P5DD , zS)' ;

subplot(311)
plot(zS , vS(:,1:m)) , axis ij
title('Ansatzfunktionen') ;

% ------------------------------------------------------------------------------------
% Aufbau des Gleichungssystems zur Bestimmung der Ansatzparameter:
K  = zeros (m , m) ;
b  = zeros (m , 1) ;

for ii = 1:m                                         % Schleife über alle Gleichungen
   Summe = (qS(1)*vS(1,ii)+qS(nS)*vS(nS,ii)) ;      % Rechte Seite ...
   faktor = 4 ;
   for k = 2:n                                      % Numerische Integration ...
     Summe = Summe + qS(k)*vS(k,ii)*faktor ;        % ... für Linienlastanteil ...
     if   (faktor == 4) faktor = 2 ;
     else               faktor = 4 ;
     end ;
   end
% 33333333333333 ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... 3333333333333333333333333
   b(ii) = Summe*dz/3 +
F*vS(nS,ii) ;               % ... + Einzelkraftanteil,
   % + weitere Einzelkraftanteile + Anteile von aeusseren Momenten
% 33333333333333 ... BIS HIER ... 333333333333333333333333333333333333333333333333333333
  for jj = 1:m                                       % ii-te Zeile (Koeffizienten A)
   Summe = EIS(1)*vddS(1,ii)*vddS(1,jj) + EIS(nS)*vddS(nS,ii)*vddS(nS,jj) ;
   faktor = 4 ;                                     % Numerische Integration ...
   for k = 2:n                                      % ... nach Simpsonscher Regel ...
     Summe = Summe + EIS(k) * vddS(k,ii) * vddS (k,jj) * faktor ;
     if   (faktor == 4) faktor = 2 ;                % ... für Anteil aus der ...
     else               faktor = 4 ;                % ... Biegesteifigkeit ...
     end ;
   end
% 44444444444444 ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... 4444444444444444444444444
  
K(ii,jj) = Summe*dz/3 + cT*vdS(1,ii)*vdS (1,jj) ; % ... + Drehfederanteil
   % + weitere Anteile von Federn
% 44444444444444 ... BIS HIER ... 444444444444444444444444444444444444444444444444444444
  end 
end

% Lösen des Gleichungssystems
ai = zeros (5,1) ;              % Maximal moegliche Zahl der Ansatzfznktionen: 5
aim = K\b ;
ai(1:m) = aim(1:m) ;

% Biegelinie graphisch darstellen:
vSchlange = vS * ai ;
subplot(312)
plot(zS , vSchlange) , axis ij
title('Durchbiegung') , ylabel('mm') ;

% Biegemoment graphisch darstellen:
v2Strich = vddS * ai ;
Mb = - EIS .* v2Strich;
subplot(313)
plot(zS , Mb) , title('Biegemoment')
ylabel('Nmm') ;

% Spezielle Werte:
iMitte               = round (n/2 + 1.1) ;
DurchbiegungMitte    = vSchlange (iMitte)
MaximaleDurchbiegung = max(abs(vSchlange))
vStrich = vdS * ai ;
BiegewinkelLinks     = vStrich (1)
BiegewinkelRechts    = vStrich (nS)
BiegemomentLinks     = Mb (1)
BiegemomentMitte     = Mb (iMitte)
BiegemomentRechts    = Mb (nS)
MaximalesBiegemoment = max(abs(Mb))

Das oben gelistete Script ist zum Download als RitzBiegungSimpson5.m verfügbar.

Homepage

www.D@nkert.de

D

nkert.de