CES Seminararbeit - RWTH Aachen University · 2021. 6. 8. · Fertigungsverfahren Prof. Dr.-Ing....
Transcript of CES Seminararbeit - RWTH Aachen University · 2021. 6. 8. · Fertigungsverfahren Prof. Dr.-Ing....
Lehrstuhl für Technologie der
Fertigungsverfahren Prof. Dr.-Ing. Dr.-Ing. E.h. Dr. h.c. Dr. h.c Fritz Klocke
CES Seminararbeit
Angefertigt von Joachim Stanke
Matr.-Nr.: 283273
Kurzthema: Entwicklung und Implementierung einer user
subroutine in Abaqus zur Visualisierung des
Reibkoeffizienten eines Blechumformprozesses
Betreuender Assistent: Dipl.-Ing. Dipl.-Wirt.Ing. Daniel Trauth
Aachen, den 22.12.2014
Inhalt und Ergebnis dieser Arbeit sind ausschließlich zum internen Gebrauch
bestimmt. Alle Urheberrechte liegen bei der RWTH Aachen. Ohne ausdrückliche
Genehmigung des betreuenden Lehrstuhls ist es nicht gestattet, diese Arbeit oder
Teile daraus an Dritte weiterzugeben.
Lehrstuhl für Technologie der
Fertigungsverfahren Prof. Dr.-Ing. Dr.-Ing. E.h. Dr. h.c. Dr. h.c Fritz Klocke
Aachen, 01.04.2014
D. Trauth - Tel. 0241-8027999
CES Seminararbeit
für Herrn.
Joachim Stanke
Matrikelnummer: 283273
Thema: Numerische Modellentwicklung des Tiefziehprozesses eines Türinnenblechs von
Haushaltsgeräten und Ableitung eines Analogiebauteils mit Hilfe der Finite-Elemente-
Methode.
Das Tiefziehen ist eines der bedeutendsten Fertigungsverfahren zur Herstellung von
Blechwerkstücken mit dreidimensionaler Geometrie. Die Formgebung der Blechwerkstücke wird dabei
maßgeblich von dem vorherrschenden tribologischen System, insbesondere den Reibeigenschaften,
beeinflusst. Die Reibung wiederum ist von der Kontaktnormalspannung, der Relativgeschwindigkeit
und der Temperatur abhängig.
Die Berücksichtigung der Reibung in Modellen auf Basis der Finite-Elemente-Methode (FEM) erfolgt in
der Regel anhand eines konstanten Reibkoeffizienten, weswegen es keine nativen Möglichkeiten und
Notwendigkeiten gibt, Reibkoeffizienten in einer FE-Software zu visualisieren. In vorherigen Arbeiten
am WZL wurde jedoch eine user subroutine (VFRICTION) für die FE-Software Abaqus entwickelt,
wodurch erstmalig für jeden Knotenpunkt der lokale Reibkoeffizient einer Kontaktpaarung als Funktion
der Kontaktnormalspannung, der Relativgeschwindigkeit und der Temperatur berechnet werden kann.
Hiermit wird die Visualisierung der lokalen Reibkoeffizienten in Konturplots und somit die Analyse von
kritischen Bereichen eines Blechumformbauteils erst ermöglicht.
Das Ziel dieser Arbeit ist daher die Entwicklung und Implementierung einer weiteren user subroutine in
Abaqus zur Visualisierung des Reibkoeffizienten eines Blechumformprozesses als Konturplot.
Die Vorgehensweise zur Erreichung der Zielsetzung basiert auf zwei Schritten. In einem ersten Schritt
erfolgt mithilfe eines vorhandenen FE-Modells zur numerischen Abbildung eines Streifenziehversuchs
die Entwicklung und Implementierung der user subroutine. Hierzu werden die vorhandene user
subroutine VFRICTION zur Berechnung der Reibkoeffizienten und die zu entwickelnde user
subroutine VUFIELD zur Visualisierung der berechneten Reibkoeffizienten benötigt. Die
Herausforderung besteht dabei in der Übergabe der Reibkoeffizienten von VFRICTION an VUFIELD
sowie die Verknüpfung der Reibkoeffizienten mit entsprechenden Knoteninformationen. Das FE-
Modell eines Streifenziehversuchs zur numerischen Analyse von Reibkoeffizienten erlaubt aufgrund
seiner geringen Komplexität die effiziente Entwicklung der user subroutine. In einem zweiten Schritt
erfolgt die Implementierung der entwickelten user subroutine VUFIELD in ein vorhandenes FE-Modell
eines Prozesses zur Herstellung einer Geschirrspülerinnentür. Für dieses FE-Modell liegen innerhalb
eines öffentlich geförderten Projekts umfangreiche numerische und experimentelle Vergleichsdaten
vor, welche eine Plausibilitätsprüfung der user subroutine ermöglichen. Im Fokus der
durchzuführendenen Simulation steht die Analyse des Reibkoeffizienten entlang eines definierten
Ziehweges. Hierdurch kann überprüft werden, ob der Reibkoeffizient lokalen Extremwerten ausgesetzt
ist und als Ursache für mögliche Tiefziehfehler ausfindig gemacht werden kann.
Zur Erfüllung der Aufgabenstellung erfolgt eine Einarbeitung in die Simulationssoftware Abaqus, die
vorhandenen FE-Modelle und die Programmiersprache Fortran. Die Ergebnisse werden in einer
schriftlichen Arbeit wissenschaftlich aufbereitet und zusammengefasst.
Im Einzelnen sind die folgenden Teilaufgaben zu lösen:
Einarbeitung in die numerische Simulationssoftware Abaqus, die vorhandenen FE-Modelle und die Programmiersprache Fortran
Programmierung einer user subroutine durch Implementierung der user subroutines VFRICTION und VUFIELD in einem FE-Modell eines Streifenziehversuchs für Entwicklungszwecke
Programmierung geeigneter Funktionsblöcke zur Übergabe der Reibkoeffizienten und der Knoteninformationen zwischen VFRICTION und VUFIELD
Durchführung einer Plausibilitätsprüfung der user subroutine mithilfe eines vorhandenen FE-Modells eines Streifenziehversuchs durch Abgleichen der Simulationsergebnisse mit vorhandenen Versuchsdaten
Übertragung der user subroutine auf ein FE-Modell zur Simulation einer Geschirrspülerinnentürfertigung und Analyse des Reibkoeffizienten entlang eines definierten Ziehweges
Schriftliche Dokumentation der Arbeit
Inhaltsverzeichnis i
Inhaltsverzeichnis
Inhaltsverzeichnis .................................................................................................................................. i
Formelzeichen und Abkürzungen ........................................................................................................ ii
Abbildungsverzeichnis ........................................................................................................................ iv
1 Einleitung ....................................................................................................................................... 1
2 Stand der Technik ......................................................................................................................... 2
2.1 Grundlagen des Tiefziehens .................................................................................................. 2
2.2 Grundlagen des Streifenziehversuchs ................................................................................... 4
2.3 Reibmodell nach Filzek .......................................................................................................... 5
2.4 Abaqus und User Subroutinen ............................................................................................... 5
3 Entwicklung der User Subroutinen am Streifenziehversuch ................................................... 7
3.1 User Subroutinen zum Plotten von statischen Reibkoeffizienten .......................................... 7
3.2 User Subroutinen zum Plotten von Reibkoeffizienten nach Filzeks Reibmodell ................. 14
3.3 Erweiterte Konturplot-Optionen ............................................................................................ 16
4 Übertragung der User Subroutinen auf einen industriellen Tiefziehprozess ....................... 21
4.1 Ergebnisübersicht der Reibkoeffizienten ............................................................................. 21
4.2 Aktuelle Reibkoeffizienten .................................................................................................... 21
4.3 Gemittelte Reibkoeffizienten ................................................................................................ 23
5 Zusammenfassung und Fazit .................................................................................................... 25
Literaturverzeichnis ............................................................................................................................ 26
Anhang ................................................................................................................................................. 27
Formelzeichen und Abkürzungen ii
Formelzeichen und Abkürzungen
Formelzeichen Einheit Beschreibung
𝜇, 𝜇0 Reibkoeffizient
𝑣𝑠𝑡 mm/s Stempelgeschwindigkeit
𝑣𝑍 mm/s Ziehringgeschwindigkeit
𝜑1 Hauptumformung
𝜑2 Nebenumformung
𝐾𝑓 MPa Fließspannung
𝐹 kN Kraft
𝑟, r-Werte nach Lankford
𝑡 s Zeit
𝜗, 𝜗0 °C Temperatur
𝜆 W/m∙K Wärmeleitfähigkeit
𝑐𝑝𝑚 kJ/kg∙K Mittlere spezifische Wärmekapazität
𝛼𝑚 K-1
Mittlerer linearer Ausdehnungskoeffizient
𝜌 kg/m3
Dichte
𝐸 MPa Elastizitätsmodul
𝑣𝑧 Mm/s Ziehringgeschwindigkeit
Δ𝑡 s Schrittdauer
Δ𝑠 mm Verschiebung
𝑝, 𝑝0 MPa Druck
𝑇, 𝑇0 °C Temperatur
𝑣, 𝑣0 mm/s Relativgeschwindigkeit
Formelzeichen und Abkürzungen iii
𝑛,𝑚, 𝑘 Koeffizienten im Reibmodell nach Filzek
Abkürzung Beschreibung
DIN Deutsches Institut für Normung
FE Finite-Elemente
FEM Finite-Elemente-Methode
USR User Subroutine
Abbildungsverzeichnis iv
Abbildungsverzeichnis
Abbildung 2.1 Prinzipskizze des Tiefziehprozesses [KLOC06]............................................................... 2
Abbildung 2.2 Aufbau und Prozessbeschreibung des Tiefziehmodells .................................................. 3
Abbildung 2.3 FE-Modell des Streifenziehversuchs ................................................................................ 4
Abbildung 3.1 Kontaktdefinition in Abaqus\CAE ................................................................................... 12
Abbildung 3.2 Konturplot am SZV mit statischen Reibkoeffizienten ..................................................... 13
Abbildung 3.3 Kontaktdefinition in der Abaqus\CAE für das Reibmodell nach Filzek .......................... 14
Abbildung 3.4 Konturplot am SZV mit dem Reibmodell nach Filzek ..................................................... 16
Abbildung 3.5 Ablaufplan der Subroutinen VFRICTION und VUFIELD ................................................ 19
Abbildung 3.6 Reibwertergebnisse am SZV mit 4 Visualisierungsvarianten ........................................ 20
Abbildung 4.1 Reibwertergebnisse am Tiefziehmodell mit 4 Visualisierungsvarianten ........................ 21
Abbildung 4.2 Verlauf der Reibkoeffizienten am Tiefziehmodell über die Prozessschritte ................... 22
Abbildung 4.3 Detailbetrachtung der Reibkoeffizienten zum Prozessende .......................................... 22
Abbildung 4.4 Verlauf der gemittelten Reibkoeffizienten am Tiefziehmodell über die Prozessschritte 23
Abbildung 4.5 Detailbetrachtung der gemittelten Reibkoeffizienten zum Prozessende ....................... 24
1 Einleitung 1
1 Einleitung
Auf Grund seiner weiten Verbreitung in der Groß- und Kleinserienproduktion, z.B. in der Maschinen-,
Luftfahrt-, Fahrzeug- und Haushaltsindustrie, ist der Tiefziehprozess einer der wichtigsten
Blechumformprozesse. Die gefertigten Blechteile können eine Größe von nur wenigen Millimetern, bis
hin zu mehreren Metern aufweisen. Außerdem können sie verschiedensten Anforderungen genügen,
da nahezu alle metallischen Werkstoffe als Verarbeitungsmaterial für das Blechwerkstück gewählt
werden können [DOEG10].
Die Simulation mittels Finite-Elemente-Methode (FEM) hat sich als wertvolles Werkzeug etabliert, um
zusätzliche Erkenntnisse über den Tiefziehprozess zu erhalten. Die FEM kann in der Simulation von
Blechumformprozessen unter anderem Einblicke in den Materialzustand geben, welche auf
experimentellem Wege nicht oder nur unter sehr großem Aufwand gewonnen werden können. Um
aussagekräftige Ergebnisse zu erhalten, ist aber eine genaue Abbildung der Blechumformprozesse
notwendig. Dazu gehört z.B. die Abbildung der tribologischen Bedingungen. Dies geschieht in der
Regel über die Wahl eines geeigneten Reibgesetzes. Stand der Dinge in der Simulation von
Blechumformprozessen ist die Verwendung des Coulombschen Reibgesetzes. Bei diesem ist die
Reibkraft nur Abhängig von der Normalkraft. Es wurde aber bereits in mehreren Arbeiten gezeigt,
dass die Reibkräfte von einer Vielzahl von Parametern abhängen [FILZ13].
Diese Arbeit befasst sich mit einem Finite-Elemente (FE)-Modell eines Tiefziehprozesses, welches ein
von Filzek vorgestelltes Reibgesetz verwendet. Bei diesem wird ein Reibkoeffizient in Abhängigkeit
von der Relativgeschwindigkeit, der Normalspannung und der Temperatur berechnet. Als Simula-
tionssoftware dient die FE-Software Abaqus 6.14. Ziel ist es, die nach Filzeks Reibgesetz berechneten
Reibkoeffizienten auf dem Blechwerkstück in einem Konturplot abzubilden. Die Umsetzung erfolgt
mittels Programmierung von Abaqus User Subroutinen (USR).
Im ersten Schritt wird eine USR programmiert, um einen statischen Reibkoeffizienten auf einem
Blechstreifen eines Streifenziehversuchs zu plotten. Der Streifenziehversuch stellt einen
Analogieprozess zum Tiefziehen dar, ist aber stark vereinfacht. Dadurch eignet er sich sehr gut, um
die Implementierung einer USR zu testen, bevor sie auf das industrielle Tiefziehmodell übertragen
wird. Anschließend wird die USR so modifiziert, dass die nach Filzeks Reibgesetz berechneten
Reibkoeffizienten auf den Blechstreifen geplottet werden. Diese USR wird dann auf das
Tiefziehmodell übertragen und die Resultate ausgewertet.
2 Stand der Technik 2
2 Stand der Technik
In diesem Abschnitt werden die technischen Grundlagen dieser Arbeit vermittelt. Die FE-Modelle des
Streifenziehversuchs sowie des Tiefziehversuches werden vorgestellt. Ebenso erfolgt eine kurze
Vorstellung der FE-Software Abaqus sowie den Abaqus User Subroutinen (USRs). Am Ende des
Abschnittes wird das Reibmodell nach Filzek präsentiert, welches in den verwendeten FE-Modellen
Anwendung fand.
2.1 Grundlagen des Tiefziehens
Der Tiefziehprozess besitzt eine bedeutende Position unter den Herstellungsverfahren zur Produktion
von Blechwerkstücken mit dreidimensionaler Geometrie. Dieses Herstellungsverfahren wird Industrie-
zweigen wie der Automobil-, der Flugzeug- oder der Verpackungsindustrie eingesetzt [KLOC06]. Nach
DIN 8584-3 ist das Tiefziehen folgendermaßen definiert [DIN03]: „Tiefziehen ist das
Zugdruckumformen eines Blechzuschnitts (je nach Werkstoff auch eine Folie oder Platte, eines
Ausschnittes oder Abschnitts) zu einem Hohlkörper oder eines Hohlkörpers zu einem Hohlkörper mit
kleinerem Umfang ohne eine beabsichtigte Änderung der Blechdicke.“ Der prinzipielle
Werkzeugaufbau für das Tiefziehen eines runden Napfes besteht aus einem Niederhalter, einem
Ziehring und einem Stempel. Ein Schnitt durch ein schematisches Modell dieses Prozesses ist in
Abbildung 2.1 zu sehen.
Abbildung 2.1 Prinzipskizze des Tiefziehprozesses [KLOC06]
Das Werkstück wird vor Beginn des Prozesses zwischen Niederhalter und Ziehring mit einer
definierten Kraft eingespannt, wodurch die Faltenbildung in diesem Bereich verhindert wird. Der
Stempel befindet sich im oberen Totpunkt. Skizziert ist dies in Abbildung 2.1 links. Während des
Arbeitsschrittes fährt der Stempel herunter, wobei er das Blech durch die Öffnung im Ziehring zieht.
Mit Erreichen des unteren Totpunkts ist das Tiefziehen abgeschlossen. Je nach Anforderung kann das
Werkstück einen Flansch behalten, oder es wird vom Stempel vollständig eingezogen. In Abbildung
2.1 ist ein Ziehteil / Blechwerkstück mit Flansch skizziert [KLOC06].
Die Einleitung der Umformung erfolgt durch ein Streckziehen an der Stirnfläche des Stempels, sobald
dieser das Blech durch die Öffnung im Ziehring zieht. Es findet eine Superpositionierung von Tiefzieh-
2 Stand der Technik 3
und Streckziehprozess beim Übergang vom Boden in die Zarge des Napfes statt. Der eigentliche
Tiefziehvorgang geschieht im Flansch sowie im Bereich der Ziehkante, also im Kontaktbereich von
Blech und Ziehring [KUWE07].
Das Tiefzieh-/ FE-Modell
Das Tiefziehmodell dieser Arbeit ist ein FE-Modell eines dreistufigen FE-Prozesses zur Herstellung
einer Geschirrspülerinnentür. In Abbildung 2.2 ist eine Abbildung des Modells mit einer Übersicht über
die drei grundlegenden Prozessschritte zu sehen. Das Modell besteht aus einem Stempel, einem
Ziehring, einem Niederhalter und zusätzlich einem Gegenhalter. Der Gegenhalter gibt zusammen mit
dem Stempel dem Blech seine Endkontur. Im ersten Schritt werden der Ziehring und der Niederhalter
um 13 mm gesenkt. Darauf folgt das Senken des Stempels, bis das Blech den Boden des
Gegenhalters erreicht. Im letzten Schritt werden noch einmal der Ziehring und der Niederhalter bis
zum Erreichen der Endkontur herabgesenkt. Anfangs- und Randbedingungen des FE-Modells
stammen aus Daten über den realen Prozess, der mit diesem Simulationsmodell abgebildet wird.
Abbildung 2.2 Aufbau und Prozessbeschreibung des Tiefziehmodells
Tabelle 2.1 soll einen Überblick über das verwendete Materialmodell des Bleches geben. Mit einer
Fließkurve wurden die plastischen Verformungseigenschaften modelliert. Das Grenzformänderungs-
2 Stand der Technik 4
schaubild dient als Schadensmodell. Das anisotrope Dehnungsverhalten wurde im Materialmodell
ebenfalls mittels der r-Werte berücksichtigt.
Tabelle 2.1 Materialmodell des Blechwerkstoffes, siehe BA Joachim Stanke
Weitere physikalische Eigenschaften, wie die Wärmeleitfähigkeit, die Dichte und das Elastizitätsmodul
wurden temperaturabhängig modelliert, da die Simulationen temperaturgekoppelt durchgeführt
werden. Weiterführende Informationen können der Bachelorarbeit .von Joachim Stanke [STAN13]
entnommen werden.
2.2 Grundlagen des Streifenziehversuchs
Der Streifenziehversuch dient als Analogiemodell zum Tiefziehversuch. In Experimenten wird er
beispielsweise benutzt, um Reibkoeffizienten oder die Koeffizienten eines Reibgesetzes zu
bestimmen. In dieser Arbeit wird ein FE-Modell des Streifenziehversuches genutzt, um die
programmierten USRs zu testen, bevor sie auf das Tiefziehmodell übertragen werden.
Abbildung 2.3 FE-Modell des Streifenziehversuchs
2 Stand der Technik 5
Das Modell besteht aus einem Blechstreifen, der zwischen zwei Niederhalterplatten gehalten wird,
siehe Abbildung 2.3 Dabei wird eine konstante Haltekraft über die obere Niederhalterplatte in
negative z-Richtung aufgebracht, während die untere Niederhalterplatte in alle Koordinatenrichtungen
fest ist. Der Prozess besteht darin, den Blechstreifen in negative y-Richtung zwischen den beiden
Niederhalterplatten hindurch zu ziehen. Der Blechstreifen besitzt die gleiche Dicke und das gleiche
Materialmodell wie das Blechwerkstück im Tiefziehmodell.
2.3 Reibmodell nach Filzek
Die Reibungsverhältnisse haben im Tiefziehprozess einen entscheidenden Einfluss auf die
mechanischen und geometrischen Eigenschaften, sowie die Oberflächenqualität der Bauteile. Das
macht eine realitätsnahe Abbildung in der FEM Simulation unumgänglich, wenn es darum geht, mit
numerischen Methoden aussagekräftige Ergebnisse zu erhalten. Filzek stellt deswegen ein
empirisches Reibgesetz speziell für den Tiefziehprozess vor. Ausgehend von einem Potenzansatz
sieht sein Reibgesetz bei einer reinen Druckabhängigkeit wie folgt aus:
𝜇 = 𝜇0 ∙ (𝑝
𝑝0)𝑛−1
(1)
Dabei bezeichnet 𝜇 den Reibkoeffizient und 𝑝 den Druck bzw. die Flächenpressung. 𝑝0 und 𝜇0 sind
konstante Bezugswerte. Der Exponent 𝑛 ist ebenfalls konstant und im Bereich 0 < 𝑛 ≤ 1. Um die
Relativgeschwindigkeit 𝑣 und die Temperatur 𝑇 in das Gesetz miteinzubeziehen, kann man den
Ansatz auf folgende Form erweitern:
𝜇 = 𝜇0 ∙ (𝑝
𝑝0)𝑛−1
∙ (𝑣
𝑣0)𝑚−1
∙ (𝑇
𝑇0)𝑘−1
(2)
Um dieses Gesetz praktisch nutzen zu können, müssen die konstanten Bezugswerte für die
Temperatur, die Relativgeschwindigkeit, die Flächenpressung und den Reibkoeffizienten festgelegt
werden sowie die Koeffizienten 𝑛, 𝑚, und 𝑘 aus Experimenten gewonnen werden [FILZ13].
2.4 Abaqus und User Subroutinen
Abaqus ist ein vor allem in der Industrie weit verbreitetes, kommerzielles Software-Paket zur
Modellierung und Analyse verschiedener Probleme, wie:
Festkörper-
Fluiddynamik-
Wärmeleitungsprobleme
Die Modellierung geschieht über eine grafische Benutzeroberfläche, der Abaqus\CAE. In dieser
Benutzeroberfläche stehen aber nicht alle Optionen und Möglichkeiten zur Verfügung, welche
Abaqus dem Anwender bietet. Durch die Arbeit mit der Abaqus Input File und den User Subroutinen
kann der Benutzer auf ein noch größeres Spektrum an Optionen und Funktionen zugreifen. In den
Input Files werden bestimmte Optionen mit Hilfe von Schlüsselworten aktiviert. Abaqus User
Subroutinen sind Programme, die in Fortran 77 geschrieben werden und in Abaqus miteingebunden
werden. Abaqus bietet eine Reihe von USRs für verschiedene Anwendungsmöglichkeiten. Ein paar
Beispiele für solche Subroutinen sind:
VUMAT: Erstellen von Materialmodellen
2 Stand der Technik 6
VFRICTION: Erstellen von Reibungsmodellen
VUFIELD: Modifizieren von Feldvariablen
Die USRs VFRICTION und VUFIELD werden im folgenden Kapitel genauer erläutert, da diese Teil
dieser Arbeit sind.
3 Entwicklung der User Subroutinen am Streifenziehversuch 7
3 Entwicklung der User Subroutinen am
Streifenziehversuch
In diesem Kapitel werden die Subroutinen zum Plotten der Reibkoeffizienten auf dem Blechstreifen
des Streifenziehversuches vorgestellt. Hierbei werden zwei Abaqus USRs eine Rolle spielen, die
VFRICTION zur Berechnung der Reibkoeffizienten und die VUFIELD zum Erstellen einer
benutzerdefinierten Field-Variable. Zum Einstieg wird im ersten Schritt eine Subroutine erstellt, die
lediglich einen statischen Reibkoeffizient nutzt, welcher dann auf dem Blechstreifen geplottet wird.
Anschließend wird das gleiche Prinzip angewendet, um Reibkoeffizienten zu plotten, welche nach
Filzeks Reibmodell berechnet wurden.
3.1 User Subroutinen zum Plotten von statischen Reibkoeffizienten
Die Berechnung der Reibkoeffizienten geschieht in der Abaqus USR VFRICTION. Es wird in diesem
Abschnitt nur ein konstanter Reibkoeffizient für die Berechnung der Reibung genutzt und dieser
anschließend in die VUFIELD USR übergeben, damit sich die korrekte Funktionsweise der
Subroutinen leicht überprüfen lässt. Aus demselben Grund wird als Basis für die VFRICTION
Subroutine eine von Abaqus verifizierte Subroutine aus dem Handbuch [DASS14] benutzt. Wie in
Abschnitt 2.4 erläutert, werden die Subroutinen in der Programmiersprache Fortran 77 geschrieben.
Im Folgenden wird zuerst die USR VFRICTION beschrieben, im Anschluss daran die VUFIELD USR.
Dann folgen die Beschreibungen zur Kontaktdefinition in der Abaqus\CAE und die notwendigen
Änderungen an der Input File.
c 1 c User subroutine VFRICTION to define friction forces 2 c 3 subroutine VFRICTION ( 4 c Write only - 5 * fTangential, 6 c Read/Write - 7 * state, 8 c Read only - 9 * nBlock, nBlockAnal, nBlockEdge, 10 * nNodState, nNodSlv, nNodMst, 11 * nFricDir, nDir, 12 * nStates, nProps, nTemp, nFields, 13 * jFlags, rData, 14 * surfInt, surfSlv, surfMst, 15 * jConSlvUid, jConMstUid, props, 16 * dSlipFric, fStickForce, fTangPrev, fNormal, 17 * areaCont, dircosN, dircosSl, 18 * shapeSlv, shapeMst, 19 * coordSlv, coordMst, 20 * velSlv, velMst, 21 * tempSlv, tempMst, 22 * fieldSlv, fieldMst ) 23
In den Zeilen 1 bis 23 werden der Name der Subroutine und die Aufrufparameter definiert. Es wird im
Folgenden ausschließlich eine Erklärung für die wichtigsten Variablen gegeben, welche auch in der
Programmierung im späteren Verlauf genutzt werden.
fTangential(nFricDir,nBlock): Ein Array, in welches die berechneten Werte für die Reibkraft
geschrieben werden
3 Entwicklung der User Subroutinen am Streifenziehversuch 8
nBlock: Entspricht der Anzahl der Knoten, die in diesem Funktionsaufruf bearbeitet werden
nBlockAnal: Ist in dieser Anwendung gleich nBlock
nFricDir: Anzahl der Richtungskomponenten der Reibkraft. In diesem Fall zwei
nProps: Anzahl der vom Benutzer über die Input File/CAE übergebenen Variablen
jFlags(1): Schrittnummer
jFlags(2): Inkrementnummer
rData(3): Größe des aktuellen Zeitinkrements
jConSlvUid(1,nBlock): Knotennummer des Knotens der Slave-Oberfläche
props(nprops): Vom Benutzer mittels Input File/CAE übergebene Variablen
dSlipFric(nDir,nBlock): Die Verschiebung, welche die Knoten innerhalb des aktuellen
Zeitinkrements erfahren haben
fNormal(nBlock): Array mit den Kräften in Normalenrichtung der Knoten
areaCont(nBlock): Array mit der Größe der Kontaktflächen, mit denen die Knoten in Verbindung
stehen
tempSlv(nBlock): Array mit den Temperaturen der Slave-Knoten
tempMst(nBlockAnal): Array mit den Temperaturen der Master-Knoten
c 24 c Array dimensioning variables: 25 c 26 c nBlockAnal = nBlock (non-analytical-rigid master surface) 27 c nBlockAnal = 1 (analytical rigid master surface) 28 c nBlockEdge = 1 (non-edge-type slave surface) 29 c nBlockEdge = nBlock (edge-type slave surface) 30 c nNodState = 1 (node-type slave surface) 31 c nNodState = 4 (edge-type slave surface) 32 c nNodSlv = 1 (node-type slave surface) 33 c nNodSlv = 2 (edge-type slave surface) 34 c nNodMst = 4 (facet-type master surface) 35 c nNodMst = 2 (edge-type master surface) 36 c nNodMst = 1 (analytical rigid master surface) 37 c 38 c Surface names surfSlv and surfMst are not available for 39 c general contact (set to blank). 40 c 41 include 'vaba_param.inc' 42 dimension fTangential(nFricDir,nBlock), 43 * state(nStates,nNodState,nBlock), 44 * jConSlvUid(nNodSlv,nBlock), 45 * jConMstUid(nNodMst,nBlockAnal), 46 * props(nProps), 47 * dSlipFric(nDir,nBlock), 48 * fStickForce(nBlock), 49 * fTangPrev(nDir,nBlock), 50 * fNormal(nBlock), 51
3 Entwicklung der User Subroutinen am Streifenziehversuch 9
* areaCont(nBlock), 52 * dircosN(nDir,nBlock), 53 * dircosSl(nDir,nBlock), 54 * shapeSlv(nNodSlv,nBlockEdge), 55 * shapeMst(nNodMst,nBlockAnal), 56 * coordSlv(nDir,nNodSlv,nBlock), 57 * coordMst(nDir,nNodMst,nBlockAnal), 58 * velSlv(nDir,nNodSlv,nBlock), 59 * velMst(nDir,nNodMst,nBlockAnal), 60 * tempSlv(nBlock), 61 * tempMst(nBlockAnal), 62 * fieldSlv(nFields,nBlock), 63 * fieldMst(nFields,nBlockAnal) 64 c 65 parameter( iKStep = 1, 66 * iKInc = 2, 67 * iLConType = 3, 68 * nFlags = 3 ) 69 parameter( iTimStep = 1, 70 * iTimGlb = 2, 71 * iDTimCur = 3, 72 * iFrictionWork = 4, 73 * nData = 4 ) 74 c 75 dimension jFlags(nFlags), rData(nData) 76 character*80 surfInt, surfSlv, surfMst 77 parameter( zero=0.d0 ) 78 c 79 integer ProbSize 80 parameter (ProbSize=100000) 81 c 82
In Zeile 80 beginnt der Benutzer-Code. Die Variable ProbSize dient als Variable zur
Dimensionierung einiger Arrays, die im anschließenden Teil definiert werden. Sie entspricht der
Knotenanzahl im Blechstück des Modells.
real ContMu(ProbSize) 83 common /kContact/ ContMu 84 c 85
Um die Reibkoeffizienten aus der USR VFRICTION in die USR VUFIELD zu transferieren, werden
common Blocks genutzt. Common Blocks sind in Fortran 77 geteilte Speicherbereiche, auf die jedes
Fortran 77 Unterprogramm zugreifen kann, wenn der betreffende common Block mit dem „common“
Schlüsselbegriff verfügbar gemacht wurde. In Zeile 84 enden die Variablen Deklarationen.
u = props(1) 86 c 87
Wie bereits erwähnt, können vom Benutzer über die Input File/CAE Eingangsvariablen bestimmt
werden. In Zeile 86 wird die erste dieser Variablen an u übergeben. u entspricht dem statischen
Reibkoeffizienten.
do k = 1, nBlock 88 fn = fNormal(k) 89 fs = fStickForce(k) 90 ContMu(jConSlvUid(1,k))=u 91 ft = min ( ContMu(jConSlvUid(1,k))*fn, fs ) 92 fTangential(1,k) = -ft 93
3 Entwicklung der User Subroutinen am Streifenziehversuch 10
fTangential(2,k) = zero 94 end do 95 c 96
In Zeile 88 beginnt eine Schleife über alle Knoten, die in diesem Aufruf der Subroutine bearbeitet
werden. In der Zeile 89 wird der Wert der Normalkraft an fn übergeben. fStickForce(k) in Zeile
90 enthält die Kräfte, die notwendig sind, um die Haftbedingungen zu erhalten. In der nächsten Zeile
wird der statische Reibkoeffizient an die common Block Variable ContMu übergeben. Der Wert des
Reibkoeffizienten wird im Array ContMu genau an die Stelle geschrieben, die der Knoten-Nummer
entspricht. Dafür sorgt die Verwendung von jConSlvUid(1,k). Dies macht in der VUFIELD USR
ein Suchen nach dem Reibkoeffizienten zu einem bestimmten Knoten überflüssig. Die Bestimmung
des Minimums in Zeile 92 dient der Erhaltung der Haftbedingung.
return 97 end 98 c 99
In Zeile 99 endet die VFRICTION Subroutine. Als nächstes folgt die Beschreibung der VUFIELD
Subroutine, mit der die benutzerdefinierte Field-Variable mit den Reibkoeffizienten aus der
VFRICTION aktualisiert wird, um den angestrebten Konturplot zu erhalten.
SUBROUTINE VUFIELD(FIELD, NBLOCK, NFIELD, KFIELD, NCOMP, 100 1 KSTEP, JFLAGS, JNODEID, TIME, 101 2 COORDS, U, V, A) 102 c 103
In Zeile 100 bis 102 befindet sich die Definition der Aufrufparameter der Subroutine. In dieser
Subroutine sind folgende Variablen für unsere Anwendung wichtig:
FIELD(NBLOCK, NCOMP, NFIELD): Field ist das Array, welches in dieser Subroutine upgedatet
wird
NBLOCK: Analog zur vorherigen Subroutine entspricht NBLOCK der Anzahl der in diesem Aufruf der
Subroutine zu bearbeitenden Knoten
NCOMP: Ist bei der Verwendung von Kontinuums-Elementen gleich 1
NFIELD: Entspricht der Anzahl der definierten Field-Variablen und ist in dieser Version der Subroutine
gleich 1
TIME (i_ufld_Total=4): Gesamtzeit der Simulation
JNODEID(NBLOCK): Array mit den Knotennummern die in diesem Aufruf der USR bearbeitet werden
INCLUDE 'VABA_PARAM.INC' 104 105 C indices for the time array TIME 106 PARAMETER( i_ufld_Current = 1, 107 * i_ufld_Increment = 2, 108 * i_ufld_Period = 3, 109 * i_ufld_Total = 4 ) 110 111 C indices for the coordinate array COORDS 112 PARAMETER( i_ufld_CoordX = 1, 113 * i_ufld_CoordY = 2, 114 * i_ufld_CoordZ = 3 ) 115
3 Entwicklung der User Subroutinen am Streifenziehversuch 11
116 C indices for the displacement array U 117 PARAMETER( i_ufld_SpaDisplX = 1, 118 * i_ufld_SpaDisplY = 2, 119 * i_ufld_SpaDisplZ = 3, 120 * i_ufld_RotDisplX = 4, 121 * i_ufld_RotDisplY = 5, 122 * i_ufld_RotDisplZ = 6, 123 * i_ufld_AcoPress = 7, 124 * i_ufld_Temp = 8 ) 125 126 C indices for the velocity array V 127 PARAMETER( i_ufld_SpaVelX = 1, 128 * i_ufld_SpaVelY = 2, 129 * i_ufld_SpaVelZ = 3, 130 * i_ufld_RotVelX = 4, 131 * i_ufld_RotVelY = 5, 132 * i_ufld_RotVelZ = 6, 133 * i_ufld_DAcoPress = 7, 134 * i_ufld_DTemp = 8 ) 135 136 C indices for the acceleration array A 137 PARAMETER( i_ufld_SpaAccelX = 1, 138 * i_ufld_SpaAccelY = 2, 139 * i_ufld_SpaAccelZ = 3, 140 * i_ufld_RotAccelX = 4, 141 * i_ufld_RotAccelY = 5, 142 * i_ufld_RotAccelZ = 6, 143 * i_ufld_DDAcoPress = 7, 144 * i_ufld_DDTemp = 8 ) 145 146 C indices for JFLAGS 147 PARAMETER( i_ufld_kInc = 1, 148 * i_ufld_kPass = 2 ) 149 C 150 DIMENSION FIELD(NBLOCK,NCOMP,NFIELD) 151 DIMENSION JFLAGS(2), JNODEID(NBLOCK), TIME(4), 152 * COORDS(3,NBLOCK) 153 DIMENSION U(8,NBLOCK), V(8,NBLOCK), A(8,NBLOCK) 154 C 155
In den oben stehenden Zeilen stehen für die Subroutine notwendige Dimensionierungen und
Parameterdefinitionen.
integer ProbSize 156 Parameter (ProbSize=100000) 157 c 158
In Zeile 156 beginnt der benutzerdefinierte Code der VUFIELD Subroutine. Als erstes wird wieder die
Variable „ProbSize“ deklariert und als Parameter festgelegt.
real ContMu(ProbSize) 159 common /kContact/ ContMu 160 c 161
Die common Block Variable ContMu muss auch in dieser Subroutine deklariert werden. Aus ihr erhält
die VUFIELD USR die Reibkoeffizienten aus der VFRICTION Subroutine.
if (TIME(i_ufld_Total) .eq. 0.0) then 162 do k = 1, ProbSize 163
3 Entwicklung der User Subroutinen am Streifenziehversuch 12
ContMu(k) = 0.0 164 end do 165 c 166
Zeile 162 bis 166 dienen der Initialisierung von ContMu und werden nur ausgeführt, wenn die
Simulationszeit gleich Null ist, also zu Beginn der Simulation.
else 167 do k = 1,NBLOCK 168 FIELD(k, 1, 1) = KontMu(JNODEID(k)) 169 end do 170 end if 171 c 172
In der Schleife von Zeile 168 bis 170 werden die Reibkoeffizienten von der KontMu Variable in die
Field Variable geschrieben.
RETURN 173 END 174
In Zeile 174 endet die Subroutine.
Um die Abaqus USRs nutzen zu können, sind noch Änderungen bzw. Ergänzungen sowohl in der
Abaqus Input File als auch über die Abaqus\CAE notwendig. Eine Vorgehensweise rein über die
Abaqus Input File wäre ebenso möglich, soll hier aber nicht erläutert werden. Nachdem das Abaqus
Modell in der Abaqus\CAE geöffnet ist, wird, wie in Abbildung 3.1 beschrieben, die Kontaktdefinition
aktualisiert. Zuerst öffnet man das „Edit Contact Property“ Menü durch Doppelklick mit der rechten
Maustaste auf die betreffende „Contact Property“. In zweiten Schritt wählt man mittels Drop Down
Menü „User-defined“. Damit definiert man, dass man eine benutzerdefinierte Subroutine zur
Berechnung der Reibungskräfte benutzen möchte. Im dritten und letzten Schritt geben wir einen Wert
für „Friction Properties“ an. In diesem Fall „0.1“. Die hier definierten Werte werden der Subroutine in
der zuvor erwähnten „props(n)“ Variable zur Verfügung gestellt.
Abbildung 3.1 Kontaktdefinition in Abaqus\CAE
Nachdem die Kontaktdefinition abgeschlossen ist, schreibt man die Abaqus Input File. An dieser Datei
werden jetzt noch zwei wichtige Änderungen durchgeführt. Im Folgenden ist der erste von zwei
Ausschnitten der Input File zu sehen.
3 Entwicklung der User Subroutinen am Streifenziehversuch 13
** INTERACTION PROPERTIES 1 ** 2 *Surface Interaction, name=Coulomb_0-137 3 *Friction, user=friction, depvar=0, properties=1 4 0.1 5
Zeile 1 bis 5 zeigt die Kontaktdefinition in der Input File. Wichtig ist das Hinzufügen von „=friction“ in
Zeile 4 hinter „user“. Damit wird die VFRICTION Subroutine als USR zur Berechnung der Reibkräfte
festgelegt. Als Nächstes ist die zweite Ergänzung zu sehen.
*FIELD, USER 1 BlechSet-C3D8 2
Die beiden Zeilen müssen in der Input File zu der Definition jeden Schrittes hinzugefügt werden. Mit
diesen beiden Zeilen wird das benutzerdefinierte Feld erzeugt, welches in der VUFIELD Subroutine
mit den Reibkoeffizienten upgedatet wird. Dabei bezeichnet „BlechSet-C3D8“ den Namen eines
Knoten-Sets, welches mit den Knoten des Blechstreifens in der Abaqus\CAE manuell generiert wurde.
In Abbildung 3.2 ist ein Konturplot der Field Variable zu sehen, an die die Reibkoeffizienten übergeben
wurden. Die Simulation befindet sich zu diesem Zeitpunkt in etwa der halben Simulationszeit. Die
obere Niederhalterplatte wurde ausgeblendet, um das Ergebnis besser betrachten zu können. Die
untere Niederhalterplatte wurde in der Ansicht belassen, um zu beurteilen, wo sich beide Platten
aktuell befinden. Man sieht, dass die Übergabe des statischen Reibkoeffizienten funktioniert.
Zwischen den beiden Platten besitzt die Field Variable den Wert 0,1, was dem Reibkoeffizienten
entspricht. Es wird ebenfalls festgestellt, dass der Bereich des Blechstreifens, der noch keinen
Kontakt mit den Niederhalteplatten hatte, den Wert 0 besitzt. Aber der Bereich, der keinen Kontakt
mehr mit den Niederhalterplatten hat, wird nicht wieder auf 0 zurückgesetzt, was ein unerwünschtes
Verhalten darstellt. Grund hierfür ist, dass die VFRICTION USR nur bei Kontakt aufgerufen wird und
es keinen Mechanismus zum Zurücksetzen der Werte für kontaktlose Knoten gibt. Dies wird später
noch korrigiert werden. Die aktuelle Variante ist dennoch als Analysewerkzeug verwendbar und wird
später beibehalten. Dazu in Abschnitt 3.3 mehr.
Abbildung 3.2 Konturplot am SZV mit statischen Reibkoeffizienten
Zuvor wird im nächsten Abschnitt das Reibmodell nach Filzek aus Abschnitt 2.3 in die VFRICTION
USR implementiert.
3 Entwicklung der User Subroutinen am Streifenziehversuch 14
3.2 User Subroutinen zum Plotten von Reibkoeffizienten nach
Filzeks Reibmodell
Um das Reibmodell nach Filzek zu implementieren, werden lediglich Änderungen an der VFRICTION
Subroutine vorgenommen. Die VUFIELD USR bleibt in diesem Abschnitt unverändert. Es werden nur
die Teile der Subroutine beschrieben, die bearbeitet wurden. Da nun das Reibmodell nach Filzek
genutzt wird, ändern sich aber auch die Werte, die in der Abaqus\CAE in die Kontaktdefinition
eingetragen werden. Die geänderten Koeffizienten sind aus diesem Grunde in Abbildung 3.1
dargestellt.
Abbildung 3.3 Kontaktdefinition in der Abaqus\CAE für das Reibmodell nach Filzek
Die Input File ändert sich im Vergleich zum letzten Abschnitt nicht. In der VFRICTION Subroutine
muss die Schleife über die einzelnen Knoten neu geschrieben werden.
mu0=props(1) 1 a=props(2) 2 press0=props(3) 3 b=props(4) 4 ve0=props(5) 5 c=props(6) 6 temp0=props(7) 7
Bevor die Schleife beginnen kann, werden die Eingabeparameter aus der Kontaktdefinition an lokale
Variablen übergeben. Dies findet in den Zeilen 1 bis 7 statt.
do k = 1, nBlock 8 disp=sqrt(dSlipFric(1,k)**2.0 + 9 * dSlipFric(2,k)**2.0 + 10 * dSlipFric(3,k)**2.0) 11
3 Entwicklung der User Subroutinen am Streifenziehversuch 15
In Zeile 8 beginnt die Schleife über die Knoten dieses Subroutine-Aufrufes. Als erstes wird mit disp in
der Schleife die Verschiebung des Knotens berechnet.
if (fNormal(k) .ne. 0.0 .and. 12 * disp .ne. 0.0) then 13
In Zeile 12 folgt eine „if“ Abfrage für den Wert der Normalkraft und den der Verschiebung, um
Divisionen durch Null bei der Berechnung des Reibkoeffizienten zu vermeiden.
press=fNormal(k)/areaCont(k) 14 ve=disp/rData(iDTimCur) 15 temp=(tempMst(k)+tempSlv(k))/2.0 16
In den Zeilen 14 bis 16 werden drei Berechnungen durchgeführt. Zuerst in Zeile 14 die Berechnung
des Drucks (press). Diese ergibt sich aus der Normalkraft (fNormal(k)) dividiert durch die
Kontaktfläche (areaCont(k)). Die Verschiebung (disp) geteilt durch die Größe des Zeitinkrements
(rData(iDTimCur)), ergibt in Zeile 15 die Relativgeschwindigkeit und in Zeile 16 wird die
Temperatur berechnet. Diese wird als gemittelter Wert aus der Temperatur am Master-Knoten
(tempMst(k)) und der Temperatur am Slave Knoten (tempSlv(k)) berechnet.
ContMu(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0) 17 * *(ve/ve0)**(b-1.0) 18 * *(temp/temp0)**(c-1.0) 19
In den Zeilen 17 bis 19 erfolgt die Berechnung des Reibkoeffizienten nach Filzeks Reibgesetz. Dieser
wird an das common Block Array ContMu übergeben.
ft = min ( ContMu(jConSlvUid(1,k))*fNormal(k), 20 * fStickForce(k) ) 21
In den Zeilen 20 bis 21 wird die Reibkraft analog zur Subroutine im letzten Abschnitt bestimmt. Nur
dieses Mal mit dem vorher berechneten Reibkoeffizienten, an Stelle eines statischen Koeffizienten.
fTangential(1,k) = -ft 22 fTangential(2,k) = 0.0 23
In Zeile 22 wird die berechnete Reibkraft an die erste Tangentialkomponente übergeben, die zweite
wird wieder auf 0 gesetzt, da sie orthogonal zu Bewegungsrichtung liegt.
else 24 ContMu(jConSlvUid(1,k)) = 0.0 25
Waren in der if Abfrage aus Zeile 12 die Normalkraft oder die Verschiebung gleich 0 wird der Wert 0 in
Zeile 25 an das common Block Array übergeben.
end if 26 end do 27
In Zeile 27 endet die Schleife zur Berechnung der Reibkoeffizienten nach Filzeks Reibgesetz.
Es folgt nun eine Betrachtung der Ergebnisse dieses Reibmodells am SZV. In Abbildung 3.4 ist das
Ergebnis der Berechnung, wie im vorherigen Abschnitt auch, in der Mitte der Simulation zu sehen. Es
ist erkennbar, dass einzelne Knoten mit sehr niedrigen oder sehr hohen Reibwerten in der
Berechnung entstehen. Dies liegt an den Randbedingungen des Modells und numerischen
Ungenauigkeiten bei der Berechnung. Im Laufe der Simulation treten Schwingungen an der oberen
Niederhalterplatte und am Blechstreifen auf. Weil aber das SZV Modell nur genutzt wird, um die
Subroutinen zum Erstellen eines Konturplots zu entwickeln, stellt dieses kein Problem dar. Abgesehen
von den Maximal- und Minimalwerten kann in der Simulation ein durchschnittlicher Reibwert von ca.
0,15 festgestellt werden.
3 Entwicklung der User Subroutinen am Streifenziehversuch 16
Abbildung 3.4 Konturplot am SZV mit dem Reibmodell nach Filzek
Wie bereits im ersten Abschnitt des Kapitels erwähnt, soll der Konturplot noch dahingehend modifiziert
werden, dass die Stellen, bei dem der Kontakt unterbrochen wurde, wieder auf 0 zurückgesetzt
werden. Diese und zwei weitere Optionen für den Konturplot werden im nächsten Abschnitt entwickelt.
3.3 Erweiterte Konturplot-Optionen
In diesem Abschnitt werden drei neue Konturplot-Optionen vorgestellt: Zum Einem die Möglichkeit,
nur die Reibkoeffizienten aus dem aktuellem Zeitinkrement zu plotten, zum Nächsten noch ein
Konturplot, der zeitlich gemittelten Reibkoeffizienten und schließlich ein Konturplot der maximalen
Reibzahlen in der Simulation. Zu Beginn eine notwendige Änderung an der Input File:
*FIELD, USER, NUMBER=4 1 BlechSet-C3D8 2
In Zeile 1 des Ausschnittes der Input File wurde hinter den „Field“ Schlüsselbegriff noch der Zusatz
„NUMBER=4“ hinzugefügt. Dieses ermöglicht, vier verschiedene Fieldvariablen zu nutzen. Alles
weitere bleibt an der Input File unverändert. Als nächstes müssen Modifikationen an beiden
Subroutinen vorgenommen werden. Als erstes wird die erweiterte Version der VFRICTION
beschrieben. Wie im vorherigen Abschnitt auch, ist die Schleife zur Berechnung der Reibkoeffizienten
zu sehen. Dieses Mal - wie folgt beschrieben - mit einigen Erweiterungen.
mu0=props(2) 1 a=props(3) 2 press0=props(4) 3 b=props(5) 4 ve0=props(6) 5 c=props(7) 6 temp0=props(8) 7 do k = 1, nBlock 8 disp=sqrt(dSlipFric(1,k)**2.0 + 9 * dSlipFric(2,k)**2.0 + 10 * dSlipFric(3,k)**2.0) 11 if (fNormal(k) .ne. 0.0 .and. 12 * disp .ne. 0.0) then 13 press=fNormal(k)/areaCont(k) 14
3 Entwicklung der User Subroutinen am Streifenziehversuch 17
ve=disp/rData(iDTimCur) 15 temp=(tempMst(k)+tempSlv(k))/2.0 16 ContMu(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0) 17 * *(ve/ve0)**(b-1.0) 18 * *(temp/temp0)**(c-1.0) 19 WeightedMu(jConSlvUid(1,k))=WeightedMu(jConSlvUid(1,k))+ 20 * ContMu(jConSlvUid(1,k))* 21 * rData(3) 22 TimeSum(jConSlvUid(1,k))=TimeSum(jConSlvUid(1,k))+ 23 * rData(3) 24
Bis Zeile 20 ist die Subroutine identisch zur vorherigen Version. In Zeile 20 und 23 werden zwei neue
common Block Arrays verwendet. Beide dienen als Mittel zur Berechnung eines gemittelten
Reibkoeffizienten. „WeightedMu“ ist die Summe aus den Produkten von Reibkoeffizient und
Zeitinkrement. Anders dargestellt:
𝜇𝑠𝑢𝑚,𝑘 = ∑ 𝜇𝑘,𝑖𝑖=𝑖𝑛𝑘# ∙ ∆𝑡𝑖 (3)
Wobei hier 𝜇𝑠𝑢𝑚,𝑘 die Werte des Arrays WeightedMu bezeichnen. 𝜇𝑘,𝑖 entspricht dem Reibkoeffizient
des i-ten Zeitinkrements und ∆𝑡𝑖 dem dazugehörige Zeitinkrementen. Das zweite Array, TimeSum
speichert die aufsummierten Werte der Zeitinkremente.
ft = min ( ContMu(jConSlvUid(1,k))*fNormal(k), 25 * fStickForce(k) ) 26 fTangential(1,k) = -ft 27 fTangential(2,k) = 0.0 28 else 29 ContMu(jConSlvUid(1,k)) = 0.0 30 end if 31 CurrentNodeInc(jConSlvUid(1,k))=jFlags(2) 32 CurrentNodeStep(jConSlvUid(1,k))=jFlags(1) 33 end do 34
Die zweite Änderung ist in Zeile 32 und 33 zu sehen. CurrentNodeInc und CurrentNodeStep
sind ebenfalls common Block Arrays. Sie speichern die Nummer des aktuellen Zeitinkrements und des
aktuellen Schrittes, der Simulation für die Knoten, mit denen die Subroutine in diesem Durchlauf
aufgerufen wurde. Diese werden in der VUFIELD Subroutine genutzt, um einen Konturplot nur mit
aktuellen Reibkoeffizienten zu erhalten.
In der VUFIELD Subroutine werden vier Schleifen definiert, um die vier verschiedenen Field Variablen
zu aktualisieren. Diese vier Schleifen werden nun beschrieben.
do k = 1,NBLOCK 1 FIELD(k, 1, 1) = ContMu(JNODEID(k)) 2 end do 3
Die erste Schleife der Subroutine erzeugt den schon aus dem letzten Abschnitt bekannten Konturplot.
3 Entwicklung der User Subroutinen am Streifenziehversuch 18
do k = 1,NBLOCK 4 if (CurrentNodeInc(JNODEID(k)) .eq. CurrentInc .and. 5 * CurrentNodeStep(JNODEID(k)) .eq. CurrentStep ) then 6 FIELD(k, 1, 2) = ContMu(JNODEID(k)) 7 else 8 FIELD(k, 1, 2) = 0.0 9 end if 10 end do 11
Die zweite Schleife erzeugt einen Konturplot nur mit den aktuellen Reibkoeffizienten. Knoten an denen
aktuell kein Kontakt herrscht, werden auf 0 gesetzt. Um zu gewährleisten, dass nur die aktuellen
Werte in die Field Variable übernommen werden, wird in den Zeilen 5 und 6 eine „if“ Abfrage
durchgeführt. In dieser werden die Werte aus den common Block Array CurrentNodeInc und
CurrentNodeStep mit dem aktuellen Werten des Inkrementes und Schrittes verglichen. So wird
überprüft, ob der Wert aus dem ContMu Array ein aktuell berechneter oder ein Wert aus einem
Vorherigem ist, welcher nicht im Konturplot erscheinen soll.
do k= 1,NBLOCK 12 if (FIELD(k, 1, 3) .lt. ContMu(JNODEID(k))) then 13 FIELD(k, 1, 3)=ContMu(JNODEID(k)) 14 end if 15 end do 16
In der dritten Schleife wird ein Konturplot der maximalen Reibkoeffizienten erzeugt. Hierfür wird eine
„if“ Abfrage in der 13. Zeile durchgeführt. Wenn der zuletzt berechnete Wert des Reibkoeffizienten
größer ist als der aktuell in der Field-Variable gespeicherte Wert, wird dieser auf den neuen Wert in
Zeile 14 gesetzt.
do k= 1,NBLOCK 17 if (TimeSum(JNODEID(k)) .ne. 0.0) then 18 FIELD(k, 1, 4)=WeightedMu(JNODEID(k))/TimeSum(JNODEID(k)) 19 else 20 FIELD(k, 1, 4)=0.0 21 end if 22 end do 23
In der vierten Schleife werden Mittelwerte der Reibkoeffizienten berechnet. Dies geschieht mit den
beiden, in der VFRICTION erzeugten common Block Arrays WeightedMu und TimeSum. Die
Mittelwerte ergeben sich aus der Division von WeightedMu durch TimeSum in Zeile 19. In Zeile 18
findet zuvor noch eine „if“-Abfrage statt, um sicher zu stellen, dass keine Divisionen durch Null
entstehen. Der Schematische Ablauf der beiden Subroutinen kann in der folgenden Abbildung
betrachtet werden.
3 Entwicklung der User Subroutinen am Streifenziehversuch 19
Abbildung 3.5 Ablaufplan der Subroutinen VFRICTION und VUFIELD
In Abbildung 3.6 kann das Ergebnis dieser Subroutinen betrachtet werden.
3 Entwicklung der User Subroutinen am Streifenziehversuch 20
Abbildung 3.6 Reibwertergebnisse am SZV mit 4 Visualisierungsvarianten
Oben links ist das Ergebnis zu sehen, dass schon aus dem vorherigen Abschnitt bekannt ist. Bei der
Betrachtung des Ergebnisses rechts daneben ist zu sehen, dass nun keine von 0 verschiedenen
Reibkoeffizienten in den Bereichen ohne Kontakt zu sehen sind. Das Bild links unten zeigt die
maximalen Reibkoeffizienten, die im Simulationsverlauf an den einzelnen Knoten berechnet wurden.
Im letzten Bild können die zeitlich gemittelten Reibkoeffizienten betrachtet werden. Genauere
Betrachtungen werden im nächsten Kapitel stattfinden, da der Streifenziehversuch nur zur
Entwicklung der Subroutinen diente. Im nächsten Kapitel werden die Ergebnisse untersucht, die sich
mit diesen Subroutinen am industriellen Tiefziehmodell ergeben.
4 Übertragung der User Subroutinen auf einen industriellen Tiefziehprozess 21
4 Übertragung der User Subroutinen auf einen
industriellen Tiefziehprozess
In diesem Kapitel wird die im vorherigen Kapitel entwickelte Abaqus User Subroutine auf das
Tiefziehmodell aus Abschnitt 2.1 angewendet. Zuerst wird eine Übersicht über die Ergebnisse aller
vier Konturplot-Varianten betrachtet. Darauf folgt eine genauere Betrachtung zweier ausgewählter
Varianten.
4.1 Ergebnisübersicht der Reibkoeffizienten
Werden die Ergebnisse aus den vier Field Variablen in Abbildung 4.1 verglichen, werden deutliche
Unterschiede festgestellt. In der Abbildung oben links sind die letzten Reibwerte zu sehen, d. h. die
aktuellen und die zuletzt berechneten Werte. Hier sind die Maximalwerte in der Innenfläche zu
erkennen. Dies ist der Bereich, der von Stempel und Gegenhalter im zweiten Schritt des
Tiefziehprozesses gebildet wird. Im zweiten Bild oben rechts sind nur die aktuellen Reibkoeffizienten
zu sehen. Hier sind die von Null verschiedenen Reibkoeffizienten im Innenbereich deutlich weniger.
Unten links werden die Maximalwerte aus den Reibkoeffizienten gezeigt. Im gesamten Einzugsbereich
treten die Maximalwerte 𝜇 = 0,147 auf. Dies liegt an der zunächst niedrigen Kontaktnormalspannung
im Schließprozess der Werkzeuge. Bei geringem Druck treten die höchsten Reibwerte auf. Unten
rechts werden die gemittelten Reibwerte gezeigt.
Abbildung 4.1 Reibwertergebnisse am Tiefziehmodell mit 4 Visualisierungsvarianten
Die gemittelten und die aktuellen Reibkoeffizienten werden in den nächsten beiden Abschnitten
genauer betrachtet.
4.2 Aktuelle Reibkoeffizienten
Zuerst wird der Verlauf der Reibkoeffizienten über die Enden der drei Tiefziehschritte betrachtet. In
Abbildung 4.2 ist eine Übersicht zu sehen. Es ist deutlich zu erkennen, dass die höchsten
4 Übertragung der User Subroutinen auf einen industriellen Tiefziehprozess 22
Reibkoeffizienten im Einzugsbereich im zweiten Schritt auftreten. Am Ende des ersten Schrittes treten
ebenfalls lokal hohe Reibkoeffizienten im Außenbereich des Bleches auf. Zum Ende des dritten
Schrittes, also zum Ende des Prozesses, sind die Reibwerte abgeklungen.
Abbildung 4.2 Verlauf der Reibkoeffizienten am Tiefziehmodell über die Prozessschritte
Das Ergebnis zum Prozessende soll nun noch genauer in Abbildung 4.3 betrachtet werden.
Abbildung 4.3 Detailbetrachtung der Reibkoeffizienten zum Prozessende
Im Innenbereich sind zu Prozessende nur schmale Kontaktbereiche mit Reibkoeffizienten ungleich
Null zu beobachten. Im oberen Bereich können bei genauerer Betrachtung maximale Reibwerte von
ca. 𝜇 = 0,12 beobachtet werden. Im seitlichen Bereich sind es sogar Werte von bis zu 𝜇 = 0,13. Im
Außenbereich sind breite Bereiche mit Reibwerten ungleich Null mit lokal schwankenden Werten zu
4 Übertragung der User Subroutinen auf einen industriellen Tiefziehprozess 23
sehen. Die Maximalwerte sind hier ca. 𝜇 = 0,11, wobei die größten Bereiche mit solchen Reibwerten
hier an den Ecken zu erkennen sind.
4.3 Gemittelte Reibkoeffizienten
Die Betrachtung der gemittelten Reibkoeffizienten ist besonders aufschlussreich, da sie Informationen
über die Reibwertverteilung des gesamten Prozessverlaufs liefern. Zuerst erfolgt wieder ein Überblick
über die Ergebnisse in Abbildung 4.4.
Abbildung 4.4 Verlauf der gemittelten Reibkoeffizienten am Tiefziehmodell über die
Prozessschritte
Im Außenbereich treten die höchsten Reibwerte hier am Ende des ersten Schrittes auf. Danach
sinken sie ab und besitzen am Ende des zweiten und dritten Schrittes ungefähr gleiche Werte. Im
Innenbereich bilden sich vom ersten bis zum letzten Schritt des Prozesses zunehmend größere
Kontaktbereiche mit Reibkoeffizienten ungleich Null.
4 Übertragung der User Subroutinen auf einen industriellen Tiefziehprozess 24
Abbildung 4.5 Detailbetrachtung der gemittelten Reibkoeffizienten zum Prozessende
Die Ergebnisse des letzten Schritts werden in Abbildung 4.5 detaillierter gezeigt. Im Innenbereich
zeigen sich die höchsten Reibwerte. In Ecken und Rändern werden Reibwerte von bis zu 𝜇 = 0.14 mit
lokal teils stark schwankenden Werten beobachtet. Im Einzugsbereich ist die Verteilung sehr
homogen. Im unteren und oberen Bereich sowie den seitlichen Bereichen liegen die Reibkoeffizienten
bei einem Wert von ca. 𝜇 = 0,11. In den Ecken ist der Wert mit ca. 𝜇 = 0,12.leicht höher.
5 Zusammenfassung und Fazit 25
5 Zusammenfassung und Fazit
Mit dieser Arbeit konnte gezeigt werden, dass das Plotten von benutzerdefinierten Reibkoeffizienten in
Abaqus unter Verwendung von Abaqus Subroutinen in einem Blechumformprozess möglich ist. Dabei
können nicht nur statische Reibkoeffizienten auf einem Blechstreifen visualisiert werden, wie als
Einführung im ersten Abschnitt des Kapitels 3 gezeigt wurde, auch die Visualisierung eines nach
einem dynamischen Reibgesetz berechneten Reibkoeffizienten ist möglich und wurde erfolgreich im
zweiten Abschnitt des 3. Kapitels durchgeführt. Wenn es darum geht, die Visualisierungsmöglich-
keiten für eine Analyse zu optimieren, kommen verschiedene Konturplot-Varianten in Frage. Im letzten
Abschnitt des 3. Kapitels wurden vier solcher Varianten vorgestellt. Als besonders nützliche Variante
ist hier der Konturplot der gemittelten Reibkoeffizienten zu nennen, da er unabhängig von zeitlichen
Schwankungen ist, welche auf Grund von numerischen Ungenauigkeiten lokal stark ausfallen können.
In Kapitel Fehler! Verweisquelle konnte nicht gefunden werden. wurde gezeigt, dass sich die an
einem Streifenziehversuch entwickelten Subroutinen ebenso auf ein komplexes Tiefziehmodell
anwenden lassen. Zur Analyse wurden aus den vier Konturplot-Varianten die der aktuellen
Reibkoeffizienten und die der gemittelten Reibkoeffizienten gewählt. Mit deren Hilfe konnten Bereiche
mit besonders hohen Reibwerten ausgemacht werden.
Diese Ergebnisse zeigen, dass Konturplots von Reibkoeffizienten bei der Simulation von
Blechumformprozessen nicht nur möglich sind, sondern auch wertvolle Informationen für
Prozessoptimierungen liefern können. Vorausgesetzt, das Simulationsmodell und natürlich besonders
das Reibungsmodell sind validiert und verhalten sich realitätsnah.
In der Zukunft werden neue Subroutinen entwickelt werden, die das Plotten der Reibkoeffizienten
direkt am Werkzeug des Simulationsmodells ermöglichen werden. Diese Methode wird die Prozedur,
Bereiche auszumachen, in denen Optimierungen am Werkzeug den Prozess positiv beeinflussen
werden, weiter erleichtern.
Literaturverzeichnis 26
Literaturverzeichnis
[DASS14] Dassault Systèmes: Abaqus 6.14 Documentation, 2014.
[DIN03] DIN 8584-3: Fertigungsverfahren Zugdruckumformen, Teil 3 Tiefziehen.
Normenausschuss Technische Grundlagen, September 2003.
[DOEG10] Doege, E.; Behrens, B.-A.: Handbuch Umformtechnik. 2. Aufl., Berlin, Heidelberg:
Springer-Verlag Berlin Heidelberg, 2010.
[FILZ13] Filzek, J.: Berücksichtigung der Reibung in der FEM-Simulation. In (Europäische
Forschungsgesellschaft für Blechverarbeitung e.V. Hrsg.): Umformen, Schneiden,
Verbinden im Leichtbau. Machbarkeit-Produktivität-Leichtbau, 2013; S. 359–370.
[KLOC06] Klocke, F.: Fertigungsverfahren 4. 5. Aufl., Berlin, Heidelberg: Springer-Verlag Berlin
Heidelberg, 2006.
[KUWE07] Kuwer, C.-J.: Verschleißreduktion beim Tiefziehen von X 5 CrNi 18-10. Diss. RWTH
Aachen, 2007.
[STAN13] Stanke, J.: Numerische Modellentwicklung des Tiefziehprozesses eines Türinnen-
blechs von Haushaltsgeräten und Ableitung eines Analogiebauteils mit Hilfe der Finite-
Elemente-Methode. BA. RWTH-Aachen, 2013.
Anhang 27
Anhang
A.1 Subroutinen für die Verwendung mit Continuumselementen
c
c User subroutine VFRICTION to define friction forces
c
subroutine vfriction (
c Write only -
* fTangential,
c Read/Write -
* state,
c Read only -
* nBlock, nBlockAnal, nBlockEdge,
* nNodState, nNodSlv, nNodMst,
* nFricDir, nDir,
* nStates, nProps, nTemp, nFields,
* jFlags, rData,
* surfInt, surfSlv, surfMst,
* jConSlvUid, jConMstUid, props,
* dSlipFric, fStickForce, fTangPrev, fNormal,
* areaCont, dircosN, dircosSl,
* shapeSlv, shapeMst,
* coordSlv, coordMst,
* velSlv, velMst,
* tempSlv, tempMst,
* fieldSlv, fieldMst )
c
c Array dimensioning variables:
c
c nBlockAnal = nBlock (non-analytical-rigid master surface)
c nBlockAnal = 1 (analytical rigid master surface)
c nBlockEdge = 1 (non-edge-type slave surface)
c nBlockEdge = nBlock (edge-type slave surface)
c nNodState = 1 (node-type slave surface)
c nNodState = 4 (edge-type slave surface)
c nNodSlv = 1 (node-type slave surface)
c nNodSlv = 2 (edge-type slave surface)
c nNodMst = 4 (facet-type master surface)
c nNodMst = 2 (edge-type master surface)
c nNodMst = 1 (analytical rigid master surface)
c
c Surface names surfSlv and surfMst are not available for
c general contact (set to blank).
c
include 'vaba_param.inc'
dimension fTangential(nFricDir,nBlock),
* state(nStates,nNodState,nBlock),
* jConSlvUid(nNodSlv,nBlock),
* jConMstUid(nNodMst,nBlockAnal),
* props(nProps),
* dSlipFric(nDir,nBlock),
* fStickForce(nBlock),
* fTangPrev(nDir,nBlock),
* fNormal(nBlock),
* areaCont(nBlock),
* dircosN(nDir,nBlock),
* dircosSl(nDir,nBlock),
* shapeSlv(nNodSlv,nBlockEdge),
* shapeMst(nNodMst,nBlockAnal),
* coordSlv(nDir,nNodSlv,nBlock),
Anhang 28
* coordMst(nDir,nNodMst,nBlockAnal),
* velSlv(nDir,nNodSlv,nBlock),
* velMst(nDir,nNodMst,nBlockAnal),
* tempSlv(nBlock),
* tempMst(nBlockAnal),
* fieldSlv(nFields,nBlock),
* fieldMst(nFields,nBlockAnal)
c
parameter( iKStep = 1,
* iKInc = 2,
* iLConType = 3,
* nFlags = 3 )
parameter( iTimStep = 1,
* iTimGlb = 2,
* iDTimCur = 3,
* iFrictionWork = 4,
* nData = 4 )
c
dimension jFlags(nFlags), rData(nData)
character*80 surfInt, surfSlv, surfMst
c
ccc-----------------------------------------------------------------ccc
ccc Benutzer Code zur Definition der VFRICTION Subroutine ccc
ccc-----------------------------------------------------------------ccc
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc \/ Variablen Deklarationen \/ ccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
ccc ProbSize bestimmt die Groeße der common Blocks.
ccc ProbSize muss mindestens 2* Anz. der Knoten im zu untersuchenden
ccc Node Set sein (Bspw. 2* Knoten Anz. des Blechstücks).
c
integer ProbSize
parameter (ProbSize=1100000)
c
ccc KontMu ist das Array zur Uebergabe der Reibungskoeffizienten aus
ccc der VFRICTION in die VUFIELD
c
real KontMu(ProbSize)
c
ccc WeightedMu dient als Array zur Berechnung eines Durchschnitts-
ccc wertes der Reibungskoeffizienten
c
double precision WeightedMu(ProbSize)
c
ccc Zur Berechnung des Mittelwertes wird die Summe ueber die Zeit-
ccc inkremente genutzt, bei der ein Reibungskoeffizient ungleich
ccc 0 berechnet wurde
c
real TimeSum(ProbSize)
c
ccc Einige lokale Variablen zur Berechnung und Definition der
ccc physikalischen Parameter: Druck, Temperatur,
ccc Exponenten: a, b, c, etc.
c
real temp0, ve0, press0, mu0, a, b, c, press, ve, temp, ft, disp
c
ccc KurentNodeInc und KurentNodeStep zeigen von welchem Inkrement und
ccc Schritt die Reibungskoeffizienten in KuntMu stammen.
ccc KurentInc und KurentStep beinhalten die aktuelle Inkrement Nr.
ccc und Step Nr..
c
Anhang 29
integer KurentNodeInc(ProbSize), KurentNodeStep(ProbSize),
* KurentInc, KurentStep
c
ccc Definition der lokalen Variable Abh, dient der Bestimmung welche
ccc Berechnungsvorschrift für den Reibungskoeffizient genutzt werden
ccc soll.
integer Abh
c
ccc KontMu, KurentNodeInc, KurentNodeStep, KurentInc, KurentStep und
ccc HistoryOpt werden als common Blocks definiert.
c
common /KontaktReal/ KontMu, WeightedMu, TimeSum
common /KTime/ KurentNodeInc, KurentInc, KurentNodeStep,
* KurentStep
c
ccc save befehl stellt sicher, dass die common Blocks nach verlassen
ccc der Subroutine erhalten bleibt.
c
save /KontaktReal/
save /KTime/
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc /\ Ende der Variablen Deklarationen /\ ccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
ccc Die aktuelle Inkrement Nr. und Schritt Nr. wird Upgedatet
c
KurentInc=jFlags(2)
KurentStep=jFlags(1)
c
ccc Die gewaehlte Berechnungsvorschrift für die
ccc Reibungskoeffizientenwird werden aus der Abaqus Inputfile an
ccc Abh in Form eines numerischen Wertes uebergeben.
c
Abh=props(1)
c
if ( Abh .eq. 1) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird druckabheangig berechnet ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
a=props(3)
press0=props(4)
do k = 1, nBlock
if (fNormal(k) .ne. 0.0) then
press=fNormal(k)/areaCont(k)
KontMu(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0)
WeightedMu(jConSlvUid(1,k))=WeightedMu(jConSlvUid(1,k))+
* KontMu(jConSlvUid(1,k))*
* rData(3)
TimeSum(jConSlvUid(1,k))=TimeSum(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMu(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
KontMu(jConSlvUid(1,k)) = 0.0
end if
KurentNodeInc(jConSlvUid(1,k))=jFlags(2)
KurentNodeStep(jConSlvUid(1,k))=jFlags(1)
end do
else if( Abh .eq. 2) then
Anhang 30
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird rel. geschwindigkeitsabhaengig ccc
ccc berechnet ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
b=props(3)
ve0=props(4)
do k = 1, nBlock
disp=sqrt(dSlipFric(1,k)**2.0 +
* dSlipFric(2,k)**2.0 +
* dSlipFric(3,k)**2.0)
if (fNormal(k) .ne. 0.0 .and.
* rData(iDTimCur) .ne. 0.0 .and.
* disp .ne. 0.0) then
ve=disp/rData(iDTimCur)
KontMu(jConSlvUid(1,k))=mu0*(ve/ve0)**(b-1.0)
WeightedMu(jConSlvUid(1,k))=WeightedMu(jConSlvUid(1,k))+
* KontMu(jConSlvUid(1,k))*
* rData(3)
TimeSum(jConSlvUid(1,k))=TimeSum(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMu(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
KontMu(jConSlvUid(1,k)) = 0.0
end if
KurentNodeInc(jConSlvUid(1,k))=jFlags(2)
KurentNodeStep(jConSlvUid(1,k))=jFlags(1)
end do
else if( Abh .eq. 3) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird temperaturabhaengig berechnet. ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
c=props(3)
temp0=props(4)
do k = 1, nBlock
if (fNormal(k) .ne. 0.0) then
temp=(tempMst(k)+tempSlv(k))/2.0
KontMu(jConSlvUid(1,k))=mu0*(temp/temp0)**(c-1.0)
WeightedMu(jConSlvUid(1,k))=WeightedMu(jConSlvUid(1,k))+
* KontMu(jConSlvUid(1,k))*
* rData(3)
TimeSum(jConSlvUid(1,k))=TimeSum(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMu(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
KontMu(jConSlvUid(1,k)) = 0.0
end if
KurentNodeInc(jConSlvUid(1,k))=jFlags(2)
KurentNodeStep(jConSlvUid(1,k))=jFlags(1)
end do
else if( Abh .eq. 4) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird druck- und rel ccc
ccc geschwindigkeitsabhaengig berechnet. ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
Anhang 31
a=props(3)
press0=props(4)
b=props(5)
ve0=props(6)
do k = 1, nBlock
disp=sqrt(dSlipFric(1,k)**2.0 +
* dSlipFric(2,k)**2.0 +
* dSlipFric(3,k)**2.0)
if (fNormal(k) .ne. 0.0 .and.
* areaCont(k) .ne. 0.0 .and.
* rData(iDTimCur) .ne. 0.0 .and.
* disp .ne. 0.0) then
press=fNormal(k)/areaCont(k)
ve=disp/rData(iDTimCur)
KontMu(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0)
* *(ve/ve0)**(b-1.0)
WeightedMu(jConSlvUid(1,k))=WeightedMu(jConSlvUid(1,k))+
* KontMu(jConSlvUid(1,k))*
* rData(3)
TimeSum(jConSlvUid(1,k))=TimeSum(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMu(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
KontMu(jConSlvUid(1,k)) = 0.0
end if
KurentNodeInc(jConSlvUid(1,k))=jFlags(2)
KurentNodeStep(jConSlvUid(1,k))=jFlags(1)
end do
else if(Abh .eq. 5) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird geschwindigkeits- und ccc
ccc temperaturabhaengig berechnet. ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
a=props(3)
press0=props(4)
c=props(5)
temp0=props(6)
do k = 1, nBlock
if (fNormal(k) .ne. 0.0) then
press=fNormal(k)/areaCont(k)
temp=(tempMst(k)+tempSlv(k))/2.0
KontMu(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0)
* *(temp/temp0)**(c-1.0)
WeightedMu(jConSlvUid(1,k))=WeightedMu(jConSlvUid(1,k))+
* KontMu(jConSlvUid(1,k))*
* rData(3)
TimeSum(jConSlvUid(1,k))=TimeSum(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMu(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
KontMu(jConSlvUid(1,k)) = 0.0
end if
KurentNodeInc(jConSlvUid(1,k))=jFlags(2)
KurentNodeStep(jConSlvUid(1,k))=jFlags(1)
end do
else if(Abh .eq. 6) then
Anhang 32
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird geschwindigkeits- und ccc
ccc temperaturabhaengig berechnet. ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
b=props(3)
ve0=props(4)
c=props(5)
temp0=props(6)
do k = 1, nBlock
disp=sqrt(dSlipFric(1,k)**2.0 +
* dSlipFric(2,k)**2.0 +
* dSlipFric(3,k)**2.0)
if (fNormal(k) .ne. 0.0 .and.
* rData(iDTimCur) .ne. 0.0 .and.
* disp .ne. 0.0) then
ve=disp/rData(iDTimCur)
temp=(tempMst(k)+tempSlv(k))/2.0
KontMu(jConSlvUid(1,k))=mu0*(ve/ve0)**(b-1.0)
* *(temp/temp0)**(c-1.0)
WeightedMu(jConSlvUid(1,k))=WeightedMu(jConSlvUid(1,k))+
* KontMu(jConSlvUid(1,k))*
* rData(3)
TimeSum(jConSlvUid(1,k))=TimeSum(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMu(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
KontMu(jConSlvUid(1,k)) = 0.0
end if
KurentNodeInc(jConSlvUid(1,k))=jFlags(2)
KurentNodeStep(jConSlvUid(1,k))=jFlags(1)
end do
else if(Abh .eq. 7) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird druck-, rel. geschwindigkeits- und ccc
ccc temperaturabhaengig berechnet. ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
a=props(3)
press0=props(4)
b=props(5)
ve0=props(6)
c=props(7)
temp0=props(8)
do k = 1, nBlock
disp=sqrt(dSlipFric(1,k)**2.0 +
* dSlipFric(2,k)**2.0 +
* dSlipFric(3,k)**2.0)
if (fNormal(k) .ne. 0.0 .and.
* rData(iDTimCur) .ne. 0.0 .and.
* disp .ne. 0.0) then
press=fNormal(k)/areaCont(k)
ve=disp/rData(iDTimCur)
temp=(tempMst(k)+tempSlv(k))/2.0
KontMu(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0)
* *(ve/ve0)**(b-1.0)
* *(temp/temp0)**(c-1.0)
WeightedMu(jConSlvUid(1,k))=WeightedMu(jConSlvUid(1,k))+
* KontMu(jConSlvUid(1,k))*
* rData(3)
Anhang 33
TimeSum(jConSlvUid(1,k))=TimeSum(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMu(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
KontMu(jConSlvUid(1,k)) = 0.0
end if
KurentNodeInc(jConSlvUid(1,k))=jFlags(2)
KurentNodeStep(jConSlvUid(1,k))=jFlags(1)
end do
end if
return
end
c
SUBROUTINE VUFIELD(FIELD, NBLOCK, NFIELD, KFIELD, NCOMP,
1 KSTEP, JFLAGS, JNODEID, TIME,
2 COORDS, U, V, A)
C
INCLUDE 'VABA_PARAM.INC'
C indices for the time array TIME
PARAMETER( i_ufld_Current = 1,
* i_ufld_Increment = 2,
* i_ufld_Period = 3,
* i_ufld_Total = 4 )
C indices for the coordinate array COORDS
PARAMETER( i_ufld_CoordX = 1,
* i_ufld_CoordY = 2,
* i_ufld_CoordZ = 3 )
C indices for the displacement array U
PARAMETER( i_ufld_SpaDisplX = 1,
* i_ufld_SpaDisplY = 2,
* i_ufld_SpaDisplZ = 3,
* i_ufld_RotDisplX = 4,
* i_ufld_RotDisplY = 5,
* i_ufld_RotDisplZ = 6,
* i_ufld_AcoPress = 7,
* i_ufld_Temp = 8 )
C indices for the velocity array V
PARAMETER( i_ufld_SpaVelX = 1,
* i_ufld_SpaVelY = 2,
* i_ufld_SpaVelZ = 3,
* i_ufld_RotVelX = 4,
* i_ufld_RotVelY = 5,
* i_ufld_RotVelZ = 6,
* i_ufld_DAcoPress = 7,
* i_ufld_DTemp = 8 )
C indices for the acceleration array A
PARAMETER( i_ufld_SpaAccelX = 1,
* i_ufld_SpaAccelY = 2,
* i_ufld_SpaAccelZ = 3,
* i_ufld_RotAccelX = 4,
* i_ufld_RotAccelY = 5,
* i_ufld_RotAccelZ = 6,
* i_ufld_DDAcoPress = 7,
* i_ufld_DDTemp = 8 )
Anhang 34
C indices for JFLAGS
PARAMETER( i_ufld_kInc = 1,
* i_ufld_kPass = 2 )
C
DIMENSION FIELD(NBLOCK,NCOMP,NFIELD)
DIMENSION JFLAGS(2), JNODEID(NBLOCK), TIME(4),
* COORDS(3,NBLOCK)
DIMENSION U(8,NBLOCK), V(8,NBLOCK), A(8,NBLOCK)
c
ccc-----------------------------------------------------------------ccc
ccc Benutzer Code zur definition der VUFIELD Subroutine ccc
ccc-----------------------------------------------------------------ccc
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc \/ Variablen Deklarationen \/ ccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
ccc ProbSize bestimmt die Groeße der common Blocks.
ccc ProbSize muss mindestens 2* Anz. der Knoten im zu untersuchenden
ccc Node Set sein (Bspw. 2* Knoten Anz. des Blechwerkstücks).
c
integer ProbSize
parameter (ProbSize=1100000)
c
ccc KontMu ist das Array zur Uebergabe der Reibungskoeffizienten aus
ccc der VFRICTION in die VUFIELD
c
real KontMu(ProbSize)
c
ccc WeightedMu dient als Array zur Berechnung eines Durchschnitts-
ccc wertes der Reibungskoeffizienten
c
double precision WeightedMu(ProbSize)
c
ccc Zur Berechnung des Mittelwertes wird die Summe ueber die Zeit-
ccc inkremente genutzt, bei der ein Reibungskoeffizient ungleich
ccc 0 berechnet wurde
c
real TimeSum(ProbSize)
c
ccc KurentNodeInc und KurentNodeStep zeigen von welchem Inkrement und
ccc Schritt die Reibungskoeffizienten in KuntMu stammen.
ccc KurentInc und KurentStep beinhalten die aktuelle Inkrement Nr.
ccc und Step Nr..
c
integer KurentNodeInc(ProbSize), KurentNodeStep(ProbSize),
* KurentInc, KurentStep
c
ccc KontMu, KurentNodeInc, KurentNodeStep, KurentInc, KurentStep und
ccc WeightedMu werden als common Blocks definiert.
c
common /KontaktReal/ KontMu, WeightedMu, TimeSum
common /KTime/ KurentNodeInc, KurentInc, KurentNodeStep,
* KurentStep
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc /\ Ende der Variablen Deklarationen /\ ccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
ccc Initialisierung der Array KontMu, WeightedMu und KurentNodeInc
ccc zur Zeit des ersten Inkrements der Rechnung
c
if (TIME(i_ufld_Total) .eq. 0.0) then
Anhang 35
do k = 1, ProbSize
KontMu(k) = 0.0
KurentNodeInc(k) = 0
KurentNodeStep(k) = 1
WeightedMu(k) = 0.0
TimeSum(k) = 0.0
end do
end if
c
ccc Update der Field Variable zur Erzeugung des Konturplots
c
if (TIME(i_ufld_Total) .ne. 0.0) then
c
ccc In der ersten Field Variable werden die Reibungskoeffizienten
ccc mit Historie uebergeben
c
do k = 1,NBLOCK
FIELD(k, 1, 1) = KontMu(JNODEID(k))
end do
c
ccc In der zweiten Field Variable werden die Reibungskoeffizienten
ccc ohne Historie uebertragen
c
do k = 1,NBLOCK
if (KurentNodeInc(JNODEID(k)) .eq. KurentInc .and.
* KurentNodeStep(JNODEID(k)) .eq. KurentStep ) then
FIELD(k, 1, 2) = KontMu(JNODEID(k))
else
FIELD(k, 1, 2) = 0.0
end if
end do
c
ccc in der dritten Field Variable werden die maximalen
ccc Reibungskoeffizienten gespeichert
c
do k= 1,NBLOCK
if (FIELD(k, 1, 3) .lt. KontMu(JNODEID(k))) then
FIELD(k, 1, 3)=KontMu(JNODEID(k))
end if
end do
c
ccc in der vierten Field Variable werden die Reibungs-
ccc koeffizienten nur über den zeitbereich gemittelt, in den sie
ccc ungleich 0 sind, also Kontakt besteht
c
do k= 1,NBLOCK
if (TimeSum(JNODEID(k)) .ne. 0.0) then
FIELD(k, 1, 4)=WeightedMu(JNODEID(k))/TimeSum(JNODEID(k))
else
FIELD(k, 1, 4)=0.0
end if
end do
end if
RETURN
END
A.2 Subroutinen für die Verwendung mit Shell Elementen
c
c User subroutine VFRICTION to define friction forces
c
subroutine vfriction (
c Write only -
Anhang 36
* fTangential,
c Read/Write -
* state,
c Read only -
* nBlock, nBlockAnal, nBlockEdge,
* nNodState, nNodSlv, nNodMst,
* nFricDir, nDir,
* nStates, nProps, nTemp, nFields,
* jFlags, rData,
* surfInt, surfSlv, surfMst,
* jConSlvUid, jConMstUid, props,
* dSlipFric, fStickForce, fTangPrev, fNormal,
* areaCont, dircosN, dircosSl,
* shapeSlv, shapeMst,
* coordSlv, coordMst,
* velSlv, velMst,
* tempSlv, tempMst,
* fieldSlv, fieldMst )
c
c Array dimensioning variables:
c
c nBlockAnal = nBlock (non-analytical-rigid master surface)
c nBlockAnal = 1 (analytical rigid master surface)
c nBlockEdge = 1 (non-edge-type slave surface)
c nBlockEdge = nBlock (edge-type slave surface)
c nNodState = 1 (node-type slave surface)
c nNodState = 4 (edge-type slave surface)
c nNodSlv = 1 (node-type slave surface)
c nNodSlv = 2 (edge-type slave surface)
c nNodMst = 4 (facet-type master surface)
c nNodMst = 2 (edge-type master surface)
c nNodMst = 1 (analytical rigid master surface)
c
c Surface names surfSlv and surfMst are not available for
c general contact (set to blank).
c
include 'vaba_param.inc'
dimension fTangential(nFricDir,nBlock),
* state(nStates,nNodState,nBlock),
* jConSlvUid(nNodSlv,nBlock),
* jConMstUid(nNodMst,nBlockAnal),
* props(nProps),
* dSlipFric(nDir,nBlock),
* fStickForce(nBlock),
* fTangPrev(nDir,nBlock),
* fNormal(nBlock),
* areaCont(nBlock),
* dircosN(nDir,nBlock),
* dircosSl(nDir,nBlock),
* shapeSlv(nNodSlv,nBlockEdge),
* shapeMst(nNodMst,nBlockAnal),
* coordSlv(nDir,nNodSlv,nBlock),
* coordMst(nDir,nNodMst,nBlockAnal),
* velSlv(nDir,nNodSlv,nBlock),
* velMst(nDir,nNodMst,nBlockAnal),
* tempSlv(nBlock),
* tempMst(nBlockAnal),
* fieldSlv(nFields,nBlock),
* fieldMst(nFields,nBlockAnal)
c
parameter( iKStep = 1,
* iKInc = 2,
* iLConType = 3,
Anhang 37
* nFlags = 3 )
parameter( iTimStep = 1,
* iTimGlb = 2,
* iDTimCur = 3,
* iFrictionWork = 4,
* nData = 4 )
c
dimension jFlags(nFlags), rData(nData)
character*80 surfInt, surfSlv, surfMst
c
ccc-----------------------------------------------------------------ccc
ccc Benutzer Code zur Definition der VFRICTION Subroutine ccc
ccc-----------------------------------------------------------------ccc
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc \/ Variablen Deklarationen \/ ccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
ccc ProbSize bestimmt die Groeße der common Blocks.
ccc ProbSize muss mindestens 2* Anz. der Knoten im zu untersuchenden
ccc Node Set sein (Bspw. 2* Knoten Anz. des Blechstücks).
c
integer ProbSize
parameter (ProbSize=1100000)
c
ccc KontMu ist das Array zur Uebergabe der Reibungskoeffizienten aus
ccc der VFRICTION in die VUFIELD
c
real KontMuO(ProbSize), KontMuU(ProbSize)
c
ccc WeightedMu dient als Array zur Berechnung eines Durchschnitts-
ccc wertes der Reibungskoeffizienten
c
double precision WeightedMuO(ProbSize), WeightedMuU(ProbSize)
c
ccc Zur Berechnung des Mittelwertes wird die Summe ueber die Zeit-
ccc inkremente genutzt, bei der ein Reibungskoeffizient ungleich
ccc 0 berechnet wurde
c
real TimeSumO(ProbSize), TimeSumU(ProbSize)
c
ccc Einige lokale Variablen zur Berechnung und Definition der
ccc physikalischen Parameter: Druck, Temperatur,
ccc Exponenten: a, b, c, etc.
c
real temp0, ve0, press0, mu0, a, b, c, press, ve, temp, ft, disp
c
ccc KurentNodeInc und KurentNodeStep zeigen von welchem Inkrement und
ccc Schritt die Reibungskoeffizienten in KuntMu stammen.
ccc KurentInc und KurentStep beinhalten die aktuelle Inkrement Nr.
ccc und Step Nr..
c
integer KurentNodeIncO(ProbSize), KurentNodeStepO(ProbSize),
* KurentInc, KurentStep,
* KurentNodeIncU(ProbSize), KurentNodeStepU(ProbSize)
c
ccc Definition der lokalen Variable Abh, dient der Bestimmung welche
ccc Berechnungsvorschrift für den Reibungskoeffizient genutzt werden
ccc soll.
integer Abh
c
ccc KontMu, KurentNodeInc, KurentNodeStep, KurentInc, KurentStep und
ccc HistoryOpt werden als common Blocks definiert.
Anhang 38
c
common /KontaktReal/ KontMuO, WeightedMuO, TimeSumO, KontMuU,
* WeightedMuU, TimeSumU
common /KTime/ KurentNodeIncO, KurentInc, KurentNodeStepO,
* KurentStep, KurentNodeIncU, KurentNodeStepU
c
ccc save befehl stellt sicher, dass die common Blocks nach verlassen
ccc der Subroutine erhalten bleibt.
c
save /KontaktReal/
save /KTime/
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc /\ Ende der Variablen Deklarationen /\ ccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
ccc Die aktuelle Inkrement Nr. und Schritt Nr. wird Upgedatet
c
KurentInc=jFlags(2)
KurentStep=jFlags(1)
c
ccc Die gewaehlte Berechnungsvorschrift für die
ccc Reibungskoeffizientenwird werden aus der Abaqus Inputfile an
ccc Abh in Form eines numerischen Wertes uebergeben.
c
Abh=props(1)
c
if ( Abh .eq. 1) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird druckabheangig berechnet ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
a=props(3)
press0=props(4)
do k = 1, nBlock
if (fNormal(k) .ne. 0.0) then
press=fNormal(k)/areaCont(k)
if (surfInt .eq. 'CONTECT-O') then
KontMuO(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0)
WeightedMuO(jConSlvUid(1,k))=
* WeightedMuO(jConSlvUid(1,k))+
* KontMuO(jConSlvUid(1,k))*
* rData(3)
TimeSumO(jConSlvUid(1,k))=TimeSumO(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuO(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
else
KontMuU(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0)
WeightedMuU(jConSlvUid(1,k))=
* WeightedMuU(jConSlvUid(1,k))+
* KontMuU(jConSlvUid(1,k))*
* rData(3)
TimeSumU(jConSlvUid(1,k))=TimeSumU(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuU(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
end if
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
if (surfInt .eq. 'CONTECT-O') then
KontMuO(jConSlvUid(1,k)) = 0.0
Anhang 39
else
KontMuU(jConSlvUid(1,k)) = 0.0
end if
end if
if (surfInt .eq. 'CONTECT-O') then
KurentNodeIncO(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepO(jConSlvUid(1,k))=jFlags(1)
else
KurentNodeIncU(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepU(jConSlvUid(1,k))=jFlags(1)
end if
end do
else if( Abh .eq. 2) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird rel. geschwindigkeitsabhaengig ccc
ccc berechnet ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
b=props(3)
ve0=props(4)
do k = 1, nBlock
disp=sqrt(dSlipFric(1,k)**2.0 +
* dSlipFric(2,k)**2.0 +
* dSlipFric(3,k)**2.0)
if (fNormal(k) .ne. 0.0 .and.
* rData(iDTimCur) .ne. 0.0 .and.
* disp .ne. 0.0) then
ve=disp/rData(iDTimCur)
if (surfInt .eq. 'CONTECT-O') then
KontMuO(jConSlvUid(1,k))=mu0*(ve/ve0)**(b-1.0)
WeightedMuO(jConSlvUid(1,k))=
* WeightedMuO(jConSlvUid(1,k))+
* KontMuO(jConSlvUid(1,k))*
* rData(3)
TimeSumO(jConSlvUid(1,k))=TimeSumO(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuO(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
else
KontMuO(jConSlvUid(1,k))=mu0*(ve/ve0)**(b-1.0)
WeightedMuO(jConSlvUid(1,k))=
* WeightedMuO(jConSlvUid(1,k))+
* KontMuO(jConSlvUid(1,k))*
* rData(3)
TimeSumO(jConSlvUid(1,k))=TimeSumO(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuO(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
end if
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
if (surfInt .eq. 'CONTECT-O') then
KontMuO(jConSlvUid(1,k)) = 0.0
else
KontMuU(jConSlvUid(1,k)) = 0.0
end if
end if
if (surfInt .eq. 'CONTECT-O') then
KurentNodeIncO(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepO(jConSlvUid(1,k))=jFlags(1)
else
KurentNodeIncU(jConSlvUid(1,k))=jFlags(2)
Anhang 40
KurentNodeStepU(jConSlvUid(1,k))=jFlags(1)
end if
end do
else if( Abh .eq. 3) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird temperaturabhaengig berechnet. ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
c=props(3)
temp0=props(4)
do k = 1, nBlock
temp=(tempMst(k)+tempSlv(k))/2.0
if (fNormal(k) .ne. 0.0 .and.
* temp .ne. 0.0) then
if (surfInt .eq. 'CONTECT-O') then
KontMuO(jConSlvUid(1,k))=mu0*(temp/temp0)**(c-1.0)
WeightedMuO(jConSlvUid(1,k))=
* WeightedMuO(jConSlvUid(1,k))+
* KontMuO(jConSlvUid(1,k))*
* rData(3)
TimeSumO(jConSlvUid(1,k))=TimeSumO(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuO(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
else
KontMuU(jConSlvUid(1,k))=mu0*(temp/temp0)**(c-1.0)
WeightedMuU(jConSlvUid(1,k))=
* WeightedMuU(jConSlvUid(1,k))+
* KontMuU(jConSlvUid(1,k))*
* rData(3)
TimeSumU(jConSlvUid(1,k))=TimeSumU(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuU(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
end if
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
if (surfInt .eq. 'CONTECT-O') then
KontMuO(jConSlvUid(1,k)) = 0.0
else
KontMuU(jConSlvUid(1,k)) = 0.0
end if
end if
if (surfInt .eq. 'CONTECT-O') then
KurentNodeIncO(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepO(jConSlvUid(1,k))=jFlags(1)
else
KurentNodeIncU(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepU(jConSlvUid(1,k))=jFlags(1)
end if
end do
else if( Abh .eq. 4) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird druck- und rel ccc
ccc geschwindigkeitsabhaengig berechnet. ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
a=props(3)
press0=props(4)
b=props(5)
ve0=props(6)
do k = 1, nBlock
Anhang 41
disp=sqrt(dSlipFric(1,k)**2.0 +
* dSlipFric(2,k)**2.0 +
* dSlipFric(3,k)**2.0)
if (fNormal(k) .ne. 0.0 .and.
* areaCont(k) .ne. 0.0 .and.
* rData(iDTimCur) .ne. 0.0 .and.
* disp .ne. 0.0) then
if (surfInt .eq. 'CONTECT-O') then
press=fNormal(k)/areaCont(k)
ve=disp/rData(iDTimCur)
KontMuO(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0)
* *(ve/ve0)**(b-1.0)
WeightedMuO(jConSlvUid(1,k))=
* WeightedMuO(jConSlvUid(1,k))+
* KontMuO(jConSlvUid(1,k))*
* rData(3)
TimeSumO(jConSlvUid(1,k))=TimeSumO(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuO(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
print *,'CONTECT-O'
else
press=fNormal(k)/areaCont(k)
ve=disp/rData(iDTimCur)
KontMuU(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0)
* *(ve/ve0)**(b-1.0)
WeightedMuU(jConSlvUid(1,k))=
* WeightedMuU(jConSlvUid(1,k))+
* KontMuU(jConSlvUid(1,k))*
* rData(3)
TimeSumU(jConSlvUid(1,k))=TimeSumU(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuU(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
print *,'CONTECT-U'
end if
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
if (surfInt .eq. 'CONTECT-O') then
KontMuO(jConSlvUid(1,k)) = 0.0
else
KontMuU(jConSlvUid(1,k)) = 0.0
end if
end if
if (surfInt .eq. 'CONTECT-O') then
KurentNodeIncO(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepO(jConSlvUid(1,k))=jFlags(1)
else
KurentNodeIncU(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepU(jConSlvUid(1,k))=jFlags(1)
end if
end do
else if(Abh .eq. 5) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird geschwindigkeits- und ccc
ccc temperaturabhaengig berechnet. ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
a=props(3)
press0=props(4)
c=props(5)
temp0=props(6)
Anhang 42
do k = 1, nBlock
temp=(tempMst(k)+tempSlv(k))/2.0
if (fNormal(k) .ne. 0.0 .and.
* temp .ne. 0.0) then
if (surfInt .eq. 'CONTECT-O') then
press=fNormal(k)/areaCont(k)
KontMuO(jConSlvUid(1,k))=
* mu0*(press/press0)**(a-1.0)
* *(temp/temp0)**(c-1.0)
WeightedMuO(jConSlvUid(1,k))=
* WeightedMuO(jConSlvUid(1,k))+
* KontMuO(jConSlvUid(1,k))*
* rData(3)
TimeSumO(jConSlvUid(1,k))=TimeSumO(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuO(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
else
press=fNormal(k)/areaCont(k)
KontMuU(jConSlvUid(1,k))=
* mu0*(press/press0)**(a-1.0)
* *(temp/temp0)**(c-1.0)
WeightedMuU(jConSlvUid(1,k))=
* WeightedMuU(jConSlvUid(1,k))+
* KontMuU(jConSlvUid(1,k))*
* rData(3)
TimeSumU(jConSlvUid(1,k))=TimeSumU(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuU(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
end if
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
if (surfInt .eq. 'CONTECT-O') then
KontMuO(jConSlvUid(1,k)) = 0.0
else
KontMuU(jConSlvUid(1,k)) = 0.0
end if
end if
if (surfInt .eq. 'CONTECT-O') then
KurentNodeIncO(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepO(jConSlvUid(1,k))=jFlags(1)
else
KurentNodeIncU(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepU(jConSlvUid(1,k))=jFlags(1)
end if
end do
else if(Abh .eq. 6) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird geschwindigkeits- und ccc
ccc temperaturabhaengig berechnet. ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
b=props(3)
ve0=props(4)
c=props(5)
temp0=props(6)
do k = 1, nBlock
disp=sqrt(dSlipFric(1,k)**2.0 +
* dSlipFric(2,k)**2.0 +
* dSlipFric(3,k)**2.0)
temp=(tempMst(k)+tempSlv(k))/2.0
Anhang 43
if (fNormal(k) .ne. 0.0 .and.
* rData(iDTimCur) .ne. 0.0 .and.
* disp .ne. 0.0 .and.
* temp .ne. 0.0) then
if (surfInt .eq. 'CONTECT-O') then
ve=disp/rData(iDTimCur)
KontMuO(jConSlvUid(1,k))=mu0*(ve/ve0)**(b-1.0)
* *(temp/temp0)**(c-1.0)
WeightedMuO(jConSlvUid(1,k))=
* WeightedMuO(jConSlvUid(1,k))+
* KontMuO(jConSlvUid(1,k))*
* rData(3)
TimeSumO(jConSlvUid(1,k))=TimeSumO(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuO(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
else
ve=disp/rData(iDTimCur)
KontMuU(jConSlvUid(1,k))=mu0*(ve/ve0)**(b-1.0)
* *(temp/temp0)**(c-1.0)
WeightedMuU(jConSlvUid(1,k))=
* WeightedMuU(jConSlvUid(1,k))+
* KontMuU(jConSlvUid(1,k))*
* rData(3)
TimeSumU(jConSlvUid(1,k))=TimeSumU(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuU(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
end if
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
if (surfInt .eq. 'CONTECT-O') then
KontMuO(jConSlvUid(1,k)) = 0.0
else
KontMuU(jConSlvUid(1,k)) = 0.0
end if
end if
if (surfInt .eq. 'CONTECT-O') then
KurentNodeIncO(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepO(jConSlvUid(1,k))=jFlags(1)
else
KurentNodeIncU(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepU(jConSlvUid(1,k))=jFlags(1)
end if
end do
else if(Abh .eq. 7) then
ccc-----------------------------------------------------------------ccc
ccc Reibungskoeffizient wird druck-, rel. geschwindigkeits- und ccc
ccc temperaturabhaengig berechnet. ccc
ccc-----------------------------------------------------------------ccc
mu0=props(2)
a=props(3)
press0=props(4)
b=props(5)
ve0=props(6)
c=props(7)
temp0=props(8)
do k = 1, nBlock
disp=sqrt(dSlipFric(1,k)**2.0 +
* dSlipFric(2,k)**2.0 +
* dSlipFric(3,k)**2.0)
temp=(tempMst(k)+tempSlv(k))/2.0
Anhang 44
if (fNormal(k) .ne. 0.0 .and.
* rData(iDTimCur) .ne. 0.0 .and.
* disp .ne. 0.0 .and.
* temp .ne. 0.0) then
if (surfInt .eq. 'CONTECT-O') then
press=fNormal(k)/areaCont(k)
ve=disp/rData(iDTimCur)
KontMuO(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0)
* *(ve/ve0)**(b-1.0)
* *(temp/temp0)**(c-1.0)
WeightedMuO(jConSlvUid(1,k))=
* WeightedMuO(jConSlvUid(1,k))+
* KontMuO(jConSlvUid(1,k))*
* rData(3)
TimeSumO(jConSlvUid(1,k))=TimeSumO(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuO(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
else
press=fNormal(k)/areaCont(k)
ve=disp/rData(iDTimCur)
KontMuU(jConSlvUid(1,k))=mu0*(press/press0)**(a-1.0)
* *(ve/ve0)**(b-1.0)
* *(temp/temp0)**(c-1.0)
WeightedMuU(jConSlvUid(1,k))=
* WeightedMuU(jConSlvUid(1,k))+
* KontMuU(jConSlvUid(1,k))*
* rData(3)
TimeSumU(jConSlvUid(1,k))=TimeSumU(jConSlvUid(1,k))+
* rData(3)
ft = min ( KontMuU(jConSlvUid(1,k))*fNormal(k),
* fStickForce(k) )
end if
fTangential(1,k) = -ft
fTangential(2,k) = 0.0
else
if (surfInt .eq. 'CONTECT-O') then
KontMuO(jConSlvUid(1,k)) = 0.0
else
KontMuU(jConSlvUid(1,k)) = 0.0
end if
end if
if (surfInt .eq. 'CONTECT-O') then
KurentNodeIncO(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepO(jConSlvUid(1,k))=jFlags(1)
else
KurentNodeIncU(jConSlvUid(1,k))=jFlags(2)
KurentNodeStepU(jConSlvUid(1,k))=jFlags(1)
end if
end do
end if
return
end
c
SUBROUTINE VUFIELD(FIELD, NBLOCK, NFIELD, KFIELD, NCOMP,
1 KSTEP, JFLAGS, JNODEID, TIME,
2 COORDS, U, V, A)
C
INCLUDE 'VABA_PARAM.INC'
C indices for the time array TIME
PARAMETER( i_ufld_Current = 1,
Anhang 45
* i_ufld_Increment = 2,
* i_ufld_Period = 3,
* i_ufld_Total = 4 )
C indices for the coordinate array COORDS
PARAMETER( i_ufld_CoordX = 1,
* i_ufld_CoordY = 2,
* i_ufld_CoordZ = 3 )
C indices for the displacement array U
PARAMETER( i_ufld_SpaDisplX = 1,
* i_ufld_SpaDisplY = 2,
* i_ufld_SpaDisplZ = 3,
* i_ufld_RotDisplX = 4,
* i_ufld_RotDisplY = 5,
* i_ufld_RotDisplZ = 6,
* i_ufld_AcoPress = 7,
* i_ufld_Temp = 8 )
C indices for the velocity array V
PARAMETER( i_ufld_SpaVelX = 1,
* i_ufld_SpaVelY = 2,
* i_ufld_SpaVelZ = 3,
* i_ufld_RotVelX = 4,
* i_ufld_RotVelY = 5,
* i_ufld_RotVelZ = 6,
* i_ufld_DAcoPress = 7,
* i_ufld_DTemp = 8 )
C indices for the acceleration array A
PARAMETER( i_ufld_SpaAccelX = 1,
* i_ufld_SpaAccelY = 2,
* i_ufld_SpaAccelZ = 3,
* i_ufld_RotAccelX = 4,
* i_ufld_RotAccelY = 5,
* i_ufld_RotAccelZ = 6,
* i_ufld_DDAcoPress = 7,
* i_ufld_DDTemp = 8 )
C indices for JFLAGS
PARAMETER( i_ufld_kInc = 1,
* i_ufld_kPass = 2 )
C
DIMENSION FIELD(NBLOCK,NCOMP,NFIELD)
DIMENSION JFLAGS(2), JNODEID(NBLOCK), TIME(4),
* COORDS(3,NBLOCK)
DIMENSION U(8,NBLOCK), V(8,NBLOCK), A(8,NBLOCK)
c
ccc-----------------------------------------------------------------ccc
ccc Benutzer Code zur definition der VUFIELD Subroutine ccc
ccc-----------------------------------------------------------------ccc
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc \/ Variablen Deklarationen \/ ccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
ccc ProbSize bestimmt die Groeße der common Blocks.
ccc ProbSize muss mindestens 2* Anz. der Knoten im zu untersuchenden
ccc Node Set sein (Bspw. 2* Knoten Anz. des Blechwerkstücks).
c
integer ProbSize
parameter (ProbSize=1100000)
c
Anhang 46
ccc KontMu ist das Array zur Uebergabe der Reibungskoeffizienten aus
ccc der VFRICTION in die VUFIELD
c
real KontMuO(ProbSize), KontMuU(ProbSize)
c
ccc WeightedMu dient als Array zur Berechnung eines Durchschnitts-
ccc wertes der Reibungskoeffizienten
c
double precision WeightedMuO(ProbSize), WeightedMuU(ProbSize)
c
ccc Zur Berechnung des Mittelwertes wird die Summe ueber die Zeit-
ccc inkremente genutzt, bei der ein Reibungskoeffizient ungleich
ccc 0 berechnet wurde
c
real TimeSumO(ProbSize), TimeSumU(ProbSize)
c
ccc KurentNodeInc und KurentNodeStep zeigen von welchem Inkrement und
ccc Schritt die Reibungskoeffizienten in KuntMu stammen.
ccc KurentInc und KurentStep beinhalten die aktuelle Inkrement Nr.
ccc und Step Nr..
c
integer KurentNodeIncO(ProbSize), KurentNodeStepO(ProbSize),
* KurentInc, KurentStep,
* KurentNodeIncU(ProbSize), KurentNodeStepU(ProbSize)
c
ccc KontMu, KurentNodeInc, KurentNodeStep, KurentInc, KurentStep und
ccc WeightedMu werden als common Blocks definiert.
c
common /KontaktReal/ KontMuO, WeightedMuO, TimeSumO, KontMuU,
* WeightedMuU, TimeSumU
common /KTime/ KurentNodeIncO, KurentInc, KurentNodeStepO,
* KurentStep, KurentNodeIncU, KurentNodeStepU
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc /\ Ende der Variablen Deklarationen /\ ccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
ccc Initialisierung der Array KontMu, WeightedMu und KurentNodeInc
ccc zur Zeit des ersten Inkrements der Rechnung
c
if (TIME(i_ufld_Total) .eq. 0.0) then
do k = 1, ProbSize
KontMuO(k) = 0.0
KurentNodeIncO(k) = 0
KurentNodeStepO(k) = 1
WeightedMuO(k) = 0.0
TimeSumO(k) = 0.0
KontMuU(k) = 0.0
KurentNodeIncU(k) = 0
KurentNodeStepU(k) = 1
WeightedMuU(k) = 0.0
TimeSumU(k) = 0.0
end do
end if
c
ccc Update der Field Variable zur Erzeugung des Konturplots
c
if (TIME(i_ufld_Total) .ne. 0.0) then
c
ccc In der ersten Field Variable werden die Reibungskoeffizienten
ccc mit Historie uebergeben
c
do k = 1,NBLOCK
Anhang 47
FIELD(k, 1, 1) = KontMuO(JNODEID(k))
end do
c
ccc In der zweiten Field Variable werden die Reibungskoeffizienten
ccc ohne Historie uebertragen
c
do k = 1,NBLOCK
if (KurentNodeIncO(JNODEID(k)) .eq. KurentInc .and.
* KurentNodeStepO(JNODEID(k)) .eq. KurentStep ) then
FIELD(k, 1, 2) = KontMuO(JNODEID(k))
else
FIELD(k, 1, 2) = 0.0
end if
end do
c
ccc in der dritten Field Variable werden die maximalen
ccc Reibungskoeffizienten gespeichert
c
do k= 1,NBLOCK
if (FIELD(k, 1, 3) .lt. KontMuO(JNODEID(k))) then
FIELD(k, 1, 3)=KontMuO(JNODEID(k))
end if
end do
c
ccc in der vierten Field Variable werden die Reibungs-
ccc koeffizienten nur über den zeitbereich gemittelt, in den sie
ccc ungleich 0 sind, also Kontakt besteht
c
do k= 1,NBLOCK
if (TimeSumO(JNODEID(k)) .ne. 0.0) then
FIELD(k, 1, 4)=WeightedMuO(JNODEID(k))/
* TimeSumO(JNODEID(k))
else
FIELD(k, 1, 4)=0.0
end if
end do
c
ccc In der ersten Field Variable werden die Reibungskoeffizienten
ccc mit Historie uebergeben
c
do k = 1,NBLOCK
FIELD(k, 1, 5) = KontMuU(JNODEID(k))
end do
c
ccc In der zweiten Field Variable werden die Reibungskoeffizienten
ccc ohne Historie uebertragen
c
do k = 1,NBLOCK
if (KurentNodeIncU(JNODEID(k)) .eq. KurentInc .and.
* KurentNodeStepU(JNODEID(k)) .eq. KurentStep ) then
FIELD(k, 1, 6) = KontMuU(JNODEID(k))
else
FIELD(k, 1, 6) = 0.0
end if
end do
c
ccc in der dritten Field Variable werden die maximalen
ccc Reibungskoeffizienten gespeichert
c
do k= 1,NBLOCK
if (FIELD(k, 1, 7) .lt. KontMuU(JNODEID(k))) then
FIELD(k, 1, 7)=KontMuU(JNODEID(k))
end if
Anhang 48
end do
c
ccc in der vierten Field Variable werden die Reibungs-
ccc koeffizienten nur über den zeitbereich gemittelt, in den sie
ccc ungleich 0 sind, also Kontakt besteht
c
do k= 1,NBLOCK
if (TimeSumU(JNODEID(k)) .ne. 0.0) then
FIELD(k, 1, 8)=WeightedMuU(JNODEID(k))/
* TimeSumU(JNODEID(k))
else
FIELD(k, 1, 8)=0.0
end if
end do
end if
RETURN
END