TMCu, AnhB - FEM-Welle, RitzMATLAB

Anhang B: Berechnung einer konischen Welle, Verifizieren der Ergebnisse mit dem RITZschen Verfahren mit MATLAB

Ein MATLAB-Beispiel
wird hier ausführlich
beschrieben

Die auf Seite 712 angegebenen Ergebnisse für die Berechnung der Durchbiegung der stückweise konischen Welle (Elastizitätsmodul: E = 210000 N/mm²) nach der Finite-Elemente-Methode sollen mit dem RITZschen Verfahren verifiziert werden.

Die nachfolgend gelistete M-Datei lässt die Berechnung mit maximal 5 Polynom-Ansätzen zu, die jeweils für die gesamte Welle gelten:

v = a₁ [- (z/l)²+ z/l] + a₂ [-(z/l)³+ (z/l)²] + a₃ [(z/l)⁴- 2(z/l)³+ (z/l)²]
+ a₄ [3(z/l)⁵- 7(z/l)⁴+ 4(z/l)³] + a₅ [4(z/l)⁶- 9(z/l)⁵+ 5(z/l)⁴] .

Jede Funktion erfüllt für sich die geometrischen Randbedingungen (l ist die Gesamtlänge der Welle), die beiden letzten Funktionen auch die Bedingung der Momentenfreiheit an den Rändern.

% Ritzsches Verfahren für Biegeträger mit Polynom-Ansätzfunktionen,
% hier: Konische Welle auf Seite 711

clear all

% Parameter:
tl1 = 100 ;
tl2 = 250 ;
tl3 = 250 ;
tl4 = 400 ;
tl5 = 100 ;
tl = tl1 + tl2 + tl3 + tl4 + tl5 ;
E = 210000 ;
d1 = 20 ;
d2 = 50 ;
F1 = 2000 ;
F2 = 4000 ;

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

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 = 110 ; % 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

for i=1:nS % Durchmesser
   if (zS(i) < tl1) d(i) = d1 ;
   elseif (zS(i) < tl1+tl2) d(i) = d1 + (d2-d1)/tl2 * (zS(i)-tl1) ;
   elseif (zS(i) < tl1+tl2+tl3) d(i) = d2 ;
   elseif (zS(i) < tl1+tl2+tl3+tl4) d(i) = d2 - (d2-d1)/tl4 * (zS(i)-tl1-tl2-tl3) ;
   else d(i) = d1 ;
   end
end

EI = E * pi*d.^4/64 ; % Biegesteifigkeiten

% 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)' ;

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

% Punkte der beiden Einzellasten:
pos2 = n * (tl1+tl2)/tl + 1 ;
pos3 = n * (tl1+tl2+tl3)/tl + 1 ;

for ii = 1:m % Schleife über alle Gleichungen
   b(ii) = F1*vS(pos2,ii) + F2*vS(pos3,ii) ; % Rechte Seite (Einzellastanteile)
for jj = 1:m % ii-te Zeile (Koeffizienten A)
   Summe = EI(1)*vddS(1,ii)*vddS(1,jj) + EI(nS)*vddS(nS,ii)*vddS(nS,jj) ;
   faktor = 4 ; % Numerische Integration ...
   for k = 2:n % ... nach Simpsonscher Regel ...
   Summe = Summe + EI(k) * vddS(k,ii) * vddS (k,jj) * faktor ;
   if (faktor == 4) faktor = 2 ; % ... für Anteil aus der ...
   else faktor = 4 ; % ... Biegesteifigkeit ...
   end ;
   end
   A(ii,jj) = Summe*dz/3 ;
end
end

% Lösen des Gleichungssystems
ai = A\b ;

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

% Spezielle Werte:
DurchbiegungLinks = vSchlange (1)
DurchbiegungRechts = vSchlange (nS)
MaximaleDurchbiegung = max(abs(vSchlange))
pos1 = n * (tl1)/tl + 1 ;
pos4 = n * (tl1+tl2+tl3+tl4)/tl + 1 ;
vPos1 = vSchlange(pos1)
vPos2 = vSchlange(pos2)
vPos3 = vSchlange(pos3)
vPos4 = vSchlange(pos4)

Die am Ende der Datei farblich hervorgehobenen Zeilen sorgen für die graphische Ausgabe der Biegelinie (Bild unten) und die Ausgabe der Durchbiegung an den markanten Punkten (in das “Command Window”, Bild rechts), die auch auf Seite 712 angegeben sind (jeweils die Stellen des Übergangs vom konstanten Querschnitt auf einen veränderlichen Querschnitt bzw. umgekehrt).

Die Ergebnisse zeigen, dass die Qualität sich trotz der 5 Ansatzfunktionen in Grenzen hält, als Verifikation der mit der FEM berechneten Ergebnisse sind sie allerdings durchaus brauchbar.

Die oben angegebene M-Datei kann durch Modifikation der am Anfang farblich gekennzeichneten Zeile (m = 5) die Berechnung mit weniger Ansatzfunktionen ausführen. Die nachfolgende Tabelle zeigt, dass sich (wie zu erwarten) die Ergebnisse, die von “unbrauchbar” (eine oder zwei Ansatzfunktionen) bis zu “Abweichungen um etwa 10%” (fünf Ansatzfunktionen) reichen, mit zunehmender Anzahl der Ansatzfunktionen verbessert:

	RITZsches Verfahren					FEM
	m = 1	m = 2	m = 3	m = 4	m = 5	35 Elemente
vPos1 [mm]	1,4250	1,2455	2,2200	2,3731	2,3397	2,8192
vPos2 [mm]	3,7406	3,6351	4,9494	5,0325	5,0470	5,5267
vPos3 [mm]	4,2750	4,5723	5,8808	5,9336	5,9577	6,4245
vPos4 [mm]	1,4250	1,7469	2,9314	2,8067	2,7982	3,1986

Die oben gelistete M-Datei ist als RitzBiegungKonWell.m zum Download verfügbar.

Homepage

www.D@nkert.de

nkert.de