Tribologie + Schmierungstechnik 63. Jahrgang 5/ 2016 Einführung Im Schmierfilm eines Gleitlagers kann es aufgrund dynamischer Prozesse zu einem derart schnellen Druckabfall unterhalb des Dampfdrucks kommen, dass das Schmiermittel mit seinem Übergang in die Dampfphase nicht direkt folgen kann. In diesem Zeitraum des Nichtgleichgewichts treten Zugspannungen im Lager auf, die allmählich durch das Blasenwachstum abgebaut werden. Das von Geike und Popov [3] entwickelte Kavitationsmodell für Schmierfilmkontakte basiert auf einem dynamischen Blasenmodell nach Sauer [6] und ist in der Lage, die auftretenden Zugspannungen während des dynamischen Prozesses zu berechnen. Aufgrund des rechenzeitintensiven Lösungsalgorithmus mittels Linienmethode (Method of Lines) ist das Kavitationsmodell bisher nur für einen radialsymmetrischen (1D) Testfall validiert worden. Im Rahmen dieser Arbeit wird eine Finite-Elemente-Diskretisierung des Modells vorgestellt und anhand des radialsymmetrischen Testfalls validiert. Anschließend wird der Algorithmus zur Berechnung eines zweidimensionalen dynamischen Gleitschuhtests angewendet und mit dem Kavitationsmodell von Kumar [4] verglichen. Physikalische Modellbildung Aus den Erhaltungsgleichungen für Masse, Impuls und Energie kann mit den in der Lagertechnik üblichen physikalischen Vereinfachungen die Reynoldssche Differentialgleichung (DGL) für die Druckverteilung p bei zeitlich variabler Spaltweite h in folgender Form hergeleitet werden [4] (1) Die kompressible Form der Gleichung ist hier gewählt worden, da die Dichte ρ und die Viskosität η vom Volumenanteil α des Dampfs in dem ansonsten flüssigen Schmiermittel abhängig sind. Für die Zweiphasenströmung gilt: 25 Aus Wissenschaft und Forschung * Lukas van Husen, M.Sc. Prof. Dr. rer. nat. Siegfried Müller Institut für Geometrie und Praktische Mathematik, RWTH Aachen, 52062 Aachen Dr.-Ing. Birgit Reinartz Prof. Dr.-Ing. habil.Gunter Knoll IST GmbH, Schloss-Rahe-Str. 12, 52072 Aachen Berechnung von Kavitationsvorgängen in hydrodynamischen Gleitlagern auf der Basis eines blasendynamischen Modells L. v. Husen, S. Müller, B. Reinartz, G. Knoll* Eingereicht: 25. 9. 2015 Nach Begutachtung angenommen: 30. 10. 2015 Ein Finite Elemente Algorithmus zur Berechnung von Kavitation in Gleitlagern wird vorgestellt. Der dynamische Kavitationsvorgang wird beschrieben durch die zweidimensionale Reynoldssche Differentialgleichung für den Druck im Schmierfilm-Dampf Gemisch. Die Reynoldssche Gleichung wird gekoppelt mit einer partiellen Differentialgleichung für den Dampfanteil, welche auf einem dynamischen Blasenmodell basiert. Das Verfahren wird anhand einer eindimensionalen Vergleichsrechnung validiert und für eine zweidimensionale Gleitschuhsimulation mit dem Modell nach Kumar verglichen. Schlüsselwörter Kavitation, Gleitlager, Reynolds Gleichung, Blasendynamik, FEM, Zugspannungen A Finite Element algorithm for the simulation of cavitation in lubrication flows is presented. The dynamic cavitation model is based on the two-dimensional Reynolds’ differential equation and coupled to a relation for vapor change due to bubble growth. The method is validated using a one-dimensional test case from literature and computed results of a two-dimensional lubrication flow are compared to the cavitation model by Kumar. Keywords Cavitation, lubrication flows, Reynolds equation, bubble growth, FEM, tension Kurzfassung Abstract ( ∇ · ( ρ(α) h 3 12 η(α) ∇p ) = ∇ · (ρ(α) hU) + ∂ ∂t (ρ(α) h) . T+S_5_16 29.07.16 11: 27 Seite 25 26 Tribologie + Schmierungstechnik 63. Jahrgang 5/ 2016 (2) Als Zwischenergebnis der Herleitung der Reynoldsschen Gleichung ergibt sich aus der Impulserhaltung die Gleichung für die über die Spalthöhe gemittelte Geschwindigkeit in folgender Form [4]: (3) wobei U sich aus der Relativgeschwindigkeit der Lagerseiten zueinander ergibt. Im Rahmen einer reduzierten Modellbildung werden die für den Kavitationsprozess benötigten Keime als Mikroblasen idealisiert, die einen mittleren Radius R 0 und eine homogene Verteilungsdichte n 0 aufweisen [5]. Dabei ist die Stoffkonstante n 0 definiert als Anzahl Kavitationskeime pro Volumen Flüssigkeit. Bei Unterschreiten des Kavitationsdrucks wachsen alle Blasen gleichförmig und unabhängig voneinander aus den Mikroblasen. Diese vereinfachte Modellannahme gilt nur zu Beginn des Blasenwachstums und ist für große Kavitationsgebiete nicht zulässig, da sich hier Dampfblasen verbinden [3]. Um ein dynamisches Modell für den Kavitationsdruck zu Beginn der Ablösung zu erhalten, wird nach Sauer [5] eine Transportgleichung für den Dampfanteil α formuliert, die sowohl die Änderung von α aufgrund von Konvektion als auch durch Blasenwachstum bzw. -kollaps berücksichtigt. Um die Ableitung von α in Strömungsrichtung zu beschreiben, wird die substantielle Ableitung von α mit der konvektiven Geschwindigkeit u aus (3) formuliert: (4) Die substantielle Zeitableitung kann durch die Änderung des Blasenradius R ausgedrückt werden. Ausgehend von der Definition von α als Dampfanteil und der Stoffkonstanten n 0 findet sich: (5) Die zeitliche Ableitung ergibt für die Änderung des Dampfanteils nur eine Abhängigkeit von der zeitlichen Änderung des Blasenradius: (6) Geike und Popov [2] folgend wird das Blasenwachstum durch die Rayleigh-Plesset Gleichung angenähert, die für die hier analysierte Initialphase des Blasenwachstums gültig ist [6]. Vereinfachend wird der Einfluss der Oberflächenspannung sowie die zeitliche Ableitung zweiter Ordnung vernachlässigt und nur die Wechselwirkung von Trägheit, Viskosität und Druck berücksichtigt. Aus der so vereinfachten quadratischen Rayleigh-Plesset Gleichung ergibt sich die Änderung des Blasenradius zu [7]: (7) wobei p v der Kavitationsdruck ist und der letzte Term die Annährung an den Kavitationsdruck und den Abbau der Zugspannungen aufgrund des Blasenwachstums ausdrückt. Anschließend werden die Gleichungen (3), (6) und (7) nacheinander in die Transportgleichung für den Dampfanteil (4) eingesetzt und somit eine zweite Differentialgleichung erzeugt, die gemeinsam mit der Reynoldsschen Gleichung (1) sowie Anfangs- und Randbedingungen für die zeitliche Entwicklung des Drucks im kavitierenden Lager gelöst werden kann: (8) Hierin ist η(α) durch (2) bestimmt und aus (5) erhält man Damit ergeben sich folgende Abhängigkeiten aus den beiden DGLen (1) und (8): Numerisches Verfahren Die beiden DGLen (1) und (8) werden vor der Diskretisierung noch in eine dimensionslose Form gebracht. Dabei werden als Referenzgrößen die charakteristische Länge L, die Referenzzeit T ref = η liq / p amb , der Umgebungsdruck p amb sowie die Stoffwerte der Flüssigkeit ρ liq und η liq herangezogen. Weitere dimensionslose Variablen sind: Für Gleichung (7) findet sich folgende dimensionslose Darstellung: Die zu diskretisierenden Gleichungen in dimensionsloser Form ergeben sich wie folgt: Aus Wissenschaft und Forschung ( η(α) = αη vap + (1 − α)η liq . ( ρ(α) = αρ vap + (1 − α)ρ liq ( u = − h 2 12 η(α) ∇p + U , ( dα dt = ∂α ∂t + u · ∇α . . ( α = V vap V = n 0 V liq 43 πR 3 V = n 0 (1 − α) 4 3 πR 3 . ( dα dt = 3α(1 − α) 1 R dR dt , ( dR dt = − 8η liq 6ρ liq R + √ 64η 2 liq 36ρ 2liq R 2 + 2(p v − p) 3ρ liq ∂α ∂t = 3α(1 − α) 1 R(α) [ − 8η liq 6ρ liq R(α) + √ 64η 2 liq 36ρ 2 liq R(α) 2 + 2(p v − p) 3ρ liq ] + h 2 12η(α) ∇p · ∇α − U · ∇α . . R(α) = ( 3 4π n 0 · α 1 − α ) 1/ 3 L p (α, p; h, U, η liq , η vap , ρ liq , ρ vap ) = 0 L α (α, p; h, U, η liq , η vap , ρ liq , p v , n 0 ) = 0 . ¯ u = u T ref L , ¯ n 0 = n 0 L 3 , φ : = 8η 2 liq 3ρ liq p amb L 2 . ψ(¯ p, α) : = 1 ¯ R(α) d ¯ R(α)) d¯ t . = − φ 2 ¯ R(α) 2 + 1 2 ¯ R(α) √ φ 2 ¯ R(α) 2 + ( ¯ p v − ¯ p)φ T+S_5_16 29.07.16 11: 27 Seite 26 Tribologie + Schmierungstechnik 63. Jahrgang 5/ 2016 (9a) (9b) Dieses Problem muss durch geeignete Anfangs- und Randbedingungen geschlossen werden: (10a) (10b) (10c) wobei ∂ Ω_ den Einströmrand beschreibt und Ω das Berechnungsgebiet. Das Problem (9), (10) ist ein zeitabhängiges Problem, bei dem in jedem Zeitschritt der Dampfanteil berechnet werden muss. Da bei jeder Auswertung der rechten Seite in (9a) das Druckfeld p- benötigt wird, muss dazu die DGL (9b) gelöst werden. Da dabei aber wiederum der Dampfanteil eingeht, ist das Problem nichtlinear gekoppelt. Zum Zeitpunkt t - = 0 wird ausgehend von der Anfangsbedingung (10a) für den Dampfanteil das Druckfeld durch Lösen von (9b) und (10c) bestimmt. Zur Durchführung eines Zeitschritts werden nachfolgend die Teilprobleme (9a), (10a, b) und (9b), (10c) alternierend gelöst. Dampfanteil: Für die Diskretisierung der zeitlichen Ableitung von α in (9a) wird ein implizites Crank-Nicolson Verfahren verwendet. Dabei wird die daraus resultierende zeitlich semi-diskrete nichtlineare Transportgleichung in jedem Zeitschritt mittels einer Picard-Iteration (Fixpunktiteration mit dem Index k) gelöst und die Druckverteilung während der Iteration konstant gehalten. Als Startwert für die Picard-Iteration wird im ersten Zeitschritt ein konstanter Dampfanteil α 0 nahe Null vorgegeben und danach wird jeweils die α-Verteilung des vorangegangenen Zeitschritts gewählt. Als Abbruchkriterium wird eine Toleranz für die Differenz von α in der L 2 - Norm vorgegeben: Druckfeld: Die nichtlineare Druckgleichung wird mit einem Newton-Verfahren gelöst. Als Startwert für das Newton-Verfahren bietet sich im ersten Zeitschritt der konstante Umgebungsdruck p amb an, während in den folgenden Zeitschritten jeweils die Funktion des vorherigen Zeitschritts als Startwert verwendet wird. Als Abbruchkriterium wird analog zur Picard-Iteration eine Toleranz Tol new für p- in der L 2 -Norm vorgegeben: Diese beiden Iterationsprozesse werden solange alternierend durchgeführt, bis die Differenz der Druckfelder zweier aufeinanderfolgender alternierender Iterationen in der L 2 -Norm unter eine vorgegebene Toleranz Tol ges fällt. Zur räumlichen Diskretisierung der DGLen für den Druck und den Dampfanteil wird die Finite Elemente Methode angewendet. Zu diesem Zweck werden die Gleichungen in ihre schwache Formulierung überführt. Als Ansatzraum wird der Raum der stückweise linearen Elemente verwendet und die Lösung als Linearkombination der Elementknotenwerte konstruiert. Eine detaillierte Darlegung des Algorithmus findet sich bei van Husen [7]. Ergebnisse Das instationäre Kavitationsmodell soll nun anhand von zwei verschiedenen Konfigurationen getestet werden. Zunächst wird eine radialsymmetrische Konfiguration mit einer negativen Pressbewegung, wie sie bei Geike und Popov [3] Verwendung findet, betrachtet. Aufgrund der Symmetrie der Konfiguration können die Gleichungen unter Verwendung von Polarkoordinaten zu einem eindimensionalen Modell vereinfacht werden. Danach wird die Druckverteilung im Schmierspalt eines zweidimensionalen Gleitschuhs betrachtet. Radialsymmetrische Negative Pressung Es werden zwei Platten betrachtet, die durch einen Schmierfilm der Höhe h getrennt sind. Die untere Platte ist statisch, während die obere kreisrunde Platte (Radius L) mit konstanter Geschwindigkeit v t in Normalenrichtung bewegt wird, siehe Bild 1. 27 Aus Wissenschaft und Forschung ∂α ∂¯ t = ( 3α(1 − α)ψ(¯ p, α) + ¯ h 2 12 ¯ η(α) ∇¯ p · ∇α − ¯ U · ∇α ) ∇ · ( (α( ¯ ρ vap − 1) + 1)¯ h 3 12 ¯ η(α) ∇¯ p ) = ∇ · ( (α( ¯ ρ vap − 1) + 1)¯ h ¯ U ) ( ¯ h 2 ( ( ) ) +(α( ¯ ρ vap −1)+1) ¯ V +¯ h( ¯ ρ vap −1) ( 3α(1 − α)ψ(¯ p, α) + ¯ h 2 12 ¯ η(α) ∇¯ p · ∇α − ¯ U · ∇α ) . α(0, ¯ x) = α 0 (¯ x) , ¯ x ∈ ¯ Ω α(¯ t, ¯ x) = α R (¯ t, ¯ x) , ¯ x ∈ ∂ ¯ Ω − = {¯ x ∈ ∂ ¯ Ω − : ∇α(¯ t, ¯ x) · n(¯ x) < 0} ¯ p(¯ t, ¯ x) = p R (¯ t, ¯ x) , ¯ x ∈ ¯ Ω , . ‖α n+1 k+1 − α n+1 k ‖ L 2 (Ω) < T ol pic . ‖ ¯ p n+1 j+1 − ¯ p n+1 j ‖ L 2 (Ω) < T ol new Bild 1: Aufbau der radialsymmetrischen Konfiguration, entnommen aus [3] . T+S_5_16 29.07.16 11: 27 Seite 27 28 Tribologie + Schmierungstechnik 63. Jahrgang 5/ 2016 Aufgrund der Aufwärtsbewegung entsteht ein Unterdruck im Schmierfilm und die Kavitationskeime beginnen mit der Zeit zu wachsen. Die eindimensionale Form der Gleichungen in Polarkoordinaten ergibt: wobei das Berechnungsgebiet gewählt ist als Ω- = [0,1]. Es werden folgende Rand- und Anfangsbedingungen gestellt: und die Gleichungen mit folgenden Parametern entsprechend [3] gelöst: Ausgangspunkt sind die Ergebnisse von Geike und Popov [3], die mittels Finiter Differenzen (FD) berechnet wurden. Diese wurden hier reproduziert, diesmal mit einer äquidistanten Einteilung der Stützstellen anstatt eines Gauss-Lobatto-Gitters. Dabei wurden 50 Stützstellen verwendet und die Parameter für den ODE-Löser von Matlab aus [3] übernommen. Bei der hier vorgestellten neuen Methode wird die räumliche Diskretisierung mit der Finite Elemente Methode (FEM) durchgeführt. Es wurden 64 Gitterpunkte verwendet. Die resultierenden linearen Gleichungssysteme wurden dabei mit einem direkten Gleichungssystemlöser für dünn besetzte Matrizen basierend auf einer LU-Zerlegung gelöst. Bild 2 zeigt die zeitliche Entwicklung des radialen Druckverlaufs. In Bild 2b ist jeder zehnte Zeitschritt dargestellt, beginnend mit dem grünen Graphen zum Zeitpunkt t- = 0 bis zum roten Graphen für t- = 0.1. Man erkennt einen ähnlichen zeitlichen Verlauf wie in Bild 2a für das Finite Differenzen Verfahren nach Geike und Popov dargestellt. Abgesehen vom Randbereich befindet sich der Druck im gesamten Berechnungsgebiet unterhalb des Dampfdrucks, wodurch das Blasenwachstum angeregt wird. In Bild 3 ist der Dampfanteil α zum Endzeitpunkt t-= 0.1 dargestellt. Man erkennt, dass entsprechend der Druckverteilung in Bild 2 der Blasenanteil insbesondere zur Plattenmitte hin stark angestiegen ist. Um die beiden unterschiedlichen Lösungsverfahren zu vergleichen, ist in Bild 3b Aus Wissenschaft und Forschung ¯ h(0) = 10 −2 , ¯ η liq = 1 und ∆¯ t = 0.001 ¯ ρ vap = 0, ¯ p v = 0, α 0 = 10 −8 , ¯ V = 10 −4 , ¯ n 0 = 10 9 , φ = 3 · 10 −8 , ∂α ∂¯ t = 3α(1 − α)ψ(¯ p, α) + ¯ h 2 12¯ η(α) ∂ ¯ p ∂ ¯ r ∂α ∂ ¯ r in ¯ Ω und ( ) 1 ¯ r ∂ ∂ ¯ r ( ¯ r ¯ ρ(α)¯ h 3 12¯ η(α) ∂ ¯ p ∂ ¯ r ) = ¯ h ∂ ¯ ρ(α) ∂¯ t + ¯ ρ(α) ¯ V t , in ¯ Ω f¨ ur ¯ t ∈ (0, ¯ T end ] ¯ p(¯ t, 1) = 1 f¨ ur ¯ t ∈ [0, ¯ T end ] α(¯ t = 0, ¯ r) = α 0 ∀¯ r ∈ Ω α(¯ t, 1) = α 0 , falls − ∂α(¯ t, 1) ∂ ¯ r < 0 Bild 2: Vergleich des zeitlichen Verlaufs der eindimensionalen Druckverteilung im Zeitinterval [0; 0.1], links: FD analog zu [3], rechts: FEM Bild 3: Eindimensionaler Dampfanteil zum Zeitpunkt t- = 0.1; links: radialer Verlauf für FD und FEM, rechts: Differenz in α zwischen FD- und FEM-Lösung a b a b T+S_5_16 29.07.16 11: 27 Seite 28 Tribologie + Schmierungstechnik 63. Jahrgang 5/ 2016 die Differenz Δ α der beiden Lösungen aufgetragen. Dabei ist die FEM-Lösung an den Gitterpunkten der FD-Lösung ausgewertet. Zusammenfassend lässt sich feststellen, dass sich eine gute Übereinstimmung zwischen dem neu entwickelten Finite Elemente Verfahren und der validierten Lösung aus [3] ergibt. Zweidimensionale Gleitschuhlagerung In der zweidimensionalen Konfiguration, siehe Bild 4, werden zwei quadratische Platten der Länge L betrachtet, wobei die untere sich mit konstanter Geschwindigkeit in x-Richtung bewegt und die obere unbeweglich ist. Die Höhe des Schmierspalts ist nur abhängig vom Ort x und nicht mehr von der Zeit. Es wird ein stationäres parabelförmiges Spaltprofil angesetzt: Das dimensionslose Gleichungssystem (9) wird mit folgenden Anfangs- und Randbedingungen: sowie folgenden Parametern ergänzt: Für die Simulation ist das Berechnungsgebiet Ω = {(x,y); x ∈ [0,L],y ∈ [0,L]} in 4096 quadratische Zellen unterteilt worden. Die aus der FEM-Diskretisierung resultierenden linearen Gleichungssysteme werden mit einem iterativen Löser [1] basierend auf einem stabilisierten Gradientenverfahren biorthogonaler Richtungen (BiCGStab) mit Vorkonditionierung (SSOR) gelöst. Aufgrund des konvergenten Spaltverlaufs im Modell aus Bild 4 steigt der Druck erst an und fällt danach im divergenten Teil unter den Kavitationsdruck p v - = 0 (vgl. Bild 6). Im Unterschied zum eindimensionalen Modell gibt es Bereiche, in denen der dimensionslose Druck so groß ist, dass aufgrund des Modells der Dampfanteil α stark fällt. Physikalisch ist nur ein Dampfanteil α > 0 sinnvoll. Um dies sicherzustellen, wird die Ungleichung 29 Aus Wissenschaft und Forschung z A P h(x) = 4 h 1 − h 0 L 2 x 2 + 4 h 0 − h 1 L x + h 1 Bild 4: Parabelförmiges Gleitschuhprofil über bewegter, quadratischer Patte 1e-5 2e-5 3e-5 4e-5 1.000e-14 4.893e-05 20 40 60 -3.283e-01 7.158e+01 20 40 60 -3.283e-01 7.158e+01 1e-5 2e-5 3e-5 4e-5 1.000e-14 4.893e-05 Bild 5: Druckverteilung (unten) und Dampfanteil (oben) innerhalb des Spalts zum Zeitpunkt t- =10 Bild 6: Druckverteilung entlang der Mittellinie (y- =1⁄2) für verschiedene Zeitpunkte und im Vergleich zum Modell von Kumar [2] A P ¯ p(¯ t, ¯ x, ¯ y) = 1 f¨ ur (¯ x, ¯ y) ∈ ∂ ¯ Ω, ¯ t ∈ [0, ¯ T end ] P α(0, ¯ x, ¯ y) = α 0 f¨ ur (¯ x, ¯ y) ∈ ¯ Ω A P α(¯ t, ¯ x, ¯ y) = α 0 , falls (¯ x, ¯ y) ∈ ∂ ¯ Ω − , ¯ t ∈ [0, ¯ T end ], P mit: ∂ ¯ Ω − : = {(¯ x, ¯ y) ∈ ∂ ¯ Ω : −∇α(¯ t, ¯ x, ¯ y) · n < 0} A P α 0 = 10 −8 , ¯ U x = 2, 325 · 10 −6 , A P ¯ U y = 0, ¯ n 0 = 10 12 , φ = 3 · 10 −10 , ¯ ρ vap = 0, P ¯ p v = 0, ¯ h 0 = 10 −4 , ¯ h 1 = 2 · 10 −4 , A P p ¯ η liq = 1 und ∆¯ t = 0.01 . T+S_5_16 29.07.16 11: 27 Seite 29 30 Tribologie + Schmierungstechnik 63. Jahrgang 5/ 2016 nach jeder Picard-Iteration für alle Freiheitsgrade α i überprüft. Falls sie nicht erfüllt ist, wird lokal α > 0 gesetzt. Bild 5 zeigt die Verteilung des dimensionslosen Drucks und Dampfanteils zum Zeitpunkt t - =10 im Spalt in der Länge und Tiefe. Da die Verteilung symmetrisch ist, werden beide Größen jeweils auf der unteren bzw. oberen Hälfte dargestellt. Bild 6 zeigt die dimensionslose Druckverteilung entlang y- = 1 ⁄ 2 zu verschiedenen Zeitpunkten. Man erkennt deutlich zum Startzeitpunkt die typische Lösung der Reynoldsschen Differentialgleichung für einen konstanten Dampfanteil α 0 und somit für eine konstante Dichte. Der Druck steigt im zeitlichen Verlauf an und strebt im Kavitationsgebiet (p- < 0) gegen 0. Zum Vergleich ist eine Lösung mit dem kommerziellen Programm FIRST [2] berechnet worden. FIRST verwendet einen Kavitationsalgorithmus, der auf Kumar [5] basiert. Die zeitliche Konvergenz der Lösung des Drucks gegen die Lösung aus FIRST ist deutlich zu erkennen. Für den Dampfanteil konnte jedoch keine zeitlich stationäre Lösung nachgewiesen werden. Mit fortlaufender Zeit steigt α im Kavitationsgebiet immer weiter an, während sich der Druck kaum noch ändert. Die gewählte Transportgleichung (8) für α ist nicht konservativ formuliert; durch den konvektiven Term wird ein vom Druck unabhängiger Quellterm für α generiert. Eine mögliche Verbesserung könnte sein, dass man für die Druckdifferenz zum Kavitationsdruck eine Toleranz vorgibt, ab der p- = p v - gesetzt wird. Damit wird das Blasenwachstum gestoppt, und es findet in diesem Bereich nur noch Transport statt. Anderseits ist die Formulierung des Kavitationsmodells bewusst nur für den Inertialbereich des Schmierfilmabriss formuliert worden. Hier müsste ggf. eine Erweiterung der Modellierung erfolgen. Der Algorithmus von Kumar [5] löst im Kavitationsgebiet die Reynoldssche Differentialgleichung (1) nach der Dichte anstatt des Drucks auf und erhält damit eine Volumenfraktion für die Dichte, welche für den vorliegenden Fall konvergiert und eine stationäre Lösung liefert. Dabei muss jedoch ständig das Kavitationsgebiet neu bestimmt werden, was rechenzeittechnisch aufwendig ist. Umso bemerkenswerter ist, dass das hier vorgestellte Modell einen zu Kumar ähnlichen Druckverlauf zeigt, ohne eine derartige Einteilung vorzunehmen. Ein weiterer Vorteil des hier vorgestellten Modells ist der zeitliche Verlauf der Lösung. Während der Algorithmus nach Kumar nur die Gleichgewichtslösung für das Flüssigkeits-Dampf-Gemisch zu jeden Zeitpunkt ermittelt, ist hier eine in der zeitlichen Entwicklung physikalisch plausible Lösung vorhanden, bei der die berechneten negativen Drücke im Kavitationsgebiet als Zugspannungen interpretiert werden. Durch das zeitliche Wachsen der Blasen steigt der Druck im Kavitationsgebiet an. In dem hier vorgestellten Modell hängt die Geschwindigkeit des Wachstums maßgeblich von der Druckdifferenz ab. Aber auch komplexere Modelle, bei denen das Wachstum mit von der Oberflächenspannung oder der Umgebungstemperatur abhängen, sind prinzipiell möglich. Die benötigte Rechenzeit für das zweidimensionale Problem lag bei 300 Sekunden und entspricht der CPU-Zeit der Vergleichsrechnung von FIRST. Allerdings kann der vorgestellte FEM-Algorithmus ohne großen Aufwand durch Einführung von Schrittweitensteuerung im Iterationsalgorithmus noch weiter beschleunigt werden. Des Weiteren ist die Berechnung des ersten Druckfelds, ausgehend von α 0 , noch sehr zeitintensiv im Vergleich zu den weiteren Iterationen. Dieses muss aber nicht, wie bisher geschehen, auskonvergiert werden, da es sich nur um eine Startlösung handelt. Zusammenfassung Im Rahmen dieser Arbeit wurde ein Verfahren zur Simulation von kavitierenden Strömungen vorgestellt und die dazugehörigen Gleichungen wurden mittels der Finite Elemente Methode diskretisiert. Das Verfahren wurde anhand eines eindimensionalen Modells von Geike und Popov [3] verifiziert. Für eine zweidimensionale Gleitschuhlagerung wurden Berechnungen durchgeführt und mit Ergebnissen des Algorithmus von Kumar [4] verglichen. Für die Lösung des Drucks konnte eine gute Übereinstimmung mit der Lösung von Kumar festgestellt werden, während für die Volumenfraktion die Ergebnisse nicht mit denen der Dichte aus Kumar vergleichbar waren. Mögliche Ursachen dafür ergeben sich aus der physikalischen Modellierung, insbesondere der nicht-konservativen Formulierung der Transportgleichung für den Dampfanteil sowie der Verwendung einer reduzierten Rayleigh-Plesset Gleichung für das Blasenwachstum. Der vorgestellte FEM-Algorithmus bietet jedoch die Möglichkeit entsprechende Modellerweiterungen einfach zu integrieren. Damit steht ein effizientes und flexibles Verfahren zur Verfügung, das für diverse Anwendungen eingesetzt werden kann. Im nächsten Schritt soll das Verfahren anhand von lagerspezifischen Belastungen, wie z. B. zeitlich veränderlichen Lasten und umlaufenden Bohrungen, bezüglich seiner physikalischen Modellierung und seiner rechenzeittechnischen Leistungsfähigkeit analysiert werden. Literaturverzeichnis [1] Bangerth, W., Heister, T., Heltai, L., Kanschat, G., Kronbichler, M., Maier, M., Turchsin, B., Young T.: „The deal.II library, Version 8.2“, Archive of Numerical Software, 3, 2015. [2] FIRST Simulationstool für elasto-hydrodynamische Mehrkörpersysteme, Benutzerhandbuch, IST Ingenieurgesellschaft für Strukuranalyse und Tribologie mbH (www.ist-aachen.com), 2015. [3] Geike, T., Popov, V.: „A Bubble Dynamics Based Approach to the Simulation of Cavitation in Lubricated Contacts“, Journal of Tribology, 131: 011704, 2009. Aus Wissenschaft und Forschung T+S_5_16 29.07.16 11: 27 Seite 30 Tribologie + Schmierungstechnik 63. Jahrgang 5/ 2016 [4] Kumar A.,: „Mass Conservation in Cavitating Bearings: A Finite Element Approach“, PhD Thesis, Cornell University, 1991. [5] Lang, R., Steinhilper, W.: „Gleitlager“, Springer-Verlag, 1. Edition, 1978. [6] Sauer, J.: „Instationär kavitierende Strömungen - Ein neues Modell, basierend auf front Capturing (VoF) und Blasendynamik“, Doktorarbeit, Universität Karlsruhe, 2000. [7] Van Husen, L.: „Finite Elemente Diskretisierung der Reynoldschen Differentialgleichung unter Berücksichtigung eines blasendynamischen Kavitationsmodells“, Masterarbeit, IGPM, RWTH Aachen, 2015. 31 Aus Wissenschaft und Forschung Impressum expert verlag GmbH: Wankelstr. 13, 71272 Renningen Postfach 20 20, 71268 Renningen Tel. (0 71 59) 92 65 - 0, Fax (0 71 59) 92 65 -20 E-Mail expert@expertverlag.de Vereinigte Volksbank AG, Sindelfingen BIC GENODES1 BBV, IBAN DE51 6039 0000 0032 9460 07 Postbank Stuttgart BIC PBNKDEFF, IBAN DE87 6001 0070 0022 5467 07 USt.-IdNr. DE 145162062 Anzeigen: Sigrid Hackenberg, expert verlag Tel. (0 71 59) 92 65 -13, Fax (0 71 59) 92 65-20 E-Mail anzeigen@expertverlag.de Informationen und Mediendaten senden wir Ihnen gerne zu. Vertrieb: Rainer Paulsen, expert verlag Tel. (0 71 59) 92 65 -16, Fax (0 71 59) 92 65-20 E-Mail paulsen@expertverlag.de Die zweimonatlich erscheinende Zeitschrift kostet bei Vorauszahlung im Jahresvorzugspreis für incl. Versand im Inland 189,- 7 (incl. 7 % MwSt.), im Ausland 198,- 7 * , Einzelheft 39,- 7; * (in der EU bei fehlender UID-Nr. zzgl. MwSt.); Studenten und persönliche Mitglieder der GfT erhalten gegen Vorlage eines entsprechenden Nachweises einen Nachlass von 20 % auf das Abo-Netto. Für Mitglieder der ÖTG ist der Abonnementspreis im Mitgliedschaftsbeitrag enthalten. Die Abonnementsgebühren sind jährlich im Voraus bei Rechnungsstellung durch den Verlag ohne Abzug zahlbar; kürzere Rechnungszeiträume bedingen einen Bearbeitungszuschlag von 3,- 7 pro Rechnungslegung. Abbestellungen müssen spätestens sechs Wochen vor Ende des Bezugsjahres schriftlich vorliegen. Der Bezug der Zeitschriften zum Jahresvorzugspreis verpflichtet den Besteller zur Abnahme eines vollen Jahrgangs. Bei vorzeitiger Beendigung eines Abonnementauftrages wird der Einzelpreis nachbelastet. Bei höherer Gewalt keine Lieferungspflicht. Erfüllungsort und Gerichtsstand: Leonberg expert verlag, 71272 Renningen ISSN 0724-3472 5/ 16 Tribologie und Schmierungstechnik Organ der Gesellschaft für Tribologie | Organ der Österreichischen Tribologischen Gesellschaft | Organ der Swiss Tribology Heft 5 September/ Oktober 2016 63. Jahrgang Herausgeber und Schriftleiter: Prof. Dr.-Ing. Dr. h.c. Wilfried J. Bartz Mühlhaldenstr. 91, 73770 Denkendorf Tel./ Fax (07 11) 3 46 48 35 E-Mail wilfried.bartz@tribo-lubri.de www.tribo-lubri.de Redaktion: Dr. rer. nat. Erich Santner, Bonn Tel. (02 28) 9 61 61 36 E-Mail esantner@arcor.de Redaktionssekretariat: expert verlag Tel. (0 71 59) 92 65 - 0, Fax (0 71 59) 92 65 -20 E-Mail: expert@expertverlag.de Beiträge, die mit vollem Namen oder auch mit Kurzzeichen des Autors gezeichnet sind, stellen die Meinung des Autors, nicht unbedingt auch die der Redaktion dar. Unverlangte Zusendungen redaktioneller Beiträge auf eigene Gefahr und ohne Gewähr für die Rücksendung. Die Einholung des Abdruckrechtes für dem Verlag eingesandte Fotos obliegt dem Einsender. Die Rechte an Abbildungen ohne Quellenhinweis liegen beim Autor oder der Redaktion. Ansprüche Dritter gegenüber dem Verlag sind, wenn keine besonderen Vereinbarungen getroffen sind, ausgeschlossen. Überarbeitungen und Kürzungen liegen im Ermessen der Redaktion. Die Wiedergabe von Gebrauchsnamen, Warenbezeichnungen und Handelsnamen in dieser Zeitschrift berechtigt nicht zu der Annahme, dass solche Namen ohne Weiteres von jedermann benutzt werden dürfen. Vielmehr handelt es sich häufig um geschützte, eingetragene Warenzeichen. Die Zeitschrift und alle in ihr enthaltenen Beiträge und Abbildungen sind urheberrechtlich geschützt. Mit Ausnahme der gesetzlich zugelassenen Fälle ist eine Verwertung ohne Einwilligung des Verlags strafbar. Dies gilt insbesondere für Vervielfältigungen, Übersetzungen, Mikroverfilmungen und die Einspeicherung und Verarbeitung in elektronischen Systemen. Entwurf und Layout: Ludwig-Kirn Layout, 71638 Ludwigsburg Impressum T+S_5_16 29.07.16 11: 27 Seite 31