<?page no="0"?> Bernd Acker Simulationstechnik Grundlagen und praktische Anwendungen 2. Auflage <?page no="1"?> Simulationstechnik <?page no="3"?> Bernd Acker Simulationstechnik -2., erweiterte Auflage <?page no="4"?> Bibliografische Information der Deutschen Nationalbibliothek Die Deutsche Nationalbibliothek verzeichnet diese Publikation in der Deutschen Nationalbibliografie; detaillierte bibliografische Daten sind im Internet über http: / / dnb. dnb.de abrufbar. © 2019 · expert verlag GmbH Dischingerweg 5 · D-72070 Tübingen Das Werk einschließlich aller seiner Teile ist urheberrechtlich geschützt. Jede Verwertung außerhalb der engen Grenzen des Urheberrechtsgesetzes ist ohne Zustimmung des Verlages unzulässig und strafbar. Das gilt insbesondere für Vervielfältigungen, Übersetzungen, Mikroverfilmungen und die Einspeicherung und Verarbeitung in elektronischen Systemen. Internet: www.expertverlag.de E-Mail: info@verlag.expert CPI books GmbH, Leck ISBN 978-3-8169-2999-4 <?page no="5"?> Zusammenfassung In der Einführung wird nach der Definition der Begriffe System und Modell auf Simulation im Rahmen von Systemstudien eingegangen. Anhand vieler gezeigter Beispiele aus Wissenschaft und Technik wird die große Bedeutung der Simulationstechnik in verschiedenen Anwendungsgebieten herausgestellt. Die Simulation zeitkontinuierlicher Systeme kann auf dem Digitalrechner durchgeführt werden. Dazu werden systemorientierte Sprachen - Beispiel ACSL oder Matlab Simulink - eingeführt. Beide Simulationsmethoden werden ausführlich behandelt. Dabei veranschaulichen Beispiele den Lehrinhalt. Nach einer Einführung zeitdiskreter Systeme, in der auch die Grundlagen der Statistik und Warteschlangentheorie vermittelt werden, wird die Simulation dieser Systemklasse mit der Sprache GPSS vorgestellt. Durch ein Beispiel, das stufenweise erweitert wird, kann die Struktur von GPSS in verständlicher Form dargestellt werden. Zur Vertiefung der Kenntnisse gibt es - Angaben über weiterführende Literatur - im Download-Bereich des Expert Verlags unter http: / / www.expertverlag.de/ (Kennwort: BAC_SIMTEC) die im Buch behandelten Matlab Simulink-Beispiele und weitere einfache Modelle als mdl-Files. Großen Dank schulde ich meiner Tochter Isabel für ihre wiederholte Lektüre und die vielen Kommentare und Anregungen. <?page no="7"?> Inhaltsverzeichnis 1 Einführung............................................................................................................ 1 1.1 System................................................................................................................... 1 1.2 Modell .................................................................................................................... 4 1.3 Systemsimulation..................................................................................................16 2 Simulation zeitkontinuierlicher Systeme...........................................................41 2.1 Einführung ............................................................................................................41 2.2 Simulation auf dem Digitalrechner ........................................................................44 2.3 ACSL Advanced Continuous Simulation Language...............................................49 2.4 Matlab Simulink.....................................................................................................62 3 Simulation zeitdiskreter Systeme ......................................................................91 3.1 Einführung ............................................................................................................91 3.2 Zufallsvariable.......................................................................................................94 3.3 Betriebsstrategien ...............................................................................................106 3.4 Beispiel Warte-Verlust-System............................................................................107 3.5 Simulation auf dem Digitalrechner (Beispiel GPSS) ............................................112 1.1.1 Systemstudie ......................................................................................................... 1 1.1.2 Der Systembegriff .................................................................................................. 2 1.1.3 Systembeschreibung.............................................................................................. 3 1.2.1 Physikalische Modelle............................................................................................ 5 1.2.2 Mathematische Modelle ........................................................................................10 1.2.3 Vorgehensweise bei der Modellbildung .................................................................12 1.3.1 Definition und Beispiele.........................................................................................16 1.3.2 Zusammenfassung ...............................................................................................37 2.2.1 Numerische Integration .........................................................................................44 2.2.2 Simulationssprachen für kontinuierliche Systeme .................................................47 2.3.1 Einführung ............................................................................................................49 2.3.2 Beispiel Gedämpftes nichtlineares Pendel ............................................................58 2.4.1 Beispiel Spannungsmessinstrument .....................................................................63 2.4.2 Beispiel Gedämpftes nichtlineares Pendel ............................................................69 2.4.3 Vertikaldynamik eines Zweimassen-Fahrzeugs.....................................................70 3.1.1 Unterscheidung zeitdiskreter und kontinuierlicher Systeme ..................................91 3.1.2 Konzept der diskreten Simulation..........................................................................92 3.1.3 Beispiel Fernsprechvermittlung .............................................................................92 3.2.1 Diskrete Zufallsvariable .........................................................................................94 3.2.2 Statistische Kenngrößen diskreter Zufallsvariablen ...............................................95 3.2.3 Stetige Zufallsvariable...........................................................................................98 3.2.4 Beispiel: Messung einer stetigen Zufallsvariablen .................................................99 3.2.5 Typische Verteilungsfunktionen ..........................................................................101 3.5.1 Allgemeines ........................................................................................................112 3.5.2 Erzeugung und Transport der Transaktionen durch das Modell ..........................112 3.5.3 Reihenfolge der Ereignisse .................................................................................114 3.5.4 Verzweigungen ...................................................................................................115 3.5.5 Einfaches Beispiel Fertigung und Prüfung von Teilen (Modell 1).........................116 3.5.6 Eine oder mehrere gleichartige Bedienstationen (Modelle 2 und 3).....................118 3.5.7 Block-Typen für weitere Messungen (Modell 4)...................................................122 3.5.8 Bedingte Verzweigung (Modell 5)........................................................................125 3.5.9 Simulation Modell 4.............................................................................................130 <?page no="8"?> 4 Anhang: ACSL...................................................................................................135 4.1 ACSL Model Definition Statements (Examples) ..................................................135 4.2 ACSL Runtime Executive Commands (Examples) ..............................................136 4.3 ACSL Simulationen mit dem nichtlinearen Pendel ..............................................137 5 Anhang: Aktive Federung mit Matlab Simulink ..............................................149 6 Weiterführende Literatur ..................................................................................157 7 Sachregister ......................................................................................................159 <?page no="9"?> Abbildungsverzeichnis Abbildung 1: Studie eines realen Systems....................................................................... 1 Abbildung 2: Von der Wirklichkeit zur Abbildung.............................................................. 4 Abbildung 3: Modelltypen ................................................................................................ 4 Abbildung 4: Beispiel für ein Analogie-Modell .................................................................. 5 Abbildung 5: Mechanischer Analogrechner ..................................................................... 7 Abbildung 6: Baumartiges Netzwerk ................................................................................ 7 Abbildung 7: Wasser-Analog-Computer........................................................................... 8 Abbildung 8: Mechanischer Analogrechner für statische Tragwerke................................ 9 Abbildung 9: Volkswirtschaftsmodell...............................................................................10 Abbildung 10: Blockdiagramm eines Fabrikationssystems ...............................................12 Abbildung 11: Nichtlineares Fahrzeugmodell....................................................................13 Abbildung 12: Nutzfahrzeugmodell ...................................................................................13 Abbildung 13: PKW Rohbau und Finite Element (FE) Ersatzsystem.................................14 Abbildung 14: NEWEUL Programmsystem, Generierung von Bewegungsgleichungen ....15 Abbildung 15: Beispiele dynamischer Modelle..................................................................16 Abbildung 16: Großer Windkanal in Stuttgart-Untertürkheim ............................................18 Abbildung 17: Geräuschversuche im Akustiklabor............................................................19 Abbildung 18: Crash-Simulation .......................................................................................19 Abbildung 19: Fahrdynamik-Simulation, Ergebnis-Animation............................................21 Abbildung 20: NFZ-Simulation, Ergebnis-Animation .........................................................21 Abbildung 21: Kosimulation Gesamtfahrzeugmodell.........................................................22 Abbildung 22: Dauerfestigkeit-Simulation mit Animation...................................................22 Abbildung 23: Ergonomisches Modell eines Menschen ....................................................24 Abbildung 24: Ergonomische Simulation ..........................................................................24 Abbildung 25: Ergonomie in der Fabrikmontage ...............................................................25 Abbildung 26: Humanoid-Dummies ..................................................................................27 Abbildung 27: Insassenklima-Simulation ..........................................................................27 Abbildung 28: Simulation von Eisbildung in der Ostsee ...................................................28 Abbildung 29: Hydropuls-Anlage ......................................................................................29 Abbildung 30: Regensimulation ........................................................................................30 Abbildung 31: Getriebe-Simulator.....................................................................................31 Abbildung 32: Simulation der Schwerelosigkeit ................................................................31 Abbildung 33: Fahrsimulator Gesamtansicht ....................................................................33 Abbildung 34: Fahrsimulator Blockschaltbild ....................................................................33 Abbildung 35: Fahrsimulator Beschleunigungssimulation .................................................34 Abbildung 36: Fahrsimulator Kontrollstation .....................................................................36 Abbildung 37: Fahrsimulator Bildsystem...........................................................................36 Abbildung 38: Ablauf einer Systemstudie mit Hilfe der Simulation ....................................39 Abbildung 39: Lösungen einfacher linearer Differenzialgleichungen .................................42 Abbildung 40: Radaufhängung .........................................................................................43 Abbildung 41: Bestimmtes Integral ...................................................................................44 Abbildung 42: Integration mit Rechteckformel ..................................................................45 Abbildung 43: Integration mit Trapezformel ......................................................................46 Abbildung 44: Integration mit Simpson-Regel...................................................................46 Abbildung 45: Ablauf einer Simulation nach CSSL-Standard............................................49 Abbildung 46: Ablaufdiagramm zur Erzeugung eines lauffähigen ACSL-Programms .......53 Abbildung 47: ACSL Simulationsablauf ............................................................................54 Abbildung 48: Gedämpftes nichtlineares Pendel ..............................................................58 Abbildung 49: Verhalten des nichtlinearen Pendels ..........................................................59 Abbildung 50: Spannungsmessgerät ................................................................................63 Abbildung 51: Wirkungsplan des Spannungsmessers ......................................................64 Abbildung 52: Matlab Command Window .........................................................................65 Abbildung 53: Aufbau des Spannungsmesser-Programms in MatlabSimulink ..................66 <?page no="10"?> Abbildung 54: Beispiel eines geöffneten Gain-Blocks .......................................................66 Abbildung 55: Simulink-Programm des Spannungsmessers.............................................67 Abbildung 56: Ergebnisdarstellung Sprungantwort Spannungsmesser .............................68 Abbildung 57: Nichtlineares Pendel mit Anfangsbedingung TH0=0.0 und THp0=10.0 ......69 Abbildung 58: Modell des Zweimassenschwingers ...........................................................71 Abbildung 59: Simulink-Modell des Zweimassenschwingers ............................................73 Abbildung 60: Create Subsystem .....................................................................................73 Abbildung 61: Simulink Modell des Zweimassenschwingers mit Subsystems...................74 Abbildung 62: Möglichkeiten des Signal-Routing (Goto From) .....................................75 Abbildung 63: Endgültige Fassung des Modells Zweimassenschwinger...........................75 Abbildung 64: Manöver 1: Reaktion mit schwacher Dämpfung d_A = 1400N/ (m/ s) ..........76 Abbildung 65: Manöver 1: Dynamische Radlaständerung bei Sprungantwort ...................77 Abbildung 66: Manöver 1: Sprungantwort mit starker Dämfung´d_A = 3000 N/ (m/ s) ........77 Abbildung 67: Manöver 2: Reaktion auf eine Zusatzlast F_ZL = 1000 N...........................78 Abbildung 68: Manöver 3: Reaktion auf Zusatzlast + 10 Hz-Anregung der Straße............78 Abbildung 69: Ersatz konstanter Werte durch wählbare Variablen ...................................80 Abbildung 70: m-File zur Zuweisung gewünschter Werte an alle Systemparameter .........80 Abbildung 71: Auszug aus dem Matlab Command Window..............................................81 Abbildung 72: Erzeugung der Ausgabegröße Radaufstandskraft .....................................82 Abbildung 73: Grafische Ausgabe der Radaufstandskraft.................................................82 Abbildung 74: Sprungantwort des Aufbaus mit unterschiedlicher Dämpfung ....................83 Abbildung 75: Möglichkeit von Radabheben .....................................................................84 Abbildung 76: Sprungantwort auf einen Sprung von -0,15 m ............................................85 Abbildung 77: Modellerweiterung für Dämpfer mit Zug- und Druckstufe ...........................85 Abbildung 78: Dämpfungskennlinie mit Zug- und Druckstufe............................................86 Abbildung 79: Sprungantwort des Aufbaus mit drei Dämpfungsvarianten.........................87 Abbildung 80: Komfort als Funktion der Aufbaubeschleunigung xpp_A ............................87 Abbildung 81: Absenkung des Aufbaus bei Sinus-Anregung ............................................88 Abbildung 82: Beispiel Fernsprechvermittlung ..................................................................92 Abbildung 83: Typischer Ablauf der Simulation ................................................................94 Abbildung 84: Beispiel einer diskreten Zufallsvariablen ....................................................95 Abbildung 85: Histogramm der absoluten Häufigkeiten ....................................................96 Abbildung 86: Histogramm der absoluten Summenhäufigkeit...........................................96 Abbildung 87: Wahrscheinlichkeitsdichte und Verteilungsfunktion ....................................98 Abbildung 88: Auswertung von Gesprächsdauern ..........................................................100 Abbildung 89: Exponential-Verteilung .............................................................................102 Abbildung 90: Gleichverteilung: Wahrscheinlichkeitsdichte und Verteilungsfunktion.......103 Abbildung 91: Erlang-k-Verteilung ..................................................................................104 Abbildung 92: Wahrscheinlichkeitsdichte der Normalverteilung ......................................105 Abbildung 93: Abfertigungsstrategien bei Warteschlangen .............................................106 Abbildung 94: Beispiel Frisörladen .................................................................................107 Abbildung 95: Zeitlicher Abstand der Ankünfte mit Exponentialverteilung.......................109 Abbildung 96: Verteilungsfunktion für ein getaktetes Ereignis.........................................109 Abbildung 97: Beispiel eines Simulationslaufs ................................................................110 Abbildung 98: Programmstruktur ....................................................................................111 Abbildung 99: GPSS-Blockschaltbild Modell 1................................................................116 Abbildung 100: GPSS-Blockschaltbild Modell 2................................................................119 Abbildung 101: GPSS-Blockschaltbild Modell 3................................................................121 Abbildung 102: GPSS-Blockschaltbild Modell 4................................................................124 Abbildung 103: Klasseneinteilung in Tabelle TABLE ........................................................125 Abbildung 104: GPSS-Blockschaltbild Modell 5................................................................126 Abbildung 105: Ersatzschaltbild des Zweimassenschwingers mit aktiver Federung .........150 Abbildung 106: Modell des Zweimassenschwinger mit aktiver Federung..........................151 Abbildung 107: Federbein des aktiven Systems ...............................................................152 Abbildung 108: Subsystem zur Berechnung der Federlänge ............................................152 Abbildung 109: Manöver 5: Mit der Zeit wachsende Zusatzlast F_ZL ..............................153 <?page no="11"?> Abbildung 110: Manöver 5: Kompensation der Federverkürzung .....................................153 Abbildung 111: Manöver 5: Absolutbewegung von Rad und Aufbau.................................154 Abbildung 112: Manöver 5: Druck im und Ölstrom in den Hydraulikzylinder .....................154 Abbildung 113: Manöver 5: Vom System benötigte hydraulische Leistung .......................155 <?page no="13"?> Tabellenverzeichnis Tabelle 1: Systembeispiele ............................................................................................ 3 Tabelle 2: Systemunterscheidung nach Zustandsänderung ..........................................38 Tabelle 3: Entwicklung der Software-Tools zur Simulation ............................................48 Tabelle 4: ACSL Systemvariable (Beispiele) .................................................................57 Tabelle 5: ACSL Schrittweitensteuerung .......................................................................57 Tabelle 6 Integrationsverfahren....................................................................................67 Tabelle 7: Absolute und relative Häufigkeit ...................................................................95 Tabelle 8: Absolute und relative Summenhäufigkeit ......................................................96 Tabelle 9: Klassenhäufigkeit von Gesprächsdauern......................................................99 Tabelle 10: Systemunterscheidung als Funktion möglicher Warteplätze .......................106 Tabelle 11: Festlegung der Betriebsstrategie des Beispiels ..........................................107 Tabelle 12: Kennung der Belegung einer Station des Beispiels ....................................107 Tabelle 13: Darstellung aller Speicherzellen des Beispiels............................................108 Tabelle 14: Simulationsablauf des Beispiels .................................................................110 Tabelle 15: Typen von Bedienstationen ........................................................................118 Tabelle 16: Klassenhäufigkeit der Durchlaufzeiten in Modell 4 ......................................133 Tabelle 17: ACSL Model Definition Statements (Examples) ..........................................135 Tabelle 18: ACSL Runtime Executive Commands (Examples)......................................136 <?page no="15"?> 1 Reales System Theoretische Untersuchung Analytische Berechnung Simulation Modell Experiment 1 Einführung Sehr vielfältig ist in allen Bereichen der • Wissenschaft, • Technik, • Planung, usw. die Aufgabe, ein System zu studieren, um dessen Eigenschaften und Verhalten zu verstehen, vorherzusagen oder zu verbessern. Simulationstechnik ist ein Hilfsmittel hierfür. Sie setzt ein geeignetes Modell voraus, welches das zu untersuchende System genügend genau beschreibt. Zum besseren Verständnis werden die Begriffe zunächst näher erläutert. 1.1 System 1.1.1 Systemstudie Abbildung 1 zeigt die prinzipiellen Möglichkeiten, wie die Eigenschaften eines Systems erarbeitet werden können. Man muss dabei unterscheiden, ob es sich um ein existierendes “reales“ System handelt oder um ein rein hypothetisches System. Am realen System können Experimente durchgeführt und so ein Einblick in das System gewonnen werden. Abbildung 1: Studie eines realen Systems <?page no="16"?> 2 Aber auch theoretische Untersuchungen können zum Ziel führen. Hat man sich ein Modell geschaffen, das das System genügend genau beschreibt, so kann man analytische Berechnungen durchführen oder über die Simulation Systemerkenntnisse gewinnen. Üblich ist auch eine parallele Vorgehensweise, wobei sich Experiment und theoretische Untersuchung gegenseitig beeinflussen und ergänzen. Oft sind Experimente am realen System gar nicht oder nur stark eingeschränkt durchführbar, z.B. in • der Wetterforschung, -vorhersage, (Mensch ohne Einfluss) • der Humanmedizin oder bei (Ethik) • Wirtschafts- und Produktionsprozessen. (Kosten) Deshalb ist es das Ziel, Systeme einer theoretischen Untersuchung (Berechnung, Simulation, Laborexperiment) zugänglich zu machen. Es ergibt sich dabei eine ganze Reihe von Vorteilen: Die Untersuchung existierender Systeme ist u. a. • kostengünstig Schulung im Flugsimulator • ungefährlich Fahren im Grenzbereich • reproduzierbar Unfallrekonstruktion, Reihenuntersuchungen • zeitsparend Lebensdauertest Die Untersuchung hypothetischer Systeme ist u. a. • kostengünstig Variantenrechnung • flexibel Parametervariation, -optimierung • entwicklungsgerecht frühzeitiger Funktionsnachweis 1.1.2 Der Systembegriff Unter einem System soll Folgendes verstanden werden: • Es besteht aus einer Sammlung von Elementen und Komponenten. • Die Attribute kennzeichnen die konstanten oder variablen Eigenschaften der Elemente des Systems. • Die Aktivitäten sind Vorgänge, die Änderungen dieser Eigenschaften hervorrufen. Dabei ist aus praktischen Gründen durch • Abgrenzung von System und Umwelt und • Detaillierungsgrad der Umfang eines zu definierenden Systems zu begrenzen. Abhängig von der Zielsetzung der Systembetrachtung muss es die interessierenden Elemente, die wesentlichen Attribute und die Aktivitäten beinhalten. <?page no="17"?> 3 In Tabelle 1 sind Beispiele zusammengestellt: System Elemente Attribute Aktivitäten Straßenverkehr Fahrzeuge Straßen Verkehrszeichen Geschwindigkeit Richtung Anzahl Fahren Parken Regeln Bank Kunden Konten Kassen Anzahl Kontostand Zinssatz Einnehmen Verzinsen Abfragen Telefonvermittlung Teilnehmer Gespräche Rufnummer Dauer Wählen Sprechen Supermarkt Kunden Lager Einkaufsliste Sortiment Bezahlen Verkaufen Elektrische Schaltung Leitung Schalter Energieversorgung Stromstärke Spannung Frequenz Schalten Verstärken Unterbrechen Tabelle 1: Systembeispiele 1.1.3 Systembeschreibung Der Zustand des Systems ist die Beschreibung aller Elemente, Attribute und Aktivitäten zu jedem Zeitpunkt. Verhalten und Eigenschaften eines Systems werden untersucht, indem Änderungen im Zustand des Systems verfolgt werden. Nach Art der zeitlichen Änderung werden unterschieden • (zeit-)kontinuierliche Systeme - Änderungen sind (überwiegend) kontinuierlich oder stetig, z.B. Verkehr, elektrische (analoge) Schaltungen, • (zeit-)diskrete Systeme - Änderungen erfolgen (überwiegend) zu diskreten Zeitpunkten, z.B. Bank, Telefonvermittlung, digitale Schaltungen. <?page no="18"?> 4 Modelle physikalisch mathematisch statisch dynamisch statisch dynamisch analytisch numerisch 1.2 Modell Die Darstellung in Abbildung 2 zeigt den prinzipiellen Weg von einem realen System zum Modell als vereinfachter Abbildung der Wirklichkeit. Abbildung 2: Von der Wirklichkeit zur Abbildung Das Modell ist ein physikalisches oder formales (mathematisches) Ersatzsystem, welches das betrachtete (wirkliche) System mehr oder weniger vereinfacht wiedergibt und gestattet, dessen charakteristisches Verhalten zu studieren. Es kann auf verschiedene Art und Weise erzeugt werden, siehe Abbildung 3: Abbildung 3: Modelltypen Die Modellbildung umfasst die beiden Schritte • Festlegung der Struktur • Versorgung mit Daten (Messung, Berechnung, Schätzung) Im Folgenden soll näher auf die verschiedenen Modellarten eingegangen werden. Vereinfachung Systembeschreibung System Umwelt Elemente Attribute Aktivitäten Modell Wirklichkeit Abbildung <?page no="19"?> 5 1.2.1 Physikalische Modelle Das Modell wird nach bekannten Ähnlichkeitsgesetzen aus der Wirklichkeit hergeleitet: • Die Attribute der Systemelemente werden durch physikalische Größen (z.B. elektrische Spannungen, Kräfte) dargestellt. • Die Aktivitäten werden durch physikalische Gesetze ersetzt. Man unterscheidet • statische physikalische Modelle Maßstabsmodelle (Fahrzeugmodell im Windkanal, Schiffsmodell) Anschauungsmodelle (Molekülmodell) Die Auswertung der statischen physikalischen Modelle erfolgt durch Messungen und Transformation (Ähnlichkeitsgesetze). • dynamische physikalische Modelle z.B. Nachbildung eines mechanischen Systems durch ein elektrisches physikalisches Modell, siehe Abbildung 4: Abbildung 4: Beispiel für ein Analogie-Modell Mechanisches System: 1-Massen-Schwinger (vereinfachte Fahrzeugfederung) Elektrisches physikalisches Modell: Reihenschwingkreis i u(t) Rad Masse M Feder K Dämpfer D Radachse Zeitabhängige Kraft Die Bewegung der Masse M wird durch die folgende Differenzialgleichung 2. Ordnung in x beschrieben: M x + D x + K x = F( t ) Aufbau R L C q = i dt Das Verhalten des elektrischen Schwingkreises wird durch die folgende Differenzialgleichung 2. Ordnung in q beschrieben: L q + R q + 1/ C q = u( t ) x <?page no="20"?> 6 Man erkennt, dass sich die beschreibenden Differenzialgleichungen ihrer Struktur nach genau entsprechen. Wählt man für die Koeffizienten vor den Summanden der Gleichungen gleiche Zahlenwerte, ist also z.B. die Dämpfungskonstante D des mechanischen Systems vom Wert gleich dem ohmschen Widerstand R des elektrischen Systems, usw., dann werden sich die beiden Systeme über die Zeit genau gleich verhalten, demnach analog ablaufen. Man kann sagen, dass es eine Analogie gibt zwischen den beiden Systemen: Die Bewegung x der Masse im mechanischen Modell entspricht im elektrischen Modell der Ladung q des Kondensators, und die Geschwindigkeit der Masse entspricht dem elektrischen Strom i. Dies bedeutet, dass das Verhalten des einen Systems mit dem jeweils anderen untersucht werden kann. Zur Auswertung der dynamischen physikalischen Modelle: Parametervariationen lassen sich leicht durchführen und qualitativ beobachten sowie quantitativ bestimmen (messen). „Analogiemaschinen“ nutzen also die mathematische Gleichheit zweier Systeme. Beispiele: • Der elektrische Analogrechner, der bis Ende der 60er Jahre des vorigen Jahrhunderts vielfach eingesetzt wurde, um das Verhalten dynamischer Systeme zu berechnen. Er wurde in erster Linie eingesetzt für Systeme, die durch gewöhnliche Differenzialgleichungen beschrieben werden können. Der Analogrechner wurde sukzessive ersetzt durch die in den 60er Jahren zunehmende Anzahl von Digitalrechnern, den ersten Computern der heutigen Zeit. • Hydraulische Analogiemodelle Einsatzgebiet: Strömungsvorgänge • Weitere, eher „akademische“ Beispiele von mechanischen Analogiemaschinen zeigen Abbildung 5 und Abbildung 6 auf Seite 7, sowie Abbildung 7 auf Seite 8. Anfang des 20. Jahrhunderts allerdings, bevor der Analogrechner entwickelt war, wurde eine mechanische Analogiemaschine zumindest einmal erfolgreich zur Minimierung von Lasten in der Baustatik eingesetzt, siehe Abbildung 8 auf Seite 9. Das Bild zeigt eine Rekonstruktion eines Modells, das der spanische Architekt Antonio Gaudi benutzte, um an ihm die optimale Gestaltung der Gewölbebögen der Kirche Sagrada Familia in Barcelona abzulesen. Das Bild steht auf dem Kopf, um den Eindruck eines Kirchengewölbes zu vermitteln. <?page no="21"?> 7 Tatsächlich hingen die Seile von der Decke nach unten; sie waren mit stetig verteilten kleinen Sandsäckchen bestückt, um die Gewichtsbelastung der Steinbögen zu simulieren. Abbildung 5: Mechanischer Analogrechner Abbildung 6: Baumartiges Netzwerk y 0 Länge l Durchbiegung y Länge l Durchbiegung y M g Durchbiegung y ~ l 4 bzw. l ~ y 1/ 4 Durchbiegung y ~ y 0 + const . M 3 y 0 <?page no="22"?> 8 Abbildung 7: Wasser-Analog-Computer y = a x 3 + b x² + c x + d Analoggerät zum Finden einer Ausgleichsgeraden für eine Reihe von Datenpunkten Analoggerät zum Finden der Lösung einer kubischen Gleichung Analoggerät zum Finden des Mittelwertes von Daten x³ x x² d <?page no="23"?> 9 Abbildung 8: Mechanischer Analogrechner für statische Tragwerke 1 1 Graefe, R. u.a. (1983): Ein verschollenes Modell und seine Rekonstruktion, in: Bauwelt- Sonderdruck: Konstruktionen fürs Museum - Museum für Technik, Bertelsmann Fachzeitschriften GmbH, Berlin. Massebehaftetes Seil: konstante Seilkraft, wenn die Form einer cosh-Kurve entspricht. Torbögen in cosh-Form somit auch optimal für gleichmäßige Belastung. Prinzip wurde angewandt bei der Auslegung der Kirche Sagrada Familia, Barcelona (Antonio Gaudi, ca. 1920). Das Modell war ca. 3 x 3 x 4 m groß, die Abbildung zeigt eine Rekonstruktion. <?page no="24"?> 10 1.2.2 Mathematische Modelle Die Modellbildung erfolgt durch mathematische Beschreibung des Systems: Elemente, Attribute mathematische Variable Aktivitäten mathematische Funktionen Mathematische Modelle sind die Voraussetzung dafür, dass Systeme einer (theoretischen) rechnerischen Untersuchung zugänglich gemacht werden. Man unterscheidet • statische mathematische Modelle Sie beschreiben das System in Gleichgewichtszuständen, nicht aber den (zeitabhängigen) Übergang. Ein Beispiel ist in Abbildung 9 aufgeführt: Abbildung 9: Volkswirtschaftsmodell Das einfache Modell der amerikanischen Volkswirtschaft beschreibt den Zusammenhang verschiedener Kenngrößen wie Investition, Steuern, nationales Einkommen, usw. Es besteht aus einem überbestimmten linearen Gleichungssystem und erlaubt somit durch geeignete Wahl der Modellparameter die Maximierung z.B. des nationalen Einkommens. Volkswirtschaftsmodell Verbrauch V = 20 + 0,7 ( E - S ) Investition I = 2 + 0,1 E Steuern S = 0,2 E Nationales Einkommen E = V + I + A Staatsausgaben A (alle Größen in Mrd. $) Gleichungssystem mit 4 Gleichungen und 5 Unbekannten • Eindeutige analytische Lösung, wenn 1 Variable vorgegeben • Numerische Näherungslösung erforderlich, wenn z.B. A und S so zu bestimmen sind, dass E maximal wird <?page no="25"?> 11 Allgemein: • analytische Lösung bei linearen Gleichungssystemen immer möglich bei nichtlinearen Gleichungssystemen teilweise durch Einsetzen: z.B. x = y² Lösungen: x 1 = 0 , y 1 = 0 y = - 2 x und x 2 = 1/ 4 , y 2 = - 1/ 2 • numerische Lösung bei Einsatz von Digitalrechnern für große lineare Gleichungssysteme - Näherungsverfahren für nichtlineare Gleichungssysteme - Optimierungsverfahren, Extremwertsuche - Lösung impliziter Gleichungen z.B. x = sin x bzw. x - sin x = Extr. = 0 • dynamische mathematische Modelle Hier werden die Änderungen der Systemattribute durch zeitabhängige Variablen beschrieben. So führt z.B. der mechanische Schwinger aus Kapitel 1.2.1 auf Seite 5 auf ) t ( F x K x D x M = + + , also eine lineare Differenzialgleichung 2. Ordnung mit konstanten Koeffizienten. Hier ist • die analytische Auswertung möglich, z.B. die Berechnung von Zeitantworten bzw. Frequenzgängen, • aber auch die numerische Lösung bei beliebigem F(t) sinnvoll. Allgemein: • Abhängig von der Größe und Komplexität (Nichtlineritäten, partielle Differenzialgleichungen) des Modells sind zur Auswertung meist numerische Näherungsverfahren erforderlich, u. a. Simulation. • Methoden und Verfahren sind unterschiedlich für kontinuierliche oder (diskrete) Abtastsysteme. <?page no="26"?> 12 1.2.3 Vorgehensweise bei der Modellbildung Die folgenden Richtlinien sind bei der Modellbildung hilfreich: • Blockstruktur Zerlegung des zu untersuchenden Systems in Teilsysteme (Blöcke). Durch Verknüpfung der Blöcke untereinander entsteht ein einfaches, übersichtliches Blockdiagramm. Als Beispiel hierfür ist in Abbildung 10 die Struktur eines Fertigungsbetriebes zu sehen. Abbildung 10: Blockdiagramm eines Fabrikationssystems • Detaillierung Die Berücksichtigung der für die Untersuchung wichtigen Systemeigenschaften führt zu einer Begrenzung der Komplexität. Dies führt zu größerer Transparenz besserer Interpretierbarkeit der Ergebnisse geringerem Rechenaufwand (Rechnerkapazität, Rechenzeit, …) Je nach Aufgabenstellung können für ein System unterschiedliche Modelle entstehen. Als Beispiel sei ein Flugzeug-System betrachtet: - Modellierung als Punktmasse zur Berechnung von Kurs und Kraftstoffverbrauch - Modellierung als Starrkörper zur Berechnung der Flugeigenschaften - Modellierung als elastischer Körper zur Berechnung der Struktureigenschaften (Schwingungen, Komfort, Geräusch, Festigkeit) Weitere Beispiele eines Starrkörpermodells sind in Abbildung 11 und in Abbildung 12 auf Seite 13 zu sehen, ein Modell einer elastischen PKW- Karosserie als Finite-Element-Modell (FE-Struktur) in Abbildung 13 auf Seite 14. Einkauf Fertige Produkte Fertigung Montage Versand Fertigungssteuerung Kundenaufträge Rohmaterial <?page no="27"?> 13 Abbildung 11: Nichtlineares Fahrzeugmodell Abbildung 12: Nutzfahrzeugmodell x 1 y 1 z 1 x 0 y 0 Inertialsystem z 0 y S z S x S Motor Achsmassen Motorlager z 6 Elastisch gelagertes Dämpferelement x 6 y 6 y R1 x R1 z R1 Karosserie Reifenfeder Elastisch gelagertes Dämpferelement Karosseriefestes Bezugs- System (Referenz1) Federelement Getriebe Mot Hinterräder (angetrieben) Hinterwagen (Körper 6) Hinterachse (Körper 7) Rad vr (Körper 5) Rad vl (Körper 4) Vorderachsträger (Körper 3) Kupplung Z Z Z Z A D A X B Y A A C B A Rad vl m Rad vr Rad hl Rad hr Vertikaldynamik Längsdynamik Vorderachswagen (Körper 1) Fahrerhaus (Körper 2) Vorderräder (nicht angetrieben) <?page no="28"?> 14 Abbildung 13: PKW Rohbau und Finite Element (FE) Ersatzsystem • Vereinfachung Die Möglichkeiten einer vereinfachten Darstellung sollten genutzt werden, u.a. - Linearisierung (Einsatz leistungsfähiger Methoden) - Beschreibung funktionaler Zusammenhänge durch gemessene Kennlinien und Kennfelder - Formale Methoden der Modell-Reduktion • Rechnergestützte Gleichungsgenerierung Rechnergestützte Verfahren zur Generierung mathematischer Modelle machen die Modellbildung einfacher, schneller und sicherer. Hierzu gehören u. a. - Pre-Prozessoren zur Netzgenerierung für Finite-Element-Verfahren - CAD-Verfahren Catia V4 und V5 - Symbolische Gleichungserstellung (z.B. NEWEUL) Abbildung 14 auf Seite 15 zeigt die zugehörige File-Struktur. Durch grafische Eingabe z.B. einzelner mechanischer Komponenten wie Massen, Federn, Dämpfer usw. wird ein mechanisches Modell aufgebaut, zu dessen Simulation ein Differenzialgleichungssystem benötigt wird. Diese Gleichungen müssen nun nicht von Hand mit Methoden von z.B. d’Alembert oder Lagrange aufgestellt werden, sondern werden mit diesem Tool automatisch erstellt. Rohbau PKW Karosserie PKW Finite Element - Ersatzsystem <?page no="29"?> 15 Abbildung 14: NEWEUL Programmsystem, Generierung von Bewegungsgleichungen Im diesem Kapitel wurde der Begriff Modell erläutert. Beispiele ausgewählter physikalischer und mathematischer Modelle wurden gezeigt und zum Schluss Empfehlungen gegeben, wie man bei der Modellbildung vorgeht. Wichtig ist, dass alle für die Fragestellung entscheidenden Effekte berücksichtigt werden und aus Gründen der Einfachheit und besseren Möglichkeit der Ergebnisinterpretation alles Andere weggelassen wird. Die Schwierigkeit besteht in der Unterscheidung der wichtigen von den unwichtigen Einflüssen. Dies erfordert ein wenig Feingefühl und Erfahrung bei der Modellbildung. Im nächsten Kapitel soll nun näher beschrieben werden, was man unter Systemsimulation versteht. Viele Beispiele zeigen das außerordentlich breite Spektrum an Einsatzfeldern der Simulation. Die Einsatzgebiete sind so unterschiedlich, dass nur eine Teilmenge aller Anwendungen exemplarisch dargestellt werden kann. Daher werden im Folgenden nur noch dynamische Systeme betrachtet. Interaktive Eingaben Eingabe File 09 File 08 Protokollfile Eingabefile File 36 Newton-Euler Gleichungen File 37 File 38 OUTPUT Bewegungs-Gleichungen FORTRAN-kompatibel NEWEUL Eingabe mit dem Editor <?page no="30"?> 16 1.3 Systemsimulation 1.3.1 Definition und Beispiele Der Begriff Simulieren kommt aus dem Lateinischen. simulare [lat.] = nachbilden, ähnlich machen Unter der Simulation eines Systems versteht man die Darstellung des Verhaltens eines physikalischen oder abstrakten Systems durch das Verhalten eines Ersatzsystems (Modell), mit dessen Hilfe man Verhaltensstudien vereinfacht durchführen will. Simulationsbeispiele der Technik mit Hilfe dynamischer Modelle sind im folgenden Diagramm zusammengestellt. Man unterscheidet physikalische und mathematische Modelle sowie gemischte Systeme, siehe Abbildung 15. Abbildung 15: Beispiele dynamischer Modelle Simulation mathematischer dynamischer Modelle „Software“ Simulation physikalischer dynamischer Modelle „Hardware“ Offline- Systemstudie Echtzeitsimulation Fahrdynamik Flugeigenschaften Körperfunktionen Ergonomie Wetter … Echtkomponententest Algorithmentest Synthetische Testsignale Prozesssimulation Fahr-, Flugsimulator Weltraumsimulator Parameter-Schätzverfahren … Windkanal Geräuschprüfstand Klimakammer Crash-Simulator … Gemischte Systeme Nachbildung von Prüflingen, Testumgebungen <?page no="31"?> 17 1. Physikalische dynamische Modelle Die Werte der Systemvariablen als Funktion der Zeit werden durch physikalische Messung bestimmt. Hierzu zählt u. a. die Simulation von Systemen bzw. Teilsystemen, Testumgebungen und Testbedingungen. Beispiele hierfür sind (zum Teil abgebildet): • Windkanal Abbildung 16 auf Seite 18 zeigt ein Schnittbild der Anlage. Starke Antriebsmotoren treiben ein Gebläserad an, das einen Luftstrom mit bis zu 250 km/ h durch einen Kanal treibt. Strömungsgleichrichter sorgen für eine konstante, gleichmäßige Geschwindigkeit über den gesamten Ausströmquerschnitt an einer Düse. Diesem Luftstrom werden Prüfobjekte ausgesetzt, z.B. ein Fahrzeug, die auf einer 6-Komponentenwaage stehen. Damit können die Strömungseigenschaften dieser Objekte vermessen werden. • Geräuschprüfstand Ein schalltoter Raum, Abbildung 17 auf Seite 19, ermöglicht Messungen von Geräuschen, die z.B. der Motor eines Fahrzeugs erzeugt, wenn er das Fahrzeug auf einer Rolle gegen einen Widerstand antreibt, das Getriebe dabei schaltet, usw. Durch eine Reihe von Mikrofonen rechts und links des Fahrzeugs kann das Geräuschniveau eines vorbeifahrenden Fahrzeugs detektiert werden. • Klimakammer In der Klimakammer können unterschiedliche Umweltbedingungen eingestellt werden, von sehr tiefen bis zu sehr hohen Umgebungstemperaturen, aber auch Regen, Schneefall kann simuliert werden, oder z.B. niedriger Luftdruck zur Bestimmung des Abgasverhaltens des Motors während einer Fahrt in sehr großer Höhe. • Crash-Simulation Fahrzeuge können von einem Schlitten sehr schnell auf gewünschte Geschwindigkeiten beschleunigt werden und auf einen festen Block treffen, siehe Abbildung 18 auf Seite 19. Dadurch kann das Crash-Verhalten des Gesamtfahrzeugs unter reproduzierbaren Bedingungen geprüft werden. Aber auch Teilsysteme, z.B. Gurtsysteme, können in Versuchsreihen optimiert werden, ohne jedes Mal ein komplettes Fahrzeug zu opfern. <?page no="32"?> 18 2. Mathematische dynamische Modelle Die Werte der Systemvariablen als Funktion der Zeit werden durch mathematische Berechnung bestimmt. Hierzu zählen u. a. die Offline-Durchführung von Simulationsstudien und die Simulation unter Echtzeitbedingungen. 1 Antriebsmotoren 4 Diffusor 7 6-Komponentenwaage 2 Antriebswelle 5 Strömungsgleichrichter 8 Vorbereitungswerkstatt 3 Gebläselaufrad 6 Düse 9 Trafostation, Versorgung Abbildung 16: Großer Windkanal in Stuttgart-Untertürkheim 1 2 3 4 5 7 6 8 9 Drehbarer Tisch Aerodynamisches Gleichgewicht Luftstrom <?page no="33"?> 19 Abbildung 17: Geräuschversuche im Akustiklabor Abbildung 18: Crash-Simulation - Verlauf abhängig von Form des Steuerstiftes x Schlitten Pneumatische Simulation des Kraftverlaufs (Vorbau-Deformation) Druckluft Rückhaltesysteme Hydraulische Crash-Simulation Fußgängerschutz x <?page no="34"?> 20 Beispiele für Offline-Simulation sind (zum Teil abgebildet): • Fahrdynamik-Simulation mit zugehöriger Animation Die Zustandsgrößen eines Fahrzeugs, hier eines Gespanns mit Anhänger, Abbildung 19 auf Seite 21, können grafisch in Form einzelner Kurven über der Zeit dargestellt werden. Aussagekräftig ist aber auch eine Animation in Form eines Films, hier am Beispiel eines Ausweichmanövers auf schneeglatter Fahrbahn. Das Fahrzeug ist gezwungen, einen Doppelspurwechsel als Manöver zu durchlaufen. In einem Film können alle interessierenden Zustandsgrößen in ein Fahrzeugmodell umgerechnet werden, das den dynamischen Ablauf des Manövers anschaulich zeigt. Ein weiteres Beispiel von Animation ist eine Bilderfolge oder ein Film eines Nutzfahrzeugs beim Überfahren einer einseitigen Rampe, Abbildung 20, Seite 21. • Kosimulation Gesamtfahrzeug Durch die am Markt ständig wachsende Anzahl verschiedener Simulationstools, die benutzerfreundlich aufgebaut und auf die unterschiedlichen Teilgebiete wie Mechanik, Hydraulik, Regelungstechnik speziell zugeschnitten sind, wird heute bei mechatronischen Systemen zunehmend jedes dieser Tools bei der Lösung der Differenzialgleichungen parallel eingesetzt. D.h., jedes Tool rechnet in seiner Welt ein Stückchen seiner Teillösung über der Zeit und gibt dieses Zwischenergebnis als Input an die anderen Tools. Dies ist der Grundgedanke der Ko- Simulation, Abbildung 21 auf Seite 22. Mechanik, Hydraulik und Regleralgorithmen werden mit drei unterschiedlichen Tools berechnet und miteinander vernetzt. Das Ergebnis der Simulation ist zugleich animiert: ein Fahrzeug fährt auf einer Ebene und zeigt sein Verhalten bei einem vorgegebenen Manöver. • Dauerfestigkeit-Simulation Als weiteres Beispiel sei ein Fahrzeug beschrieben, das in der Simulation über eine zuvor vermessene reale Straße mit Unebenheiten fährt, Abbildung 22, Seite 22. Die dabei auftretenden Kräfte werden unter dem Gesichtspunkt der Dauerfestigkeit von Karosserie oder anderen Bauteilen beurteilt. Wichtig ist, dass das Modell durch Messung der Kräfte im realen Fahrzeug verifiziert ist. Erst dann macht es Sinn, Komponentenvarianten rechnerisch untereinander zu vergleichen. Der Vorteil dieser Methode liegt darin, dass in der Simulation die immer gleiche Anregung der Straße in das Modell eingeleitet wird. <?page no="35"?> 21 Abbildung 19: Fahrdynamik-Simulation, Ergebnis-Animation Abbildung 20: NFZ-Simulation, Ergebnis-Animation 1 sec 2 sec 3 sec 4 sec 5 sec 6 sec 7 sec 8 sec T = 0.0 s 0.2 s 0.4 s 0.6 s 0.8 s 1.0 s 16 t - LKW Geschwind. v = 5 m/ s Rampenhöhe h = 0.4 m Überfahren einer Rampe <?page no="36"?> 22 Abbildung 21: Kosimulation Gesamtfahrzeugmodell Abbildung 22: Dauerfestigkeit-Simulation mit Animation Hydraulik Simulation Tool Chain Mechanik Regler Ausschnitt aus einem Film einer simulierten Fahrt über eine unebene Straße <?page no="37"?> 23 • Flugeigenschaften Die Flugeigenschaften von Flugzeugen, Raketen, Satelliten und Raumstationen werden bezüglich verschiedener Fragestellungen untersucht. Dabei sind die Modelle meist sehr aufwendig und rechenzeitintensiv. • Körperfunktionen Auch in der Biologie/ Medizin wird auf unterschiedlichsten Feldern das dynamische Verhalten simuliert, vom biologischen Gleichgewicht ganzer Populationen bis zu kleinsten Zellen und deren Stoffaustausch mit der Umwelt. • Ergonomie Abbildung 23, Seite 24: Die ergonomische Arbeitsplatzgestaltung, z.B. eines Fahrers in einem Fahrzeug, ist seit vielen Jahren wichtiger Bestandteil der Entwicklung. Erster Schritt zur Durchführung ist die Festlegung eines Modells eines Menschen mit allen geometrischen Maßen und der begrenzten Winkel bestimmter Gelenke. Dazu wird eine repräsentative Gruppe von Personen dreidimensional gescannt und die statistischen Eigenschaften dieser Gruppe erfasst. Abbildung 24, Seite 24: Dieses Modell eines Menschen wird in die Umgebung des Fahrzeugs eingebaut und alle Komponenten um es herum so lange verschoben, bis deren Anordnung optimal ist bezüglich Bewegungsraum, Kraftaufwand bei der Bedienung, Sichtbereich durch die Frontscheibe, …, und zwar sowohl für die kleine 5%-Frau bis hin zum sehr großen 95%-Mann. Abbildung 25, Seite 25: Die Arbeitsplatzgestaltung in einer Montagehalle kann vor dem Aufbau der Fertigungsstraße untersucht und optimiert werden. Parallel arbeitende Roboter müssen geschickt angeordnet werden, wenn sie zeitoptimal gemeinsam Arbeiten durchführen sollen. Außerdem wird die Zugänglichkeit geprüft und der Kraft-aufwand eines Menschen bestimmt für die Montage von Teilen am Band. • Humanoid-Dummies Die Computer werden ständig schneller und erlauben deshalb die Berechnung immer komplexerer Modelle trotz enormem Rechenaufwand, im Beispiel so genannter Humanoid-Dummies, deren Eigenschaften denen des realen Menschen immer ähnlicher werden, Abbildung 26 auf Seite 27. Problematisch dabei ist die Modellverifikation, das Modell soll sich bei einem Unfall realistisch verhalten bzgl. Reißfestigkeit des Gewebes, Bruchfestigkeit des Knochengerüsts, elastischem Verhalten der Organe untereinander, usw. Modelle dieser Art sind noch nicht Stand der Technik, weltweit gibt es nur wenig Forschungsinstitute, die auf diesem Feld arbeiten und Modelle wie MADYMO oder HUMOS entwickeln. <?page no="38"?> 24 Abbildung 23: Ergonomisches Modell eines Menschen Abbildung 24: Ergonomische Simulation Nachbildung der Bewegungsmöglichkeiten Die Körpermaße des Menschen werden heute statistisch erfasst durch dreidimensionales Scannen einer repräsentativen Gruppe. (Erste Untersuchungen bereits 1948, Harvard Public Health School, USA) d-Umwelt Ergebnisse: Kinematikdiagramm Kollisionsdiagramm Sichtdiagramm Kraftdiagramm Aufgaben: Automatische Bewegung Bewegungsraum Kollisionen Sichtbereiche Kraftsimulation Rechner 3D-Umwelt 3D-Körper Innere Kinematik Programme Ergonomische Simulation 3D-Verkopplung <?page no="39"?> 25 Abbildung 25: Ergonomie in der Fabrikmontage Aufgabe: - Zugänglichkeitsuntersuchungen im Bereich Rohbau und Montage - Ergonomie-Simulation <?page no="40"?> 26 • Insassenklima-Simulation Ein Klimamodell eines Fahrzeugs mit Fahrer erlaubt, den Klimakomfort eines zukünftigen Fahrzeugmodells im Voraus zu berechnen. So liegen z.B. dem thermophysiologischen Insassen-Modell TIM zahlreiche Erkenntnisse über das subjektive Komfortempfinden zugrunde, Abbildung 27 auf Seite 27. Es bildet den menschlichen Körper mit 14 Teilen nahezu vollständig nach und berücksichtigt überdies auch dessen Blutströme und Wärmeproduktion. Mit TIMs Hilfe erkennen Ingenieure in einem frühen Entwicklungsstadium, welche Leistung Heizung und Klimaanlage haben sollten, wie viele Belüftungsdüsen notwendig sind und wie groß sie sein müssen, um stets einen firmentypischen Klimakomfort zu erreichen. • Wetter Durch die Verbesserung der Klimamodelle konnte die Qualität der Wettervorhersage deutlich gesteigert werden. • Simulation von Vorgängen im Meer Komplexe Modelle mit vielen Einzelkomponenten ermöglichen die Berechnung von Meeresströmungen, von Wasserstandshöhen des Meeresspiegels bei Ebbe und Flut oder z.B. bei Sturmfluten der Nordsee durch Nordostwind in der norddeutschen Bucht. Ein anders Beispiel zeigt die Simulation von Eisbildung in der Ostsee, Abbildung 28 auf Seite 28: Dargestellt ist ein Ergebnis der Simulation zur Bildung von Eis sowie dessen Rückwirkung auf die Bewegungsvorgänge im Meer. Dies kann auf Basis eines sehr komplexen dreidimensionalen dynamisch-thermodynamischen Meereis-Ozean-Modells berechnet werden. Mit diesem letzten Beispiel wird die Serie von Simulationsbeispielen unter Verwendung mathematischer Modelle zur Berechnung dynamischer Abläufe abgeschlossen. <?page no="41"?> 27 Humanoid-Dummies 2 Starrkörper-Modelle (Multi-Body) aus Ellipsoiden, Gelenken, Muskeln für Parameterstudien, Optimierungen, Beispiel: MADYMO „Human Occupant“ FE-Menschenmodelle (extrem hoher Detaillierungsgrad, rechenzeitintensiv) mit Nachbildung innerer Organe, Muskeln, z.T. Blutgefäßen, Beispiel: HUMOS Human Model for Safety Problem: geschwindigkeitsabhängiges Verhalten von Gewebe aktives Muskelverhalten - Ethik: Modellverifikation, Messungen mit Freiwilligen, Leichen Abbildung 26: Humanoid-Dummies Insassenklima-Simulation Klimakomfort künftiger Fahrzeugmodelle im Voraus berechnen und optimieren Abbildung 27: Insassenklima-Simulation 2 Quelle: Dummy nicht gleich Dummy (Tecosim, Technische Simulation GmbH, Rüsselsheim), Automotive Materials, Heft 06/ 06, S. 40-42, Gisel Verlag. TIM Thermophysiologisches Insassen-Modell HUMOS FE- Menschenmodell MADYMO „Human Occupant“ <?page no="42"?> 28 Abbildung 28: Simulation von Eisbildung in der Ostsee 3 3. Gemischte Systeme Viele Simulationseinrichtungen verbinden die Echtzeit-Simulation mathematischer dynamischer Modelle (Software) mit realen Systemen und/ oder physikalischen dynamischen Modellen (Hardware). Dies erfolgt vor allem dann, wenn die mathematische Beschreibung einiger Teile des Gesamtsystems zu komplex ist bzw. Modellparameter nur sehr schwierig oder gar nicht zu bestimmen sind. Dann muss auf die Modellbildung dieser Teile verzichtet werden. Stattdessen werden echte Hardware-Komponenten in die Simulation eingebunden (HIL, Hardware in the Loop). Beispiele für gemischte Systeme sind (auf den folgenden Seiten zum Teil abgebildet): 3 Quelle: URL am 10.08.2008: http: / / www.io-warnemuende.de/ Projects/ Basys/ publications/ press/ press_de.htm Institute for Marine Research Kiel (IfM); Abt. Theoretische Ozeanographie, Duesternbrooker Weg 20, D-24105 Kiel; Germany. Eisdicke in Metern <?page no="43"?> 29 • Hydropuls-Anlage Auf einem Hydropuls-Prüfstand, Abbildung 29, kann ein reales Fahreug über vertikale Hydraulikzylinder dynamisch an den Rädern angeregt werden. Der Vorteil gegenüber einer Erprobung auf der Straße liegt darin, dass die Anregungen in immer gleicher Weise reproduzierbar in das Fahrzeug eingeleitet werden. Damit können unterschiedliche Varianten von Fahrzeugen oder auch Komponentenvarianten sehr genau objektiv miteinander verglichen und optimiert werden. Aber auch subjektive Eindrücke sind zu gewinnen, wenn man als Insasse die Anregungen der Straße auf sich wirken lässt. Abbildung 29: Hydropuls-Anlage Hydropulszylinder Aufbau: • 4 vertikale Hydropulszylinder Nennlast 40 kN, Amplitude max. ± 125 mm • Messkanäle maximal 72 • Messbereich 0,5 bis 1000 Hz Aufgabe: Straßensimulationen und Untersuchungen von Fahrzeugschwingungen und Geräuschen, sowie Betriebsschwingformen im Bereich bis 50 Hz an fahrfertigen Fahrzeugen und Komponenten • zum Vergleich verschiedener Fahrzeugvarianten (Eigen-/ Resonanzfrequenz und Schwingungsamplituden) • zur subjektiven Beurteilung des Fahrzeug-Schwingungsverhaltens • zur genauen Bestimmung der Eigenschwingungs-Kenngrößen (Eigenfrequenzen, Eigenschwingformen und Dämpfungen) <?page no="44"?> 30 • Regensimulation Zur Optimierung z.B. eines Scheibenwischers mit Regensensor benötigt man ebenfalls reproduzierbare Versuchsbedingungen. Nachdem für diese Aufgabe natürlicher Regen denkbar ungeeignet ist, wurde eine spezielle Regensimulationsanlage gebaut, wie Abbildung 30 zeigt. Hier gelingt es, für die unterschiedlichen Regenvari-anten Wischerstrategien zu entwickeln und deren Reglerparameter zu optimieren. Abbildung 30: Regensimulation • Getriebe-Simulator Die Modellierung einer hydraulischen Schaltung für ein Automatikgetriebe eines Fahrzeugs ist sehr komplex und wegen der schwer zu bestimmenden Modellparameter der Strömungsquerschnitte, Ventile, Drosseln, … nicht mit der nötigen Genauigkeit zu bestimmen. Deshalb bot sich in der Vergangenheit an, die Getriebeelektronik und die hydraulische Schaltplatte nicht zu berechnen, sondern in Hardware in den Modellverbund Fahrzeug - Getriebe - Motor zu integrieren, Abbildung 31 auf Seite 31. • Weltraumsimulator Versuche unter den Bedingungen von Schwerelosigkeit sind nur im Weltraum durchführbar und deshalb mit sehr hohen Kosten und Risiken verbunden. Für einen Versuchzeitraum mit Schwerelosigkeit von nur wenigen Minuten genügt allerdings ein Parabelflug eines Flugzeugs aus sehr großer Höhe, Abbildung 32 auf Seite 31. <?page no="45"?> 31 Abbildung 31: Getriebe-Simulator Abbildung 32: Simulation der Schwerelosigkeit Aufgabe für den Piloten ist, genau den Kurs in x- und y-Richtung zu fliegen, der in der Abbildung angegeben ist. Insassen des Flugzeugs erleben wenn auch nur für kurze Zeit exakt die Bedingungen, die Astronauten in Raumstationen dauerhaft vorfinden. BEDIENELEMENTE ANZEIGEELEMENTE TERMINAL GETRIEBE- ELEKTRONIK SCHALTPLATTE DRUCKER Digitalrechner: Simulieren des Fahrzeugs mit Motor, Getriebe, Wandler, … Prüfen der Funktion der Getriebeelektronik y y m g m ! << ( ) ( ) ( ) ( ) 2 2 0 0 0 2 t x v g 2 / 1 t y t v t x : t lg fo const v x : x y : mit t g 2 / 1 t y g m y m ≈ ≈ = ≈ = = x <?page no="46"?> 32 • Fahrsimulator, Abbildung 33 bis Abbildung 37, Seite 33 ff. Die Simulation des Zusammenspiels Fahrer - Fahrzeug - Umwelt ist besonders dann äußerst schwierig, wenn man versucht, das Verhalten des Fahrers und seine Reaktionen bei bestimmten Fahrsituationen mathematisch zu beschreiben. Die Informationsaufnahme und -verarbeitung im menschlichen Gehirn sind derart komplex, dass der Aufbau eines exakten Modells des Menschen von vorne herein zum Scheitern verurteilt ist. Es bietet sich deshalb an, dieses Verhalten unter Einbeziehung des echten Menschen in einem Fahrsimulator zu beurteilen, wie er in den folgenden Abbildungen dargestellt ist. Ziel ist, dem Menschen die gleichen subjektiven Reize zu vermitteln, wie er sie bei Fahrmanövern im realen Fahrzeug verspürt, d.h. Gleichgewichtssinn: er spürt die gleichen Beschleunigungen, Drehbeschleunigungen wie beim Bremsen oder der Kurvenfahrt Sehsinn: er sieht das (synthetische) Bild seiner Umgebung, wie es sich ihm normalerweise darstellt, auch Dunkelheit, Nebel, Gegenverkehr, … Hörsinn: er hört Geräusche, die sich situationsabhängig einstellen, z. B. quietschende Reifen beim Schleudern Tastsinn: er verspürt am Lenkrad die entsprechende Lenkradrückstellung bei Kurvenfahrt Da der Mensch in dieser Einrichtung mitwirkt, muss das gesamte Geschehen in Echtzeit ablaufen. Dies bedeutet aber auch, dass die Lösung der mathematischen Gleichungen in Echtzeit erfolgen muss! Fahrsimulator Gesamtansicht In einer Kabine ist ein vollständiges Fahrzeug eingebaut, Abbildung 33, Seite 33. Nicht nur PKW, sondern auch LKW- oder Bus-Fahrerhäuser können untersucht werden, um den Fahrern gewohnte Verhältnisse zu bieten. Der Dom kann über sechs Hydraulikzylinder in allen Freiheitsgraden bewegt werden, um die notwendigen Gier- und Längsbeschleunigungen auf den Fahrer einwirken zu lassen. Zur Simulation der Querbeschleunigung kann der komplette Aufbau auf einem Schlitten bewegt werden. Fahrsimulator Blockschaltbild Die Abbildung 34 auf Seite 33 zeigt das Zusammenspiel aller Komponenten. Herzstück bildet das mathematische Gesamtfahrzeugmodell, in dem das dynamische Verhalten in Echtzeit berechnet wird. Die Ergebnisse sind Grundlage erstens für die Berechnung der Bilder, die der Fahrer sieht (Bildrechner, Bildsystem) und zweitens für die Berechnung des Bewegungssystems (Bewegungsrechner/ Hydrauliksystem). Erwähnenswert ist auch ein Geräuschsystem, das fahrsituationsabhängig synthetische Geräusche in der Nähe der Räder nachbildet, wie z.B. Reifenquietschen im Falle dynamischer Manöver im Grenzbereich. <?page no="47"?> 33 Abbildung 33: Fahrsimulator Gesamtansicht Abbildung 34: Fahrsimulator Blockschaltbild Bildsystem Geräuschsystem Prozessperipherie Terminals Echtzeit Fahrdynamik- Rechner Bewegungsrechner Hydrauliksystem Video Monitore Steuerfaktoren Fahrzeugmodell Fahrzeugparameter Umweltmodell Umweltparameter Hydraulikaggregate Sicherheitssysteme Überwachung Sämtlicher Systeme Kommunikationssystem Fahrer<-> Versuchsleiter akustisch/ optisch Kontroll-Station Querschlitten Rechenzentrum Prüfstand Fzg-Anregung Mathematisches Fahrzeugmodell Bildrechner Steuerbefehle Kommunikationsrechner <?page no="48"?> 34 Fahrsimulator Beschleunigungssimulation Der Fahrer sieht von der Umwelt außerhalb des Domes nichts. Wird er nach vorne gekippt, Abbildung 35, so spürt er einen gewissen Anteil der Erdbeschleunigung, als ob er bremsen würde. Wird er nach hinten gedreht, glaubt er zu beschleunigen. Abbildung 35: Fahrsimulator Beschleunigungssimulation x y x y x y 0 x g y = − = α = α − = sin g x cos g y α − = α − = sin g x cos g y Konstant-Fahrt Beschleunigen Bremsen G = m g m g cos m g sin m g sin m g cos <?page no="49"?> 35 Fahrsimulator Kontrollstation, Von der Kontrollstation aus kann alles überwacht und gesteuert werden, Abbildung 36 auf Seite 36. Es gibt Monitore, die das Bild zeigen, das der Fahrer vor sich sieht. Gespräche zwischen Fahrer und Bedienpersonal sind ebenfalls möglich. Fahrsimulator Bildsystem Das Bild, Abbildung 37 auf Seite 36, setzt sich aus sechs Bildern zusammen, die für den Fahrer ein 210 Grad-Bild um das Fahrzeug herum aufbauen. Landschaft, Berge, andere Verkehrsteilnehmer, Gegenverkehr, Nebel, Fahrdynamikflächen, mehrspurige Autobahnen, etc. können ebenso dargestellt werden. Alles in allem muss ein recht großer Aufwand getrieben werden, der sich aber lohnt: - Untersuchungen können reproduzierbar unter immer gleichen Bedingungen und ohne Risiko für die Fahrer durchgeführt werden. Seine subjektiven Eindrücke und seine Reaktionen auf die unterschiedlichsten Komponenten, neue Fahrdynamiksysteme oder Assistenzsysteme, Anzeigevarianten in Head-up- Displays, etc. können nachträglich durch Befragung ermittelt werden. - Das objektive Verhalten des Fahrers kann messtechnisch erfasst werden. Die Stelldynamik der Hydraulikzylinder des Simulators muss so gewählt werden, dass fahrdynamische Manöver wie z.B. schnelle Spurwechsel möglich sind. Für Komfortuntersuchungen, die auch höherfrequente Anregung des Fahrzeugs in vertikaler Richtung erfordern, ist die installierte hydraulische Leistung nicht ausreichend. Diese Untersuchungen sind aber auch nicht Ziel von Tests im Fahrsimulator. Nach dieser Serie von gezeigten Simulationsbeispielen vom Typ gemischte Systeme, bei denen entweder Modellteile in Hardware integriert werden, um den Modellierungsaufwand und den Rechenaufwand zu begrenzen oder der Mensch direkt integriert wird, weil ein exaktes Modell von ihm nicht erzeugt werden kann, ist der Einführungsteil im Sinne eines Gesamtüberblicks über die breit gefächerten Anwendungsgebiete der Simulation abgeschlossen. Im folgenden Kapitel werden diese Anwendungsgebiete nochmals zusammengefasst und der prinzipielle Ablauf einer Systemstudie abschließend beschrieben. Danach wird speziell die Systemsimulation unter Verwendung mathematischer dynamischer Modelle sowie numerischer Lösungsverfahren vertieft. <?page no="50"?> 36 Abbildung 36: Fahrsimulator Kontrollstation Abbildung 37: Fahrsimulator Bildsystem <?page no="51"?> 37 1.3.2 Zusammenfassung Die Simulation von Systemen entspricht im Wesentlichen einer experimentellen Methode: • Man startet die Simulation unter bestimmten, zuvor festgelegten Randbedingungen und Anfangswerten. Dann beobachtet man das Systemverhalten über der Zeit. D.h., man schaut speziell auf die zeitlichen Änderungen aller Variablen des Systems. Ziel ist, das Zusammenspiel der einzelnen Komponenten der Systemstruktur zu verstehen und Beziehungen zwischen den Variablen abzuleiten. • Mit wachsendem Verständnis entstehen Ideen für neue Tests mit dem System. Man hat eine Vorstellung, wie sich das System unter den geänderten Bedingungen verhalten wird. Man experimentiert also mit dem System. • Diese weiteren Experimente bestätigen die Vorstellungen über das System oder zwingen dazu, die Vorstellungen zu modifizieren. Dabei entsteht ein neues Bild über das Zusammenspiel der Einzelkomponenten, das durch weitere Simulationsläufe bestätigt werden muss. Diese Methode wird angewandt bei: Systemanalyse Ziel ist das Verständnis der Funktion eines (existierenden oder geplanten) Systems sowie die Überprüfung von Modellen (Struktur und Daten) durch Vergleich mit dem realen System (Messungen). Systementwurf Ziel ist die Entwicklung eines Systems mit bestimmten Spezifikationen. Durch Modifikation des Modells (Struktur und Daten) wird iterativ das berechnete dem gewünschten Verhalten angepasst. Systempostulierung Ist das Verhalten eines Systems bekannt, seine Struktur aber unbekannt, wird ein auf Hypothesen basierendes Modell gebildet. Durch Vergleich des simulierten Verhaltens mit dem bekannten (gemessenen) und durch iterative Anpassung wird ein repräsentatives Modell entwickelt. Systemoptimierung Die Parameter eines gegebenen Modells werden so variiert, dass eine aus den Systemvariablen gebildete Zielgröße (Funktional) ein Extremum (Minimum / Maximum) annimmt. Optimierungsverfahren werden auch eingesetzt im Zusammenhang mit Parameterschätzverfahren (s. Systempostulierung). <?page no="52"?> 38 Die folgenden prinzipiellen Schritte zur Durchführung werden unterschieden: 1. Definition des Problems, Zielbeschreibung Je genauer die zu untersuchende Fragestellung beschrieben wird, umso größere Chancen gibt es für einen erfolgreichen Abschluss der Arbeit. 2. Planung der Studie, Arbeitsplan Es empfiehlt sich, die Vorgehensweise bereits im Vorfeld festzulegen. 3. Formulierung des mathematischen Modells Dabei muss zuerst entschieden werden, auf welche typische Weise Systemzustandsänderungen ablaufen. Je nach Systemtyp werden gemäß Tabelle 2 unterschiedliche Verfahren angewendet, die in den folgenden Kapiteln behandelt werden: System Zustandsänderung Beschreibung der Zustandsänderung Kapitel stetig kontinuierlich Differenzialgleichungen 2 diskret in Zeitabständen Logische Gleichungen, Ereignisse 3 Tabelle 2: Systemunterscheidung nach Zustandsänderung 4. Entwicklung eines Rechenprogramms 5. Validierung eines Modells (durch Vergleich mit Messungen am realen System und iterative Anpassung) 6. Entwurf der Simulationsläufe („Experimente“) unter Berücksichtigung von Rechnerkosten und Aussagesicherheit 7. Simulationsläufe und Analyse der Ergebnisse Diese Schritte finden sich auch im Blockdiagramm in Abbildung 38 auf Seite 39 wieder: Ein System kann untersucht werden, indem man es verschiedenen Tests unterzieht. Aus dem beobachteten Systemverhalten kann man dann Erkenntnisse über die inneren Vorgänge gewinnen. Es ist sinnvoll, das System parallel zu diesen Tests auch analytisch zu untersuchen. Ein mathematisches Modell erlaubt die Berechnung des Systemverhaltens mittels Simulation. <?page no="53"?> 39 Abbildung 38: Ablauf einer Systemstudie mit Hilfe der Simulation Beobachtetes Verhalten des Systems ok System Modell- Validierung ok? Ergebnisse der Simulation verwirklichen Ende Einführung einer geeigneten Simulationsmethode und Codierung Simulationsrechner Simulationsmodell Berechnetes Verhalten des Modells Mathematisches Modell Verbesserungen Mathematische Analyse Systemanalyse und -beschreibung Analytisches Modell Analytische Untersuchung Analytische Beschreibung des Systemverhaltens Zielbeschreibung für die Simulation Systembeobachtung Anpassung nötig <?page no="54"?> 40 Beschreibt das Modell das dynamische Verhalten des Systems nicht auf Anhieb genau genug (Modell-Validierung), muss das Modellverhalten verbessert werden: • Dabei kann es sein, dass es bereits genügt, nur das Simulationsmodell zu modifizieren, z.B. ein anderes Integrationsverfahren zu wählen oder eine andere Integrationsschrittweite, … • Eventuell muss das mathematische Modell aber auch um Teile erweitert werden, von denen man zuvor glaubte, sie vernachlässigen zu können. • Im schlimmsten Fall muss man sogar die Problemstellung hinterfragen und prüfen, ob das Ziel der Systemanalyse richtig formuliert ist. Erst nach erfolgter Validierung ist das Modell so weit abgesichert, dass mit ihm wietere Untersuchungen durchgeführt werden können ohne Notwendigkeit einer aufwendigen Versuchsbestätigung der Ergebnisse. Damit ist das erste Kapitel zur Einführung in die Simulation abgeschlossen. In den nächsten Kapiteln wird näher auf das Vorgehen der Simulation zeitkontinuierlicher und zeitdiskreter Systeme eingegangen. <?page no="55"?> 41 2 Simulation zeitkontinuierlicher Systeme 2.1 Einführung Definition: Ein zeitkontinuierliches (oder stetiges) System liegt vor, wenn die Aktivitäten des Systems kontinuierliche (stetige, glatte) Änderungen der Attribute des Systems bewirken. Für ein mathematisches Modell bedeutet dies: Die den Attributen zugeordneten Variablen des Modells werden durch stetige Funktionen bestimmt, welche die Änderungsraten in Form von Differenzialgleichungen (Dgln) beschreiben. Definition: Eine Differenzialgleichung (Dgl) ist eine Gleichung, die abhängige Variable (unbekannte Funktionen), unabhängige Variable und Ableitungen (=Differenziale) der abhängigen Variablen enthält. Gewöhnliche Dgln die abhängigen Variablen hängen nur von einer unabhängigen Variablen ab Partielle Dgln die abhängigen Variablen hängen von mehreren unabhängigen Variablen ab Ordnung einer Dgl höchste Ordnung der in einer Dgl auftretenden Differenziale Lösung einer Dgl explizite Darstellung der abhängigen Variablen als Funktion der unabhängigen Variablen durch Integration Anfangs- und zusätzliche Bedingungen zur Bestimmung der bei Randbedingungen der Integration entstehenden freien Konstanten (und Funktionen) Lineare Dgln enthalten alle (abhängigen und unabhängigen) Variablen in der ersten Potenz. Koeffizienten können konstant oder veränderlich sein. Einfachste Differenzialgleichungssysteme bestehen aus einer oder mehreren linearen Differenzialgleichungen mit konstanten Koeffizienten. Hier ist die Lösung oft analytisch möglich. <?page no="56"?> 42 Es folgen ein paar Beispiele für einfache zeitkontinuierliche Systeme, beschrieben durch lineare Differenzialgleichungen mit konstanten Koeffizienten: 1. Wachstum von Kapital durch Zins und Zinseszins bei stetiger Verzinsung: Gegeben: Dgl. x p x = mit x Kapital Gesucht: ) t ( f x = p Wachstumsrate Lösung: pt 0 e x x = mit 0 x Anfangskapital Zeitverlauf: gemäß Kurve A in Abbildung 39 Modelltyp häufig bei stetigen Systemen wie Bevölkerungswachstum, chemische Reaktionen, … 2. Wachstum einer Firma in einem Markt begrenzter Größe: Gegeben: Dgl. ( ) x x p x 1 − = mit x Marktanteil Gesucht: ) t ( f x = p Wachstumsrate Lösung: ( ) ( ) pt 0 1 0 e 1 x x x x − − − + = mit 0 x Anfangswert 1 x Maximalwert Zeitverlauf: gemäß Kurve B1 und B2 in Abbildung 39 0 0,25 0,5 0,75 1 1,25 0 5 10 15 20 25 Zeit t A B1 B2 C Abbildung 39: Lösungen einfacher linearer Differenzialgleichungen x 1 f(t) x 0 f(t) Zeit t <?page no="57"?> 43 3. Zerfall radioaktiven Materials durch Strahlung: Gegeben: Dgl. x p x − = mit x Masse Gesucht: ) t ( f x = p Zerfallsrate Lösung: pt 1 e x x − = mit 1 x Anfangswert Zeitverlauf: gemäß Kurve C in Abbildung 39 auf Seite 42 4. Radaufhängung, Einmassenschwinger (siehe Abbildung 4 auf Seite 5) Gegeben: Dgl. ( ) t F x K x D x M = + + umgeformt: ( ) t F M 1 x M K x M D x = + + bzw. ( ) t F K x x 2 x 2 2 ω = ω + ω δ + mit Dämpfungsmaß δ und Kreisfrequenz ω Gesucht: ) t ( f x = für F(t) = 0 für t = 0 und F(t) = F 0 für t > 0 Fall: plötzlich aufgebrachte Kraft (Beladung) Lösung: die „Sprungantwort“ gemäß Abbildung 40 läuft für 1 < δ als gedämpfte Schwingung mit der Frequenz π ω = 2 f 1 = δ als aperiodischer Grenzfall und 1 > δ als aperiodischer Fall (nicht schwingend) Abbildung 40: Radaufhängung x 1 , 0 = δ 25 , 0 = δ 5 , 0 = δ 0 , 1 = δ 0 , 2 = δ F 0 / K t Sprungantwort des Einmassenschwingers = aperiodischer Grenzfall <?page no="58"?> 44 Komplexere zeitkontinuierliche Systeme werden beschrieben durch • Lineare Dgln mit variablen Koeffizienten z.B. x ) t ( p x = • Nichtlineare Dgln z.B. 2 1 2 1 1 1 y y c y c y + = oder x sin c x = Hilfsmittel der Simulation zur Lösung von Dgln und Dgl-Systemen sind der Digitalrechner vor allem unter Verwendung höherer Programmiersprachen (FORTRAN, PASCAL, ACSL, MatlabSimulink, …) Analogrechner (bis in die 70er Jahre) für nicht zu große, komplexe Systeme mit wenig Nichtlinearitäten Vorteil: Ablauf immer in Echtzeit möglich 2.2 Simulation auf dem Digitalrechner 2.2.1 Numerische Integration Unter numerischer Integration versteht man die näherungsweise Berechnung eines Integrals durch schrittweises Aufsummieren von Teilflächen über der Zeitachse. Gesucht ist die Fläche unter der Kurve in Abbildung 41 ( ) dt t y F Tb Ta = . Abbildung 41: Bestimmtes Integral y T a T b t y(t) <?page no="59"?> 45 Näherungslösung: Das Intervall von T a (=T 0 ) bis T b (=T n ) wird in n gleiche Teile der Breite h zerlegt, siehe Abbildung 42. Für die Stützstellen T 0 , T 1 , … , T n werden die Funktionswerte y 0 , y 1 , … , y n berechnet. Mit Kenntnis dieser Funktionswerte lassen sich Teilflächen einfacher geometrischer Form berechnen. Der numerische Näherungswert für das Integral ergibt sich dann aus der Summe über alle Teilflächen zu ( ) = ≈ = n 1 i i Tb Ta F dt t y F . Die Berechnung der Fläche (=Integral) erfolgt bei Kenntnis der Funktionswerte an den Stützstellen z.B. nach folgenden einfachen Regeln: 1. Rechteckformel, Abbildung 42 ( ) ≈ = dt t y F Tb Ta h ( y 0 + y 1 + … + y n-1 ) F i = h y i-1 Abbildung 42: Integration mit Rechteckformel 2. Trapezformel, Abbildung 43, Seite 46 ( ) ≈ = dt t y F Tb Ta h/ 2 ( y 0 + 2 y 1 + … + 2 y n-1 + y n ) F i = h [ y i-1 + 1/ 2 (y i - y i-1 ) ] y t T 0 T a T 1 T b T n-1 T n <?page no="60"?> 46 Abbildung 43: Integration mit Trapezformel 3. Parabelformel, Simpson-Regel (n gerade! ), Abbildung 44 ( ) ≈ = dt t y F Tb Ta h/ 3 ( y 0 + 4 y 1 + 2 y 2 + 4 y 3 + … + 2 y n-2 + 4 y n-1 + y n ) Abbildung 44: Integration mit Simpson-Regel y t T 0 T a T 1 T b T n-1 T n y t T a T b i-2 i-1 i Parabel 2. Ordnung <?page no="61"?> 47 Die Genauigkeit der Integralrechnung nimmt mit der Anzahl der Teilflächen, d.h. mit der Anzahl n, zu. Die Genauigkeit wird allerdings mit wachsendem n wieder durch die zunehmenden Rundungsfehler der großen Anzahl von Rechenoperationen begrenzt. Da für diesen Fall die Rechenzeit größer wird, gilt es, die Anzahl n der Teilflächen so zu wählen, dass ein guter Kompromiss zwischen Genauigkeit und Rechenzeit erreicht wird. In Programmbibliotheken und Simulationssprachen ist eine Vielzahl zum Teil komplexer und aufwendiger Integrationsverfahren verfügbar, z.B. • Einschrittverfahren, Mehrschrittverfahren • Verfahren ohne / mit Verwendung von Ableitungen (Ableitungen sind meist nicht explizit angebbar! ) • Verfahren mit fester oder variabler Schrittweite (selbsttätige Schrittweitenanpassung innerhalb unterer / oberer Schranken) • Spezielle Verfahren für „steife“ Systeme (hohe Frequenzen) Die Auswahl geeigneter Integrationsverfahren hängt ab von • den Eigenschaften des zu simulierenden Systems und • den Anforderungen an die Genauigkeit und kann zu gravierenden Unterschieden bezüglich des Rechenzeitaufwandes führen. Vergleichsrechnungen mit unterschiedlichen Verfahren und verschiedenen Genauigkeitsschranken, aber auch Erfahrung sind bei der Auswahl hilfreich. 2.2.2 Simulationssprachen für kontinuierliche Systeme Hierunter werden problemorientierte Programmsysteme zur Untersuchung des dynamischen Verhaltens zeitkontinuierlicher Systeme verstanden. Im Unterschied zu universellen Programmiersprachen enthalten diese Simulationssprachen Sprachelemente, die eine direkte Formulierung von Simulationen gestatten. Eigenschaften: Anwendernähe Reduzierung des Programmieraufwandes Einfache Erlernbarkeit Flexibilität Komfortable Ein-/ Ausgaberoutinen Unterstützung bei Modellerstellung Erste Ansätze zur Entwicklung derartiger Sprachen gibt es seit 1955. Damals lehnten sich die Sprachen noch stark an die Programmierung von Analogrechnern an und wurden daher zunächst als „digitale Analog-Simulatoren“ aufgefasst. Tabelle 3 auf Seite 48 gibt einen (groben) Überblick über wichtige Entwicklungen der letzten 50 Jahre auf dem Gebiet der Simulationssprachen kontinuierlicher Systeme: <?page no="62"?> 48 1955 Selfridge Bildet Blockdiagramme auf Digitalrechnern nach 1959 Stein, Rose, Parker Compiler mit Eingabesprache, die sich an den Analogrechner anlehnt 1964 COBLOC Simulationssprache auf Basis CDC 1604 mit 3 Integrationsverfahren und grafischer Ausgabe 1967 CSMP Continuous System Modelling Program, IBM 1967 CSSL Normierungsversuch für Simulationssprachen 1971 CSMP III Weiterentwicklung von CSMP 1975 ACSL Advanced Continuous Simulation Language Mitchell and Gauthier Ass., MA, Basis: CSSL 1980 DSL Digital Simulation Language, IBM-Weiterentwicklung Von CSMP 1980 ISIS Versuch neuer Standards, Roy Crosbie 1985 SIMTRAN Simulationssprache für SIMSTAR, kombiniert ACSL mit ECSSL 1990 MATLAB SIMULINK Blockorientiert mit interaktiver Modellbildung Seit 1990 Zunehmend Visualisierung, Animation, parallele, verteilte Simulation, Kosimulation für die Integration unterschiedlicher Simulationssysteme Tabelle 3: Entwicklung der Software-Tools zur Simulation Die Vielzahl der Entwicklungen führte zu Standardisierungsempfehlungen, die als CSSL-Standards von 1967 und 1980 bekannt sind (CSSL = Continuous System Simulation Language). Danach besteht ein CSSL-Programm typischerweise aus drei Sektionen, der INITIAL-, der DYNAMIC- und der TERMINAL-Sektion. Diese werden nacheinander ausgeführt, siehe Abbildung 45 auf Seite 49, und haben folgende Funktionen: INITIAL Vor Beginn der Simulation: Setzen der Modellparameter, Anfangs- Parameter, Simulationsparameter (z.B. Simulationszeit, …) DYNAMIC Programm-Code, der während des Simulationslaufs ausgeführt wird TERMINAL Nach Beendigung der Simulation: Auswerten, Aufbereiten und eventuelle Weiterverarbeitung der Ergebnisse Elemente (Auszug): INTEG Integrator LINTEG begrenzter Integrator DER Ableitung IMPL Lösung impliziter Gleichungen DELAY Verzögerung HYST Hysterese RST Flip-Flop DEAD Totzeit RND Zufallszahlen-Generator Oft werden weit über 100 Elemente angeboten. <?page no="63"?> 49 Integrationsverfahren: Üblich: 5 - 12 Verfahren (z.B. Euler, Runge-Kutta, Adams, Gear, …) Empfohlen: Vergleichende Anwendung von mindestens 3 Verfahren Weitere Integrationsverfahren für Spezialanwendungen sind vom Benutzer in Form von Subroutinen, … hinzufügbar. Abbildung 45: Ablauf einer Simulation nach CSSL-Standard In den folgenden beiden Kapiteln werden exemplarisch die Programmiersprachen ACSL und Matlab Simulink betrachtet. 2.3 ACSL Advanced Continuous Simulation Language 2.3.1 Einführung „Advanced Continuous Simulation Language“, gemäß CSSL-Standard von 1967, auf FORTRAN-Basis Man unterscheidet: • Modell-Definition das eigentliche ACSL-Programm • Run-Time-Befehle interaktive Kommandos zur Laufzeit In der Modell-Definition werden die dem zu simulierenden System entsprechenden Gesetzmäßigkeiten formuliert. In der Regel sind dies die Differenzialgleichungen des Systems. Die Struktur richtet sich bei ACSL nach dem CSSL-Standard. DYNAMIC Statements TERMINAL Statements INITIAL Statements <?page no="64"?> 50 Die wichtigsten Modell-Statements sind: PROGRAM Beginn eines ACSL-Programms INITIAL Beginn der INITIAL-Sektion eines ACSL-Programms CONSTANT Definition und Vorbelegung von Variablen DYNAMIC Beginn der DYNAMIC-Sektion eines ACSL-Programms TERMT Formulierung von Endbedingungen für die Simulation DERIVATIVE Beginn der DERIVATIVE-Sektion eines ACSL-Programms INTEG Integrationsoperator TERMINAL Beginn der TERMINAL-Sektion eines ACSL-Programms END Beendet jeweils die mit PROGRAM, INITIAL, DYNAMIC, DERIVATIVE oder TERMINAL begonnenen Sektionen Mit den Runtime-Befehlen ist ein weitgehender Zugriff auf das Programm zur Laufzeit möglich. So können z.B. einzelne Modellparameter angesehen oder modifiziert werden und die Simulation mit den geänderten Werten erneut gestartet werden. Die wichtigsten Runtime-Kommandos sind: START Start eines Simulationslaufs STOP Beendigung eines Simulationslaufs PREPAR Abspeichern von Daten in einem File während der Simulation zur späteren Weiterverarbeitung PRINT Ausdrucken der Ergebnisse nach Ende des Simulationslaufs PLOT Plotten der Ergebnisse nach Ende des Simulationslaufs DISPLAY Anzeigen der Modell- oder Systemparameter SET Ändern der Modell- oder Systemparameter 2.3.1.1 ACSL Modelldefinition Struktur eines ACSL-Modells PROGRAM . Pre-INITIAL-SECTION . Anweisungen werden nur einmal je Aufruf des . Programms ausgeführt . INITIAL . INITIAL SECTION . Anweisungen werden vor jedem Simulationslauf . (START) ausgeführt END . . <?page no="65"?> 51 . DYNAMIC . DYNAMIC SECTION . Enthält . den dynamischen kontinuierlichen Modellteil . alle Anweisungen, die zur Berechnung der Zustandsvariablen benötigt werden . - Anweisungen, die zu jedem Kommunikationszeitpunkt ausgeführt werden (typisch: Umrechnungen für die Ausgabe) END TERMINAL . TERMINAL SECTION . Enthält Anweisungen, die nach der Beendigung der Simulation . ausgeführt werden . END END In der Pre-INITIAL SECTION werden Anweisungen nur einmal je Simulationslauf ausgeführt. In dieser Sektion werden also z.B. Werte von externen Dateien gelesen (große Tabellen) oder Berechnungen durchgeführt, die nur einmal zu Beginn einer Simulation notwendig sind. Die INITIAL SECTION wird jedes Mal nach Eingabe des Runtime-Kommandos START durchlaufen. In die INITIAL SECTION werden alle die Berechnungen aufgenommen, die von zur Runtime änderbaren konstanten Werten abhängen. In der DYNAMIC SECTION stehen alle Modellteile, die zur Berechnung der Zustandsvariablen benötigt werden, und die Integratoren. In jedem Durchlauf werden die speziell ausgewählten Integrationsverfahren aufgerufen. Den Abstand der Zeitpunkte, zu denen die DYNAMIC SECTION abgearbeitet wird, definiert das so genannte Kommunikationsintervall. Zu jedem Kommunikationszeitpunkt werden alle Anweisungen der DYNAMIC SECTION ausgeführt und die gewünschten Ergebnisse auf die Ausgabe- und die PREPAR-Datei geschrieben. Das Kommunikationsintervall wird durch den ACSL-Systemparameter CINT festgelegt. Die Integrationsschrittweite für Algorithmen mit fester Schrittweite ist dann CINT/ NSTP, wobei NSTP die Anzahl der Integrationsschritte innerhalb des Intervalls CINT ist und mit einem Wert von 10 vorbesetzt ist. NSTP kann allerdings vom Benutzer geändert werden. Die TERMINAL SECTION wird nach Beendigung der Integration abgearbeitet. Hier können z.B. Berechnungen aufgenommen werden, die nur einmal am Ende der Integration ausgeführt werden sollen. <?page no="66"?> 52 ACSL-Statements zur Modelldefinition Die auf Seite 135 im Anhang 4.1 befindliche Tabelle 17 gibt beispielhaft Statements und Funktionen wieder, die vom ACSL-Übersetzer akzeptiert werden. Einige davon sind identisch mit den entsprechenden FORTRAN-Statements bzw. -Functions, mit der Ausnahme, dass der ACSL-Übersetzer ein freies Format im Programmcode akzeptiert. Die Sortierung der Statements entsprechend der Kausalität der Berechnung wird dem Benutzer vom ACSL-Übersetzer abgenommen. Alle Statements werden innerhalb der DYNAMIC SECTION so sortiert, dass alle benötigten Werte zum Zeitpunkt der Berechnung bekannt sind bzw. vorher bereits berechnet wurden. Wo dies nicht möglich ist (z.B. arithmetische Schleifen), gibt der Übersetzer eine Fehlermeldung aus. Die anderen Sektionen des Programms werden - wie sonst üblich - sequenziell durchlaufen. Schritte zur Erzeugung eines lauffähigen Programms Aus dem ACSL-Programmcode wird in einem 3-stufigen Prozess ein lauffähiges ACSL-Programm erzeugt, siehe Abbildung 46 auf Seite 53. Zuerst generiert der ACSL-Übersetzer aus dem ACSL-Programm ein FORTRAN- Programm. Hierbei wird vor dem Übersetzungsvorgang das ACSL-Programm um die Makro-Aufrufe textuell ergänzt. Im zweiten Schritt übersetzt der FORTRAN-Compiler das FORTRAN-Programm und erzeugt ein maschinenspezifisches Objektprogramm. Schließlich wird im dritten Schritt durch den Linker aus dem Objektprogramm und den benötigten Routinen aus der ACSL-Library, der System-Library und eventuellen zusätzlichen User-Libraries ein ausführbares Programm erzeugt. In allen drei Schritten können Fehlermeldungen auftreten, und es ist keineswegs garantiert, dass ein vom ACSL-Übersetzer als fehlerfrei akzeptiertes Programm auch im nachfolgenden Compile- oder Linkvorgang fehlerfrei ist. 2.3.1.2 ACSL Runtime-Kommandos Zur Laufzeit (Runtime) steht ein durch den ACSL-Übersetzer und den FORTRAN- Compiler übersetztes und mit ACSLsowie FORTRAN-Bibliotheken gelinktes, ausführbares ACSL-Programm zur Verfügung. Der Benutzer steuert die Ausführung dieses Programms und die bisher noch offenen Details der Simulation durch einen Satz so genannter Runtime-Kommandos. Diese werden vom Kommando-Interpreter des ausführbaren ACSL-Programms zeilenweise zur Laufzeit gelesen und sofort ausgeführt; Verzweigungen und Schleifen sind nicht möglich. <?page no="67"?> 53 Abbildung 46: Ablaufdiagramm zur Erzeugung eines lauffähigen ACSL-Programms Durch Runtime-Kommandos können vor allem Parameterwerte des Modells geändert werden. Diese behalten ihren Wert bei, bis ein nachfolgendes Kommando ihn erneut ändert oder ein Programmstück in der Modelldefinition diesen Parameter neu berechnet. 1. Schritt 2. Schritt 3. Schritt ACSL-Modell ACSL- Übersetzer Fehler- Meldungen FORTRAN- Quelle Fehler- Meldungen FORTRAN- Übersetzer ACSL- Makros Objekt- Programm Fehler- Meldungen Linker System User- Libraries ACSL Runtime- Library Ergebnisse Simulations- Programm Runtime- Kommandos Plots und Ergebnislisten <?page no="68"?> 54 Runtime-Kommandos werden solange gelesen und ausgeführt, bis durch die Eingabe von START die Kontrolle an die Modellsektion abgegeben wird und der Ablauf nach Abbildung 47 beginnt. Die Eingabe von STOP beendet das Lesen von Kommandos und auch die Programmausführung. Abbildung 47: ACSL Simulationsablauf START INITIAL Section CONTINUE Pre-INITIAL Section Set IC DYNAMIC Section Eintrag in PREPAR- und OUTPUT-Liste Integration t= t+h ja nein TERMINAL Section END Ende? STOP? <?page no="69"?> 55 Runtime-Kommandos bestehen aus einem Kommando, z.B. PRINT, und einer mit Leerzeichen abgesetzten und untereinander durch Komma getrennten Liste von Variablennamen, in der auch Schlüsselwörter zur Detaillierung der Kommandofunktion vorkommen können. Die allgemeine Syntax lautet: Command Var1, Var2, … , ’Key1’=Value1, ’Key2’=Value2, … , Varn . Eine typische Kommando-Sequenz zur Ausführung eines ACSL-Programms mit Protokollierung der Ergebnisse in Form einfacher Liste und Plots ist wie folgt: • PREPAR ’CLEAR’, ’ALL’, <Liste> Beispiel: PREPAR T, X, XD, XDD Erklärung: Alle Variablen, die später gedruckt oder geplottet werden sollen, müssen in die PREPAR-Liste. Diese Variablen werden während der Simulation zu jedem Kommunikationszeitpunkt auf ein internes File geschrieben. • OUTPUT ’CLEAR’, ’ALL’, ’NCIOUT’=i, <Liste> Beispiel: OUTPUT T, X, ’NCIOUT’=10 Erklärung: Die Variablen, die in der OUTPUT-Anweisung stehen, erscheinen sofort während der Simulation auf dem Bildschirm. Sie werden nicht wie die Variablen der PREPAR-Liste für eine spätere Weiter verarbeitung (PRINT, PLOT) aufbewahrt. • START Dieser Befehl beendet das Lesen von Runtime-Kommandos und startet die Simulation. • PRINT ’ALL’, ’NCIPRN’=i, <Liste> Beispiel: PRINT T, X, XD, ’NCIPRN’=2 Erklärung: PRINT liest vom internen File, das während des letzten Simulationslaufs beschrieben wurde, die Werte der angegebenen Variablen und schreibt diese formatiert in die Ausgabedatei (.PRN). Die Varia blen müssen also vor dem letzten START bereits in der PREPAR-Liste gestanden haben. • PLOT ’XAXIS’=var, ’ALL’, <Liste> Beispiel: PLOT X, XD Erklärung: PLOT liest vom internen File, das während des letzten Simulationslaufs beschrieben wurde, die Werte der angegebenen Variablen und erzeugt einen Plot in der gewünschten Form, die durch Systemvariablen genauer festgelegt wird. Plot X, XD plottet X und XD über der ersten Variablen der PREPAR- Liste (in der Regel die Zeit T). Eine Änderung ist durch Angabe von ’XAXIS’ möglich. Die Variablen müssen vor dem letzten START bereits in der PREPAR-Liste gestanden haben. <?page no="70"?> 56 • DISPLY ’ALL’, <Liste> Beispiel: DISPLY A, B, C, NCIOUT, CINT, D A(5) Erklärung: DISPLY zeigt den momentanen Wert der angegebenen Variablen am Bildschirm • SET var1=val1, … , varn=valn Beispiel: SET A=5.0, CINT=3.0, PRN=9, P=3*0.5*A S CALPLT=.T. Erklärung: Wertzuweisung für Variable • STOP Beendet den Simulationslauf Weitere Beispiele zeigt die Tabelle 18 im Anhang 4.2 auf Seite 136. 2.3.1.3 Simulationsablauf Nach Eingabe des Kommandos START am Terminal, siehe Abbildung 47 auf Seite 54, durchläuft das Programm zunächst die INITIAL-Section. Danach werden die Anfangswerte der Zustandsvariablen berechnet und in deren Speicherplätze geschrieben. In der folgenden DYNAMIC-Section wird zunächst geprüft, ob das STOP-Bit gesetzt ist. Ist dies nicht der Fall, wird die DYNAMIC-Section abgearbeitet, wozu das Schreiben der in der PREPAR-Liste und OUTPUT-Liste vorhandenen Variablen gehört. Das System führt dann unter Steuerung der Integrationsroutine die vorhandenen Anweisungen aus, d.h. es wird integriert. Dies geschieht bis zum nächsten Kommunikationszeitpunkt, an dem wieder die Ende-Bedingung geprüft wird und die OUTPUT- und PREPAR-Liste erneut beschrieben werden. Bei Erreichen einer Stopp-Bedingung wird die TERMINAL-Section durchlaufen, in der abschließende Berechnungen unter Verwendung bisheriger Ergebnisse möglich sind. Danach wird die Kontrolle wieder an das Terminal übergeben. 2.3.1.4 ACSL Systemvariable Die Systemvariablen (Beispiele in Tabelle 4 auf Seite 57 ) dienen zur - Integrationssteuerung (Algorithmus, Schrittweite), - Ausgabesteuerung (Drucken, Plotten), - Integrationskontrolle (Genauigkeit, …). Sie haben einen vorbesetzten Standard-Wert und können sowohl in der Modelldefinition als auch zur Laufzeit geändert werden. <?page no="71"?> 57 Name Bedeutung Default- Wert Ändern im Programm Ändern zur Runtime CINT Kommunikationsintervall, Integrationsschrittweite H = CINT/ NSTP 0.1 CINTERVAL name=.. SET name=.. IALG Nummer des Integrations-Algorithmus 5 ALGORITHM name=.. SET name=.. MAXT Maximale Schrittweite 1.e10 MAXTERVAL name=.. SET name=.. MERROR Relativer Fehler, der bei IALG=1 oder IALG=2 erreicht werden soll 1.e-4 MERROR X=..,XD=.. MERROR X=.. MINT Minimale Schrittweite 1.e-10 MINTERVAL name=.. SET name=.. NSTP Anzahl der Integrationsschritte je Kommunikationsintervall 10 NSTEPS name=.. SET name=.. T Unabhängige Variable 0.0 VARIABLE name=.. SET name=.. XERROR Absoluter Fehler, der bei IALG=1 oder IALG=2 erreicht werden soll 1.e-4 XERROR X=.., XD=.. XERROR XD=.. Tabelle 4: ACSL Systemvariable (Beispiele) Die Ausgabeschrittweiten und die Integrationsschrittweite werden gesteuert in Abhängigkeit der zentralen Größe CINT, dem Kommunikationsintervall. Dies zeigt Tabelle 5: Berechnung T in s Schritte pro s Kommunikationsintervall CINT 0,025 40 Integrationsschrittweite CINT / NSTP 0,025 / 10 = 2,5 ms 400 Ausgabeschrittweite OUTPUT (Bildschirm) CINT * NCIOUT 0,025 * 40 = 1 1 Ausgabeschrittweite PRINT (Drucker) CINT * NCIPRN 0,025 * 20 = 0,5 2 Tabelle 5: ACSL Schrittweitensteuerung <?page no="72"?> 58 2.3.2 Beispiel Gedämpftes nichtlineares Pendel Dies ist ein Beispiel, das für Schulungszwecke gut geeignet ist. Es geht um ein gedämpftes nichtlineares Pendel gemäß Abbildung 48. Abbildung 48: Gedämpftes nichtlineares Pendel Die zugehörige Differenzialgleichung lautet: ( ) Θ + Θ − = Θ l K sin G m l 1 Die Lösung der Differenzialgleichung beschreibt das Verhalten des Pendels nicht nur für kleine Winkel hier würde eine lineare Dgl ausreichen -, sondern auch für beliebig große Winkel. Möglich ist somit auch die Überkopflage oder das umlaufende Pendel. Das zugehörige ACSL-Programm unter Verwendung anderer Namen der Variablen und anderer Einheiten lautet: Program Damped Non-linear Pendulum Initial Constant Weight = 5.0, G = 32.2 Mass = Weight / G End Dynamic Cinterval Cint = 0.025 Constant Tstp = 4.99 Termt ( T .ge. Tstp) l K Θ Gedämpftes nichtlineares Pendel m Masse des Pendels l Länge der Pendelstange K Dämpfungskonstante G Gewichtskraft m Dämpfer G = m g <?page no="73"?> 59 Algorithm Ialg = 4 Maxterval Maxt = 0.0125 Nsteps Nstp = 1 Constant Length = 2.0, Kdamp = 0.03 THDD = - (Weight * SIN(TH) + Kdamp * Length * THD) / (Length * Mass) Constant THDIC = 0.0, THIC = 1. THD = INTEG (THDD, THDIC) TH = INTEG (THD, THIC) End Terminal Logical Dump Constant Dump = .False. If (Dump) Call Debug End End Das Verhalten des Pendels für die spezielle Anfangsbedingung das Pendel hängt nach unten THIC = 0. [rad] und es erhält einen Stoß, so dass THDIC = 10. [rad/ s] ist in Abbildung 49 dargestellt. Abbildung 49: Verhalten des nichtlinearen Pendels TH THD THDD <?page no="74"?> 60 Wichtig ist, dass das Ergebnis von Abbildung 49 auf Seite 59 genau betrachtet wird, physikalisch interpretiert wird und am Ende klar ist, warum der Kurvenverlauf genau so und nicht anders aussieht. Ich möchte den Leser sensibilisieren darauf, den Kurvenverlauf jeder Größe des Pendels genau zu betrachten und zu prüfen, ob er dem entspricht, was man bei den gewählten Anfangsbedingungen erwartet. Deckt sich der Verlauf mit diesen Vorstellungen, darf man am Ende sicher sein, dass die Differenzialgleichung und das Programm zur Lösung richtig sind und man auch anderen Ergebnissen vertrauen darf. Deshalb wird nun sehr ausführlich beschrieben, was genau im Laufe der Zeit mit dem Pendel passiert. Man erkennt, dass das Pendel durch die große Anfangsgeschwindigkeit zweimal umläuft: • Die Kurve für den Winkel TH geht zum Zeitpunkt T = 0 durch 0, dem Anfangswert THIC für den Winkel TH. Der Winkel nimmt zu, d.h., das Pendel läuft wegen der großen Anfangsgeschwindigkeit THDIC = 10 rad/ s in Überkopflage und darüber hinaus, es läuft also um. Wegen der wirkenden Dämpfung verliert es an Energie und wird langsamer. Ungefähr ab dem Zeitpunkt T = 3 Sekunden ist es so langsam, dass das Pendel es nicht nochmals bis in Überkopflage schafft, sondern zu schwingen beginnt um einen Wert TH, der einem ganzzahligen Vielfachen von 2 entspricht, hier 4 . Es ist also zweimal umgelaufen, bevor es mit zunehmender Zeit beim Winkel TH = 4 zur Ruhe kommt. • Dies bestätigt sich bei der Betrachtung der Kurve für die Winkelgeschwindigkeit THD. Der Anfangswert THDIC ist 10 rad/ s. Die Geschwindigkeit nimmt auf dem Weg zur Überkopflage ab und nimmt wieder zu, wenn das Pendel auf der anderen Seite wieder nach unten läuft. Die Geschwindigkeit erreicht aber nicht mehr den Wert THD = 10 rad/ s, wenn das Pendel nach einem vollständigen Umlauf wieder den unteren Punkt durchläuft. Es hat durch die Dämpfung an Energie und somit an Geschwindigkeit verloren. Die Geschwindigkeit THD ändert in dieser Phase ihr Vorzeichen nicht, ist positiv, d.h., das Pendel hat eine konstante Umlaufrichtung im Gegenuhrzeigersinn. Dies ändert sich auch nicht im Verlauf des zweiten Umlaufs. Erst danach ist die Geschwindigkeit so klein geworden, dass ein dritter Umlauf nicht mehr möglich ist. Kurz nach T = 3 s wird die Geschwindigkeit Null, d.h., das Pendel hält an, um von da an im Uhrzeigersinn zurück zu laufen, die Geschwindigkeit wird negativ und ändert von da an ständig das Vorzeichen, d.h., das Pendel schwingt hin und her. • Am schwierigsten ist vielleicht die Interpretation der Kurve für die Winkelbeschleunigung THDD. Es wird leichter verständlich, wenn man zunächst überlegt, wann die Beschleunigung maximal ist oder Null. Die Gewichtskraft G ist verantwortlich dafür, dass sich das Pendel bewegt. Das Pendel wird maximal beschleunigt, wenn die Gewichtskraft tangential in Richtung der Bahnkurve zeigt. Dies ist der Fall, wenn das Pendel rechts oder links bis zur Horizontalen ausgelenkt ist (3Uhrbzw. 9 Uhr-Stellung). Immer, wenn die Kraft normal, also senkrecht zur Bahnkurve zeigt, ist die Beschleunigung Null (6 Uhrbzw. 12 Uhr-Stellung). Dies ist auch der Grund dafür, <?page no="75"?> 61 dass das Pendel ohne äußere Krafteinwirkung und ohne Anfangsgeschwindigkeit stabil in 6 Uhr-Stellung verharrt und nicht von alleine losläuft. Zwischen diesen Extremwerten des Winkels verläuft die Beschleunigung sinusförmig, wie auch die Differenzialgleichung aussagt. Der Anfangswert von THDD, rote Kurve, ist Null, denn das Pendel startet unten. Auf dem Weg zur 3 Uhr- Stellung wird es zunehmend stark verzögert und erreicht dann einen ersten Maximalwert der Verzögerung. Grob gesprochen verläuft die Beschleunigung THDD während des mehrfachen Umlaufs sinus-förmig. Dies charakterisiert die erste Phase „Umlaufendes Pendel“. Wenn das Pendel beginnt, hin und her zu schwingen, erreicht es zunächst (im Bereich T = 4 bis 6 s) auch Stellungen oberhalb der Horizontalen (9 Uhrbzw. 3 Uhr-Stellung). In dieser Phase „Pendel schwingt über die Horizontale hinaus“ zeigen sich im Verlauf der roten Kurve „Dellen“, deren Größe abnimmt mit der Zeit, weil das Pendel wegen der kleiner werdenden Geschwindigkeiten am Tiefpunkt (6 Uhr-Stellung) nicht mehr die Horizontale erreichen kann. In der dritten, abschließenden Phase „Pendel schwingt unterhalb der Horizontalen“ wird der Verlauf der Winkelbeschleunigung einer gedämpften Sinusschwingung immer ähnlicher. Das Pendel kommt nach endlicher Zeit bei einem Winkel von 4 vollständig zur Ruhe. Am Beispiel dieses Pendels kann man relativ leicht nachvollziehen, was genau physikalisch passiert. Je komplexer ein dynamisches mathematisches Modell ist bei einer zunehmenden Anzahl sich gegenseitig beeinflussender Zustandsgrößen, umso schwieriger gestaltet sich die Interpretation der einzelnen Effekte. Aber auch in diesen Fällen muss der Anwender sich dieser Herausforderung stellen, wenn er sicher gehen will, dass sein Programm richtig ist. Mit Runtime-Kommandos kann man die Parameter des Pendels und die Anfangsbedingungen modifizieren. Diese Kommandos müssen über die Tastatur eingegeben werden. Sie sind zum Teil sehr kurz und einfach, sie können aber auch aufwendig lang werden. Und Korrekturen bei Tippfehlern sind nicht möglich. Um den Umgang mit dem Programm einfacher zu gestalten, gibt es die Möglichkeit, zur Runtime Prozeduren zu vereinbaren. Komplizierte Kommandosequenzen können so vorgefertigt werden durch einmalige Eingabe über Tastatur. Der wiederholte Aufruf erfolgt sehr einfach durch Eingabe des Prozedurnamens. Diese Prozeduren gehen allerdings nach Beendigung der Runtime wieder verloren. Befehlsfolgen und komplette Prozeduren lassen sich auch vorfertigen und in speziellen Files abspeichern. Zur Ausführung aller Kommandos des Files zur Runtime muss der Systemvariablen CMD ein Wert 0 zugewiesen werden. Mit CMD = 0 (Default-Wert) erfolgt die Eingabe über Tastatur. Das Runtime-Kommando SET CMD = 5 (z.B.) bewirkt, dass das Programm nach dem Namen des Files fragt, welche der Nummer 5 zugeordnet werden soll. Gibt man nun den Namen an, so führt das Programm alle Kommandos aus, die in diesem File gelistet sind. In der letzten Zeile muss die Befehlsgewalt wieder an die Tastatur zurückgegeben werden über das Kommando SET CMD = 0. <?page no="76"?> 62 Beispiel eines Kommando-Files: PREPAR T, TH, THD, THDD OUTPUT T, TH, THD, 'NCIOUT'=20 S HVDPRN=.T. $ 'PRINT auch auf Bildschirm' PROCED PL PLOT 'XAXIS'=TH, 'XTAG'='(RAD)', THD, 'TAG'='(RAD/ SEC)' $ PAUSE END PROCED STRIP SET STRPLT=.T., CALPLT=.F. PLOT 'XAXIS'=T, 'XTAG'='(SEC)', TH, 'TYPE'=6, THD, 'TYPE'=2 $ PAUSE PLOT TH, 'TYPE'=6, THDD, 'TYPE'=2 $ PAUSE END PROCED PHASE SET STRPLT=.F., CALPLT=.T. PLOT 'XAXIS'=TH, 'XTAG'='(RAD)’, THD,'TYPE'=6, 'TAG'='(RAD/ SEC)' $ PAUSE SET ’XAXIS’=T END SET CMD=0 Nach Aufruf des Command-Files und nach Ausführung aller Kommandos ist nicht nur das PREPAR-File und das OUTPUT-File festgelegt, sondern auch drei Prozeduren PL, STRIP und PHASE zum Plotten bestimmter Ausgabegrößen. Der Benutzer kann jetzt durch kurze Befehle, z.B. des Runtime-Kommandos START die Simulation starten und über die Eingabe von PL, STRIP oder PHASE die Ergebnisse plotten: Kommando Bedeutung START Der Simulationslauf wird gestartet PHASE Das Plotbild gemäß Prozedur PHASE wird erzeugt SET THIC = 1 Ein neuer Anfangswert für TH wird gesetzt START Ein neuer Simulationlauf wird gestartet usw. Weitere Simulationsläufe des Pendels und die Arbeit mit verschiedenen Runtime- Kommandos zeigt ab Seite 137 das ACSL im Anhang 4.3. Im nächsten Kapitel wird eine weitere Simulationssprache behandelt, die sehr benutzerfreundlich aufgebaut ist. 2.4 Matlab Simulink MATLAB Simulink (The MathWorks Inc., Natick, MA, USA) bietet eine weitere Möglichkeit, dynamische Systeme zu simulieren und so das dynamische Verhalten zu beschreiben. Zerlegt man das zu untersuchende System in einzelne Blöcke, so ermöglicht Simulink die benutzerfreundliche Programmierung unter Verwendung einzelner, in Bibliotheken vorhandener Programmbausteine. <?page no="77"?> 63 Die prinzipielle Arbeitsweise kann sehr gut am Beispiel eines Spannungsmessinstrumentes gezeigt werden. 2.4.1 Beispiel Spannungsmessinstrument Das Instrument gemäß Abbildung 50 auf Seite 63 misst die am Gerät angelegte Spannung u und lässt einen Zeiger proportional zu u um einen Winkel ausschlagen. Es setzt sich aus einem elektrischen und einem mechanischen Teil zusammen. 2.4.1.1 Modellbildung und Wirkungsplan Der elektrische Teil besteht aus einem Magneten und einer Spule mit dem Ohmschen Widerstand R. Legt man eine Spannung u an das Gerät, fließt durch die Spule ein Strom i, der im Magnetfeld zu einem zu i proportionalen elektrischen Moment m führt (Faktor k ). Dieses Moment lenkt den Zeiger aus. Der Zeiger wird verdreht und dreht dabei auch die Spule durch das Magnetfeld des Magneten. Dadurch wird in der Spule eine der Spannung u entgegengesetzte Spannung u i induziert, die proportional zur Drehgeschwindigkeit α ist (Faktor k i ). Der Spannungsabfall an der Spule ist in jedem Zeitpunkt t proportional zu u u i , der Strom i entsprechend kleiner als zu Anfang bei t = 0. Dies führt dazu, dass die Zeigerbewegung gedämpft wird und der Zeiger nach einer gewissen Zeit zur Ruhe kommt, wenn die zu messende Spannung u konstant ist. Der mechanische Teil besteht aus einem Zeiger mit der Trägheitsmasse J. Wird er um den Winkel verdreht, spannt er eine Drehfeder mit der Drehfedersteifigkeit C F . Dabei entsteht ein mechanisches Moment m F , das der Bewegung entgegenwirkt. Wenn die beiden Momente m und m F im Gleichgewicht stehen und der Zeiger abgebremst ist, hat er seinen Endausschlag erreicht. Abbildung 50: Spannungsmessgerät u i Drehfeder (C F ) J Spule (R, k i , k ) Magnet <?page no="78"?> 64 Das diesen Vorgang beschreibende Gleichungssystem lautet: Mechanik Zeiger-Bewegung Rückstellmoment der Feder rad m kg 10 64 , 0 J mit m m J 2 5 F − Φ ⋅ = − = α rad m N 10 6 C mit C m 4 F F F − ⋅ = α ⋅ = Elektrik Vom Strom erzeugtes Maschenlauf und induzierte Moment: Spannung: A m N 10 8 K mit i K m 2 − Φ Φ Φ ⋅ = ⋅ = Ω ⋅ = = α ⋅ − = − = ⋅ 3 i i i 10 2 R und A s V 2 , 1 K mit K u u u i R Wirkungsplan: • Im Wirkungsplan werden grafisch die Wirkzusammenhänge der einzelnen Signale dargestellt. • Signale bezeichnen den Wert einer physikalischen Größe; im Wirkplan werden sie als Pfeile gekennzeichnet. • Verzweigungspunkte kennzeichnen einen Signalabgriff • Verknüpfungen von Signalen durch Summenpunkte (Addierer, Multiplizierer,..) • Übertragungsblöcke kennzeichnen das dynamische Verhalten des Systems; optische Kennung durch charakteristische Sprungantworten, …). Für das Beispiel des Spannungsmessers zeigt Abbildung 51 den Wirkungsplan: Abbildung 51: Wirkungsplan des Spannungsmessers 1/ R K u m i u i 1/ J C F K i 1 1 m F +- α α α <?page no="79"?> 65 Damit ist die Systemstruktur festgelegt, und im Wirkungsplan wird deutlich, wie die einzelnen Teile zusammenwirken. Im nächsten Schritt muss der Wirkungsplan umgesetzt werden in ein lauffähiges Simulink-Programm. Dies geschieht durch grafisches Editieren. 2.4.1.2 Grafisches Editieren des Schaltplans Nach Anwählen des Programms Matlab Simulink öffnet sich das Matlab-Command- Window, über das man zur Erstellung eines neuen Modell-Files gemäß Abbildung 52 vorgeht. Abbildung 52: Matlab Command Window Es öffnet sich ein Fenster, in dem man das zu untersuchende Modell aufbauen kann, Abbildung 53 auf Seite 66. Die Komponenten (Blöcke) werden mit Drag & Drop aus der Simulink Library geholt. Die Blöcke werden durch Doppelklick mit der linken Maustaste geöffnet, siehe Abbildung 54 auf Seite 66. In die dafür vorgesehen Felder werden die individuellen Blockparameter durch Überschreiben der Standardwerte eingetragen und durch Betätigen der OK-Taste bestätigt. Alle Blöcke können mit Namen versehen werden, zur optischen Hervorhebung auch mit einem Schatten oder mit wählbarer Linienfarbe bzw. Einfärbung, usw. versehen werden. Außerdem gibt es die Möglichkeit, im Arbeitsbereich durch Klicken auf eine gewünschte Position Texte zur Erläuterung einzufügen. Darüber hinaus empfiehlt es sich, bei Bedarf Blöcke zu gruppieren (Bildung von Subsystemen), um bei komplexen Systemen die Übersichtlichkeit der Struktur zu gewährleisten. <?page no="80"?> 66 Abbildung 53: Aufbau des Spannungsmesser-Programms in MatlabSimulink Abbildung 54: Beispiel eines geöffneten Gain-Blocks 0.08 <?page no="81"?> 67 Das vollständige Simulink-Programm des Pendels gemäß Wirkplan von Abbildung 51 zeigt Abbildung 55. Abbildung 55: Simulink-Programm des Spannungsmessers Die an den Spannungsmesser angelegte Spannung von z.B. 10 V wird über den Block „Step“ aus der Bibliothek „Sources“ als Einheitsprung in das Programm gegeben. Der Zeitverlauf des Ausschlags in Grad bzw. in Radiant kann unter Verwendung des Blockes „Scope“ aus der Bibliothek „Sinks“ aufgezeichnet werden. Man sieht: Die Umsetzung eines Differenzialgleichungssystems in ein lauffähiges Simulationsprogramm ist mit Matlab Simulink sehr einfach. In der Programmbibliothek, dem Simulink Library Browser, stehen eine Vielzahl unterschiedlichster Blöcke zur Verfügung. Entscheidend für die näherungsweise Berechnung des Zeitverhaltens ist die richtige Auswahl des Integrationsverfahrens, bevor die Simulation gestartet wird. 2.4.1.3 Integrationsverfahren Es stehen verschiedene Simulationsprogramme zur Verfügung. Je nach Anwendungsfall werden gemäß Tabelle 6 unterschiedliche Methoden empfohlen: Euler Runge-Kutta Gear Adams Linsim Referenz ++ + - o Linear + + o o ++ Nichtlinear - ++ ++ ++ o Diskontinuierlich + + - o Diskrete Anteile + ++ - o Steife Systeme - - ++ - + Tabelle 6 Integrationsverfahren <?page no="82"?> 68 • Es sollten immer mehrere Integrationsverfahren auf ein Modell angewandt werden. Stimmen die Ergebnisse überein, kann der Algorithmus ausgewählt werden, der optimal ist, z.B. die kürzeste Rechenzeit benötigt. • Sind Ergebnisse nicht reproduzierbar, deutet dies auf numerische Instabilitäten hin. Daher immer mehrere Simulationsläufe durchführen. • Ebenfalls empfehlenswert sind Referenzsimulationen mit unterschiedlichen Schrittweiten. Wichtig ist, dass man die Simulationsendezeit festlegt, bei der der Simulationslauf abgebrochen wird. Damit ist alles vorbereitet, und die Simulation kann gestartet werden. 2.4.1.4 Durchführung der Simulation Nach Festlegung eines Endzeitpunkts für die Simulation von z.B. 3 Sekunden kann der Simulationslauf gestartet werden. Das Ergebnis der Simulation ist in Abbildung 56 dargestellt. Abbildung 56: Ergebnisdarstellung Sprungantwort Spannungsmesser Es zeigt sich, das Voltmeter wird wie erwartet ausgelenkt und erreicht nach ungefähr 1 Sekunde gut gedämpft einen stationären Wert von ( = 0,67 rad ) = 38 Grad. Damit ließe sich eine Skala in Volt über dem Messgerät anbringen, die bei einem Winkel = 38 Grad eine Spannung von 10 Volt anzeigt. Der Umrechnungsfaktor für die vorliegende Konfiguration wäre also 10 V / 38 Grad. <?page no="83"?> 69 2.4.2 Beispiel Gedämpftes nichtlineares Pendel Das Beispiel des nichtlinearen Pendels von Abschnitt 2.3.2 soll hier nochmals mit MatlabSimulink behandelt werden. Zur Lösung der Differenzialgleichung wird das in Abbildung 57 dargestellte Programm aufgebaut. Als Anfangsbedingung wird vorausgesetzt, dass das Pendel bewegungslos nach unten hängt und angestoßen wird: Lage: TH(t=0) = TH0 = 0, Anfangsgeschwindigkeit: THp(t=0) = THp0 = 10 rad / s Der Zeitverlauf des Pendelausschlags TH, der Geschwindigkeit THp und der Beschleunigung THpp wird beispielhaft mit einem Scope aufgezeichnet. Eine weitere Darstellungsmöglichkeit ergibt sich unter Verwendung von XY-Graph: Hier ist die Geschwindigkeit THp des Pendels über dem Ausschlagwinkel TH geplottet. Unabhängig vom Simulationstool, also auch beim Einsatz von Matlab Simulink, ist erneut sehr wichtig, dass der Verlauf der Kurven in allen Details verstanden wird, dass die Grafik physikalisch interpretiert wird. Deckt sich alles mit den physikalischen Vorstellungen, wächst das Vertrauen in das Modell und die Gültigkeit der Ergebnisse. Zeigt der Verlauf der Kurven einen Widerspruch zur Physik, muss geprüft werden, ob bei der Entwicklung des Simulationsmodells Fehler gemacht wurden, die dann zu korrigieren sind. Der Simulink-Kurvenverlauf deckt sich mit dem mit ACSL berechneten. Abbildung 57: Nichtlineares Pendel mit Anfangsbedingung TH0=0.0 und THp0=10.0 <?page no="84"?> 70 Matlab Simulink bietet sehr viele unterschiedliche Bausteine zur Simulation dynamischer Systeme. Die Sprache ist benutzerfreundlich aufgebaut und findet weltweit zunehmende Verbreitung, auch in der Automobilindustrie. Für einen etwas tieferen Einblick in die Möglichkeiten von Matlab Simulink soll abschließend ein weiterens, etwas komplexeres Beispiel dargestellt werden. Es geht dabei um die Vertikaldynamik eines Fahrzeugs. Nach Beschreibung der Physik des Fahrzeugs wird in Anlehnung an die Empfehlungen von Kapitel 1 ein mathematisches dynamisches Modell entwickelt (Grenzen und Detaillierungsgrad). Zum Schluss wird das beschreibende System von Differenzialgleichungen in Matlab Simulink umgesetzt. 2.4.3 Vertikaldynamik eines Zweimassen-Fahrzeugs 1. Physik: Bei der Beschreibung der Vertikaldynamik eines Fahrzeugs während einer Fahrt über eine schlechte Straße spielen die beiden Gesichtspunkte Fahrkomfort und Fahrsicherheit eine bedeutende Rolle. Einerseits geht es darum, den Fahrkomfort zu maximieren, d.h., die auf den Fahrer wirkenden Vertikalbeschleunigungen zu minimieren. Andererseits soll das Rad nicht den Kontakt zur Straße verlieren und springen. Dies hätte sonst z.B. einen Verlust an Seitenführungskraft während Kurvenfahrt zur Folge. Die dynamischen Radlastschwankungen müssen folglich unter dem Aspekt der Fahrsicherheit klein gehalten werden. Zur Untersuchung dieser Fragestellung reicht es aus, das Fahrzeug als Zweimassenschwinger zu modellieren, also auf weitere Bewegungsfreiheitsgrade des Aufbaus zu verzichten, z.B. Nicken (beim Bremsen), Wanken (in der Kurve), Drehfreudigkeit um die Hochachse (zur Charakterisierung des fahrdynamischen Verhalten, Über- oder Untersteuern). Das entsprechend einfache Modell eines „Viertel-Fahrzeugs“ zeigt Abbildung 58 auf Seite 71. Es besteht aus zwei starren Körpern, dem Aufbau mit der Masse m_A und nur einem Rad mit der Masse m_R, die beide in vertikaler Richtung auf und ab schwingen können. Die Bewegungskoordinaten sind x_A und x_R und deren zeitlichen Ableitungen. Bindeglied zwischen Rad und Aufbau sind • eine Aufbaufeder mit der Federsteifigkeit c_A zur Übernahme der statischen und dynamischen Last des Aufbaus und • ein hydraulischer Stoßdämpfer mit einer geschwindigkeitsproportionalen Dämpfungskonstanten d_A zur Bedämpfung der Aufbau- und Radschwingungen. Auf eine in den Fahrzeugen übliche unterschiedliche Zug- und Druckstufe der Dämfung wird in diesem einfachen Modell zunächst verzichtet. Das Rad selbst ist über den Reifen zur Straße hin elastisch gefedert (Reifenfedersteifigkeit c_R), die Dämpfung im Reifen ist vernachlässigbar. Das Modell erlaubt eine Anregung durch die Straße der Größe x_S. Außerdem kann eine Zusatzlast F_ZL berücksichtigt werden, zur Simulation einer veränderlichen Beladung des Fahrzeugs oder einer erhöhten Radkraft am kurvenäusseren Rad bei Kurvenfahrt. 2. Mathematisches Modell: Zur Herleitung des mathematischen Modells werden die beiden Massen freigeschnitten und alle angreifenden (äußeren) Kräfte eingetragen, die zu einer Beschleunigung xpp_A und xpp_R der beiden Massen führt. Die Federkraft F_F ist proportional zum Relativweg (x_R - x_A), die Dämpferkraft F_D <?page no="85"?> 71 proportional zur Relativgeschwindigkeit (xp_R - xp_A) und die Reifenkraft F_R proportional zur Reifeneindrückung (x_S - x_R). Abbildung 58: Modell des Zweimassenschwingers Aus Gleichgewichtsbedingungen (nach Newton oder d’Alembert) folgt für die dynamischen Kräfte (mit xp=dx/ dt und xpp=dxp/ dt, Ableitung von x bzw.xp nach der Zeit): m_A * xpp_A = F_F + F_D - F_ZL m_R * xpp_R = - F_F - F_D + F_R Die statischen Lasten durch das Gewicht des Fahrzeugs fallen auf der rechten Seite des Gleichungssystems heraus und werden deshalb gar nicht als weitere Summanden erwähnt. F_F_stat - m_A * g = 0 und F_R_stat - (m_A * g + m_R * g) = 0. Beachte: Die Radaufstandskraft F_R_gesamt = F_R + F_R_stat kann nicht negativ werden, das Rad „klebt“ ja nicht an der Straße, sondern hebt ab, wenn die Radaufstandskraft Null wird. Setzt man die Beziehungen für die Kräfte ein, folgt m_A * xpp_A = c_A * (x_R - x_A) + d_A * (xp_R - xp_A) - F_ZL m_R * xpp_R = c_A * (x_R - x_A) d_A * (xp_R - xp_A) + c_R * (x_S - x_R) mit F_R_gesamt = F_R + F_R_stat 0. F_R Reifenkraft Rad m_R Aufbau m_A Feder c_A Dämpfer d_A x_S x_R x_A Reifen c_R F_ZL Modell Rad m_R Aufbau m_A F_ZL Zusatzlast F_F F_D Feder-/ Dämpfer- Kraft Kräfte x_A xp_A xpp_A x_R xp_R xpp_R Straßenanregung <?page no="86"?> 72 Es ergibt sich also ein gekoppeltes System von linearen Differenzialgleichungen vierter Ordnung. Für die beiden Massen Aufbau und Rad muss jeweils aus der Beschleunigung xpp durch zweimalige Integration die Geschwindigkeit xp und der Weg x berechnet werden. Um das Modell einfach zu halten, erfolgt bei der Berechnung dieser Beschleunigungen zunächst keine Berücksichtigung der Bedingung für die Radaufstandskraft F_R_gesamt 0. Wenn sich jedoch zeigen sollte, dass im Rahmen spezieller Manöver diese Bedingung verletzt wird, muss das mathematische Modell erweitert werden. In allen anderen Fällen sind die Simulationsergebnisse absolut korrekt. 3. Modellparameter: Es fehlt jetzt nur noch die Wahl geeigneter physikalischer Parameter für das Modell und Anfangsbedingungen für den Zeitpunkt t = 0 für jeden der vier später im Simulink-Modell verwendeten Integratoren. Für eine erste Variante soll gelten: Physikalische Modellparameter: m_A = 400 kg m_R = 50 kg c_A = 20000 N/ m d_A = 1400 N/ (m/ s) c_R = 200000 N/ m F_ZL = - 2000 bis + 2000 N Die Anfangsauslenkung bzw. die Anfangsgeschwindigkeit der beiden Massen sei x_A_0 = x_A (t=0) = 0, x_R_0 = x_R (t=0) = 0, xp_A_0 = xp_A (t=0) = 0, xp_R_0 = xp_R (t=0) = 0, das Fahrzeug steht also in Konsstruktionslage und befindet sich anfangs in Ruhe. Damit ist das mathematische Modell beschrieben und kann in Matlab Simulink umgesetzt werden. 4. Aufbau des Simulink-Modells: Es empfiehlt sich zunächst, jede der beiden Differenzialgleichungen nach der höchsten Ableitung aufzulösen hier also durch einfache Division durch die jeweiligen Massen nach den beiden Beschleunigungen xpp_A und xpp_R: xpp_A = 1/ m_A * ( c_A * (x_R - x_A) + d_A * (xp_R - xp_A) - F_ZL) xpp_R = 1/ m_R * (c_A * (x_R - x_A) d_A * (xp_R - xp_A) + c_R * (x_S - x_R)) Damit ist festgelegt, wie groß die Beschleunigungen der Massen zu jedem Zeitpunkt t ist. Es kann daraus durch Integration die Geschwindigkeit xp_A(t) und xp_R(t) berechnet werden und anschließend die Auslenkung x_A(t ) und x_R(t). Der Aufbau des Modells beginnt mit der Anordnung von vier Integrator-Böcken aus dem Kapitel Continuous des Simulink Library Browsers. Die Ausgänge der vier Intagratoren dienen zur Berechnung der Kräfte, die zu Beschleunigungen an den beiden Massen führen. Die berechneten Kräfte werden deshalb den beiden Differenzialgleichungen entsprechend auf die Eingänge der beiden Integratoren <?page no="87"?> 73 rückgeführt, die die Geschwindigkeiten xp_A und xp_R berechnen. So entsteht ein Simulink-Modell ZMS gemäß Abbildung 59: Abbildung 59: Simulink-Modell des Zweimassenschwingers Außerdem sind auch gleich 3 Scopes aus dem Kapitel Sinks in das Simulink-Modell eingebaut zur graphischen Ausgabe der interessierenden Größen x_Aufbau, x_Rad und des Relativweges x_Relativ = x_R - x_A zwischen Rad und Aufbau. Man sieht in Abbildung 59, wie selbst im Falle dieses einfachen Modells des Zweimassenschwingers viele Blöcke notwendig sind, um die Differenzialgleichungen zu programmieren. Komplexere Modelle mit vielen Bewegungsfreiheitsgraden ergeben sehr schnell eine etwas unübersichtliche Darstellung. Matlab Simulink bietet deshalb die Möglichkeit, Modellteile in Form so genannter Subsystems zusammenzufassen. 5. Subsystems: Für das vorliegende Beispiel bietet sich an, die Berechnungen für z.B. den Aufbau oder das Rad in einem Subsystem darzustellen. Hierzu werden die Blöcke markiert, die im Subsystem enthalten sein sollen, siehe die gestrichelt eingerahmte Umgebung links in Abbildung 60. Über den Edit-Befehl Create Subsystem wird das gewünschte Subsystem erzeugt, siehe Ausschnitt des Modells rechts: Abbildung 60: Create Subsystem Create Subsystem <?page no="88"?> 74 Alle Linien, die in Abbildung 60 auf Seite 73 vom gestrichelten Rahmen des Subsystems, linkes Bild, geschnitten werden (3 Eingänge, 2 Ausgänge), finden sich rechts als Eingänge In1 bis In3 bzw. Out1 und Out2 wieder. Werden in gleicher Weise auch Subsysteme für das Rad, das Federbein (Feder+Dämpfer) und den Reifen eingeführt, so ergibt sich nach etwas grafischer Überarbeitung folgendes Modell ZMS_Block mit den vier zugehörigen Subsystemen, Abbildung 61: Abbildung 61: Simulink Modell des Zweimassenschwingers mit Subsystems Ausgabe: Die vielen Ausgabegrößen werden in einem weiteren Subsystem Ausgabe in verschiedenen Scopes dargestellt. Im Beispiel werden diese Größen aber nicht über Signalleitungen in das Subsystem Block Ausgabe übergeben. Der Übersichtlichkeit halber erfolgt die Übergabe über Signal Routing jeweils mit Goto- <?page no="89"?> 75 (Sender) und From-Blöcken (Empfänger), z.B. die markierte Größe x_A aus dem Subsystem Aufbau in das Subsystem Ausgabe, Abbildung 62: Abbildung 62: Möglichkeiten des Signal-Routing (Goto From) 6. Manöver 1: Ein erster Simulationslauf soll das Verhalten zeigen bei einer Straßenanregung im Sinne einer Sprungfunktion der Höhe 0,01 m zum Zeitpunkt t = 1s. Ein Manual Switch aus Kapitel Signal Routing, siehe Abbildung 63, erlaubt bei Doppelklick ein Umschalten auf eine andere Anregung. Die Zusatzlast auf den Aufbau hat die Größe 1000N. Das entsprechende Modell zeigt Abbildung 63: Abbildung 63: Endgültige Fassung des Modells Zweimassenschwinger <?page no="90"?> 76 Simulation Manöver 1: Wird die Simulation mit F_ZL = 0N gestartet, rechnet das Programm das Verhalten des Zweimassenschwingers aus, bis der Simulationsendzeitpunkt erreicht ist, hier t = 5 s. Die Reaktion des Systems ist wie erwartet, die beiden Massen folgen der Straße und schwingen auf einer Höhe von 0,01 cm ein, Abbildung 64 : Abbildung 64: Manöver 1: Reaktion mit schwacher Dämpfung d_A = 1400N/ (m/ s) Das Rad folgt wegen der steifen Reifenfeder sehr schnell, der Aufbau wegen der 10fach höheren Masse und der deutlich weicheren Aufbaufeder träger, der Relativweg (x_R x_A) nimmt zu Anfang entsprechend positive Werte an, d.h., das Rad federt ein. Am Ende ist der Relativweg wieder Null. Der Übergang ist allerdings mit einer Dämpfung von d_A = 1400N/ (m/ s) nur sehr schwach gedämpft. Das Fahrzeug sollte nach einem einmaligen Überschwingen recht schnell zur Ruhe kommen. Die Simulation wird deshalb mit einem höheren Wert von d_A wiederholt, d_A = 3000N/ (m/ s). Jetzt zeigt das Fahrzeug zufriedenstellendes Verhalten, Abbildung 66 auf Seite 77. Zur Kontrolle wird zuvor die Radaufstandskraft auf Radabheben geprüft. Die dynamische Radlaständerung von maximal -920 N bleibt jedoch immer deutlich kleiner als die statische Kraft von F_stat = (m_A+m_R) * g = 4415 N, Abbildung 65 auf Seite 77. 7. Manöver 2: In einem zweiten Manöver wird das Verhalten des Zweimassenschwingers berechnet, wenn eine Zusatzlast F_ZL = 1000 N auf den Aufbau wirkt bei Fahrt über eine ebene Straße, x_S = 0. Der Block-Parameter Final Value des Step- Blockes wird deshalb auf Null gesetzt, der Blockparameter für die Konstante Zusatzlast auf den Wert F_ZL = 1000N. x_R x_A x_S x_Relativ <?page no="91"?> 77 Abbildung 65: Manöver 1: Dynamische Radlaständerung bei Sprungantwort Abbildung 66: Manöver 1: Sprungantwort mit starker Dämfung´d_A = 3000 N/ (m/ s) Durch diese Kraft werden beide Federn zusammengedrückt. Die Reifenfeder verkürzt sich um x_R = F_ZL / c_R = 1000 N / (200000 N/ m) = 0,005 m, die Aufbaufeder um F_ZL / c_A = 1000 N / (20000 N/ m) = 0,05 m. Damit senkt sich der Aufbau um x_A = 0,05 m + 0,005 m = 55 mm ab, siehe Abbildung 67 auf Seite 78, links. Auch die in der rechten Abbildung dargestellten Kräfte verhalten sich, wie es zu erwarten ist: sie wachsen proportional zur Federverkürzung auf einen stationären Endwert von 1000N an; die Dämpferkraft F_D wirkt nur so lange, wie die Massen in Bewegung sind, im eingeschwungenen Zustand ab ca. t = 1s ist sie wieder 0. x_R x_A x_S x_Relativ x_S - 920 N F_R <?page no="92"?> 78 Abbildung 67: Manöver 2: Reaktion auf eine Zusatzlast F_ZL = 1000 N 8. Manöver 3: Abschließend wird simuliert, wie das Fahrzeug auf eine rauhe Straße reagiert (Sinus-Anregung mit einer Amplitude von 2 mm und einer Erregerfrequenz von f_err= 10 Hz), bei gleichzeitiger Lasterhöhung, z.B. wegen Kurvenfahrt, von F_ZL = 1000 N. Durch Doppelklick auf den Manual Switch wird umgeschaltet auf die Sinusanregung durch den Block Sine Wave aus Sources, danach die entsprechenden Blockparameter eingestellt, Amplitude auf 0,002 und Frequency auf 2*pi*10. Das Ergebnis zeigt Abbildung 68: Der Aufbau führt einerseits eine Eigenschwingung aus und erreicht gut gedämpft seine neue Lage von x_A= - 0,055 m. Abbildung 68: Manöver 3: Reaktion auf Zusatzlast + 10 Hz-Anregung der Straße x_R x_A x_Relativ x_S F_F F_D F_R x_R x_A x_Relativ x_S x_A 10 Hz - 55 ± 0,54 mm <?page no="93"?> 79 Er schwingt einerseits mit der Eigenfrequenz f_0 = 1/ (2 pi) * WURZEL(m_A/ c_A) = 1/ (2 pi) * WURZEL(400/ 20000) = 1,13 Hz. Andererseits reagiert er auf die Anregung durch die Straße und führt eine erzwungene Bewegung aus (eine der Eigenschwingung überlagerte Schwingung). Er schwingt mit f_err=10 Hz, allerdings nicht mit 2 mm Amplitude, sondern bei einem Verhältnis der Frequenzen von f_err/ f_0=8,9 mit einer viel kleineren Amplitude. Das Verhältnis der Amplituden beträgt in diesem Fall (0,54/ 2)mm/ 2mm = 0,135, der relativ träge Aufbau kann also der Straße bei dieser Anregung nicht folgen; das hat hohen Komfort für die Passagiere zur Folge. Als nächstes soll am Beispiel des Zweimassenschwingers auf zwei weitere Möglichkeiten eingegangen werden, die das Arbeiten mit Matlab Simulink benutzerfreundlich macht: das Handling mit Modellparametern und das Abspeichern von Werten in der Matlab Workspace. 9. Simulink-Modellparameter: Bisher wurden im Modell des Zeimassenschwingers jeder Parameterwert des Modells als Konstante in den Blockparametern eingetragen, z.B. die Aufbaumasse mit einer Größe von m_A = 400kg als „400“. Bei komplexeren Modellen mit sehr vielen Blöcken und Parametern kann eine Änderung der Aufbaumasse auf beispielsweise m_A = 500N sehr umständlich sein: vielleicht kommt dieser Wert in mehreren Blöcken vor und muss deshalb überall von Hand umgesetzt werden. Um solchen Schwierigkeiten zu entgehen, erweist es sich als viel praktischer, eine Variable für die Masse des Aufbaus zu definieren, z.B. m_A, und dann den Variablennamen als Blockparameter einzutragen. Ändert man entsprechend auch alle anderen Werte um, entsteht das modifizierte Modell ZMS_Block_Parameter gemäß Abbildung 69 auf Seite 80. Die unterschiedlichen Variablen finden sich in Subsystemen wieder, wie z.B. „Aufbau“ oder „Federbein FB“, aber auch im Source-Block für die konstante Zusatzlast F_ZL, und den Function-Block Parameters für den Integrator zur Integration der Radbeschleunigung. Hier wurde der Anfangswert der Radgeschwindigkeit als xp_R_0 eingetragen. Vor dem Start der Simulation müssen allerdings diesen neu definierten Variablen die konstanten Zahlenwerte zugewiesen werden. Dies geschieht am besten durch Aufruf eines so genannten m-Files, wie er in Abbildung 70 auf Seite 80, mit dem Namen ZMS_Parameter gezeigt wird. Dieser File muss im Matlab Command Window aufgerufen werden. Nach der Durchführung sind in der Matlab Workspace alle Variablen bekannt und stehen für die Simulation zur Verfügung, Abbildung 71 auf Seite 81. Eingaben hinter >> werden von Matlab als Kommando verstanden und ausgeführt. Die Abfrage who liefert als Antwort alle Variablen, die in der Workspace vorhanden sind. Die erste Abfrage zeigt: noch gibt es keine Variablen. Die Abfrage what zeigt als Ergebnis alle vorhandenen m-Files, hier also nur den File ZMS_Parameter , aber auch alle mdl-Files, hier also • ZMS, das Modell gemäß Abbildung 59 auf Seite 73, • ZMS_Block, das Modell gemäß Abbildung 63 auf Seite 75 und • ZMS_Block_Parameter , das Modell gemäß Abbildung 69 auf Seite 80. <?page no="94"?> 80 Abbildung 69: Ersatz konstanter Werte durch wählbare Variablen Abbildung 70: m-File zur Zuweisung gewünschter Werte an alle Systemparameter <?page no="95"?> 81 Abbildung 71: Auszug aus dem Matlab Command Window Der Aufruf von ZMS_Parameter weist allen Variablen die gewünschten Zahlenwerte zu. Danach sind alle Variablen bekannt. Dies zeigt eine erneute Abfrage who. Alle zwölf Variablen des m-Files werden aufgelistet. Und die Abfrage von z.B. m_A beweist, dass der gewünschte Zahlenwert von „400“ tatsächlich der Aufbaumasse m_A zugewiesen wurde. 10. Abspeichern und Weiterverarbeiten von Ergebnissen in Simulink: Es gibt die Möglichkeit, den Zeitverlauf der Zustandsgrößen in Vektoren abzuspeichern. Damit wird es möglich, die Werte nach der Simulation weiter zu verarbeiten oder interessierende Größen in eigenen Plot-Routinen grafisch darzustellen. Anwendungsfall 1: Dies soll am Beispiel der Radaufstandskraft F_R_gesamt = F_R_statisch + F_R gezeigt werden. Zuerst muss die Ausgabegröße F_R_gesamt im Subsystem Ausgabe erzeugt werden, Abbildung 72 auf Seite 82. Zur dynamischen Radlast F_R wird über den Function-Block Add aus Math Operations die statische Last addiert und anschließend in einem Scope grafisch abgebildet. Neu ist jetzt, dass die Kraft F_R_gesamt mit Hilfe des Sink-Blocks To Workspace in ein spezifiziertes Feld geschrieben werden kann, hier in den Vektor F_R_gesamt Dieser Vektor steht nach der Simulation in der Main Workspace von Matlab zur Verfügung. Um diesen Vektor über der Zeit zu plotten, muss auch die Zeit als Vektor erzeugt werden. Dies geschieht über den Source-Block Clock. Die einzelnen Zeitwerte werden analog zu F_R_gesamt im Vektor Time abgespeichert. Nach erfolgter Simulation kann somit die Radaufstandskraft F_R_gesamt über der Zeit Time geplottet werden. Die bildliche Darstellung kann über einzelne Befehle individuell gestaltet werden. Speichert man diese Befehlsfolge in einem weiteren m-File ab, kann das Bild immer wieder neu durch Aufruf dieses m-Files erzeugt werden. Das m-File Plot und die bei <?page no="96"?> 82 Aufruf entstehende Figure1 ist in Abbildung 73 dargestellt. Farbe und Strichstärke haben die aktuellen Default-Parameter, können aber über entsprechende Parameter auch individuell gestaltet werden, genauso wie die Skalierung der Achse, etc. Abbildung 72: Erzeugung der Ausgabegröße Radaufstandskraft Abbildung 73: Grafische Ausgabe der Radaufstandskraft Anwendungsfall 2: Ein weiteres Beispiel, mit den abgespeicherten Werten zu arbeiten, wird abschließend gezeigt. So kann es interessant sein, die Ergebnisse aus verschiedenen Simulationsläufen in einem gemeinsamen Bild darzustellen. Das Manöver 1 wurde mit zwei verschiedenen Werten für die Dämpfung des Aufbaus simuliert. Es war d_A = 1400 bzw. 3000 Ns/ m. Die Ergebnisse sind in Abbildung 64 auf Seite 76 und Abbildung 66 auf Seite 77 dargestellt. Sollen die Bewegungen des <?page no="97"?> 83 Aufbaus x_A für beide Fälle gemeinsam dargestellt werden, müssen diese Größen z.B. als x_A_1400 und x_A_3000 abgespeichert werden mit dem bereits erwähnten Sinks-Block To Workspace. Danach können sie geplottet werden, Abbildung 74: Abbildung 74: Sprungantwort des Aufbaus mit unterschiedlicher Dämpfung 11. Modellerweiterung Radabheben: Bisher gab es bei der Simulation des Zweimassenschwingers eine einschränkende Annahme: Das Rad darf nicht abheben von der Straße. Dies war auch nicht der Fall bei den behandelten Manövern • Manöver 1: Sprungantwort des Systems auf einen Sprung der Straße von 0,01m Höhe bei unterschiedlicher Aufbaudämpfung d_A, • Manöver 2: Reaktion des Systems auf plötzliche Zuladung in Form von z.B. einer Zusatzlast F_ZL in Höhe von 1000N und • Manöver 3: Reaktion des Systems auf plötzliche Zuladung in Form von z.B. einer Zusatzlast F_ZL in Höhe von 1000N bei gleichzeitig wirkender rauher Straße in Form einer Sinusanregung von 10Hz mit einer Amplitude von 0,002m. Nun soll das Modell erweitert werden um den Aspekt Radabheben. Kontakt zwischen Rad und Straße bedeutet, dass die Radaufstandskraft F_R_gesamt > 0 d_A = 1400 N/ (m/ s) d_A = 3000 N/ (m/ s) <?page no="98"?> 84 ist. Im Falle Radabheben wird diese Kraft gleich Null. Die Kraft kann nicht negativ werden, denn der Reifen „klebt“ ja nicht an der Straße fest. Die Radaufstandskraft F_R_gesamt, siehe Abbildung 72 auf Seite 82, setzt sich zusammen aus der der statischen Last und der dynamischen Radlaständerung F_R. Es ist F_R_gesamt = (m_A+m_R) * g + F_R 0. mit F_R = c_R * (x_S x_R). Die Forderung nach F_R_gesamt 0 bedeutet also, dass die dynamische Radlaständerung F_R - F_R_stat = - (m_A + m_R) * g sein muss. Dies muss im Simulink Modell umgesetzt werden, um nicht weiterhin beschränkt zu sein bei der Wahl der Manöver. Im Subsystem Reifen, Abbildung 75, wird ein Saturation-Block aus Discontinuities eingebaut und die zugehörigen Block-Parameter gesetzt: Obere Beschränkung: Upper limit = inf, also ohne Beschränkung Untere Beschränkung: Lower limit = - (m_A + m_R) * g Durch diese Erweiterung ist das Modell geeignet für Manöver, bei denen das Rad zeitweise von der Straße abhebt. Dies soll in einem weiteren Simulationslauf gezeigt werden unter Verwendung des Modells ZMS_Block_Parameter_Radabheben. Abbildung 75: Möglichkeit von Radabheben 12. Manöver 4: Es wird wie bereits bei Manöver 1 erneut die Sprungantwort des Systems untersucht, dieses Mal jedoch auf einen Sprung der Straße nach unten, entprechend einer Fahrt vom Bordstein herunter. Nach 0,1s soll der Sprung erfolgen, die Sprunghöhe soll -0,15m betragen. Die Radaufstandskraft ist in Abbildung 76 auf Seite 85 (rechts) dargestellt. Sie wird zum Zeitpunkt t = 0,1s (Bordsteinabfahrt) schlagartig Null. Das Rad fällt nach unten und trifft erst zum Zeitpunkt t = 0,199s wieder auf die Straße, die Radaufstandskraft wird ab diesem Moment wieder größer Null. Das Rad (im Bild links) folgt der Straße wegen der Trägheit der Radmasse nur verzögert. Es kommt deshalb zum Abheben, da die statische Reifeneindrückung <?page no="99"?> 85 x_R_stat = (m_A+m_R)*g / c_R = 0,022m < 0,15m beträgt. Das Rad fällt nach unten und berührt die Straße erneut ab einer Höhe von x_R= -0,15m + 0,022m zum Zeitpunkt t = 0,199s. Abbildung 76: Sprungantwort auf einen Sprung von -0,15 m 13. Modellerweiterung um eine Zug-/ Druckstufe im Dämpfer: Um den Fahrkomfort zu steigern, ist es üblich, die Druckstufe des Dämpfers weicher zu wählen als die Zugstufe. Im Falle eines Hindernisses auf der Straße, das bei der Überfahrt zu einer Einfederung führt (Druckstufe), soll der Dämpfer eine relativ sanfte Reaktion aufweisen, damit der Fahrzeugaufbau nicht unkomfortabel nach oben geworfen wird. Der Aufbau wird dadurch weniger beschleunigt. In der Zugstufe da-gegen ist eine härtere Dämfung ohne Komfortverlust möglich. Die Berechnung der Dämpfungskraft wird dadurch komplizierter. Die Kraft ist zwar weiterhin proportional zur Dämpfergeschwindigkeit, die Proportionalitätskonstante hängt jedoch vom Vorzeichen der Relativgeschwindigkeit zwischen Rad und Aufbau ab. Es gibt damit zwei Werte für d_A, einen Wert d_A_D für die Druckstufe und einen Wert d_A_Z für die Zugstufe. Das bisher gültige Subsystem Federbein FB, Abbildung 61 auf Seite 74, muss entsprechend umgebaut werden, Abbildung 77: Abbildung 77: Modellerweiterung für Dämpfer mit Zug- und Druckstufe x_R x_A x_S x_Relativ F_R_ges <?page no="100"?> 86 Es wird ein eigenes Subsystem für die Berechnung eingeführt. Eingangsgrößen sind die beiden Geschwindigkeiten xp_R und xp_A. Im Subsystem Dämpfer mit Zug-/ Druckstufe wird zunächst die Relativgeschwindigkeit xp_rel berechnet. Im Anschluss wird unter Anwendung des Logik-Bausteins Compare To Zero aus Logic and Bit- Operations der gerade ungültige Summand mit Null multipliziert. (Dies ist nur eine von mehreren Möglichkeiten, diese Logik in Simulink umzusetzen.) Die Dämpfungskennlinie ist nichtlinear. Abbildung 78 zeigt die Kennlinie F_D in N über der Dämpfergeschwindigkeit xp_rel in m/ s für die beiden Dämpungskonstanten d_A_Druck=1400N/ (m/ s) und d_A_Z=3000N/ (m/ s). Das Ergebnis der Simulation ist in Abbildung 79 auf Seite 87 dargestellt. Der Aufbau ist nach einem Zeitraum von ca. 1,5s ähnlich gut beruhigt wie im Falle der harten Dämfung in Zug- und Druckrichtung mit d_A=3000N/ (m/ s). Das Ziel, den Aufbau ausreichend schnell zu bedämpfen, ist damit erreicht. Das weitere Ziel, den Komfort zu verbessern ohne einen gleichzeitigen Verlust an Dämpfungsleistung, wird auch erreicht. Der Komfort wird im Modell an der Aufbaubeschleunigung xpp_A beurteilt: kleine Beschleunigungen bedeuten hohen Komfort. So zeigt Abbildung 80 auf Seite 87, dass durch die weiche Druckstufe die ersten Beschleunigungsspitze des Aufbaus deutlich reduziert wird, um ca. 33%. Eine Dämpferauslegung mit unterschiedlicher Zug- und Druckstufe ist damit sehr sinnvoll. Abbildung 78: Dämpfungskennlinie mit Zug- und Druckstufe F_D [N] xp_rel [m/ s] Zug Druck 1400N/ (m/ s) 3000N/ (m/ s) <?page no="101"?> 87 Abbildung 79: Sprungantwort des Aufbaus mit drei Dämpfungsvarianten Abbildung 80: Komfort als Funktion der Aufbaubeschleunigung xpp_A Zum Schluss soll noch auf einen Nebeneffekt dieser Dämpferauslegung eingegangen werden. Dieser Effekt ist zwar unerwünscht, führt allerdings aber auch nicht zu einem großen Nachteil. Bei kleiner Druckstufe federt das Rad immer leichter ein als aus. Im Mittel - so lässt sich vermuten - wird das Rad Richtung Aufbau gezogen und der Aufbau wird abgesenkt, weil der Mittelwert der Einfederung > 0 ist. Dies wird durch Simulation 3000 1400 1400/ 3000 - 33% d_A_D/ d_A_Z = 3000 / 3000 d_A_D/ d_A_Z = 1400 / 3000 <?page no="102"?> 88 bewiesen, Abbildung 81. Der Aufbau verliert bei dieser Parameterabstimmung / Anregung ca. 4cm an Höhe. Abbildung 81: Absenkung des Aufbaus bei Sinus-Anregung Allgemein folgt daraus, dass Fahrzeuge bei Fahrt auf schlechten Straßen durch den von Null verschiedenen Mittelwert der Dämpfkraft einfedern und sich der Aufbau somit etwas absenkt. In diesem Buch kann Matlab Simulink nur sehr verkürzt dargestellt werden. Im Anhang 5 ab Seite 149 wird das Modell des Zweimassenschwingers erweitert um ein aktives Federbein, wie es prinzipiell Anwendung findet in dem aktiven Federungssystem ABC (Active Body Control), das in den Top-Baureihen bei Mercedes in Serie bzw. optional angeboten wird. Außerdem erfährt der Leser im Anhang 5 etwas über Demonstrationsbeispiele, die von The Mathworks auf ihrer Internetseite zum näheren Kennenlernen der Sprache Matlab Simulink angeboten werden. Der Leser hat auch die Möglichkeit, auf der Internetseite des Expert-Verlags die im Buch behandelten Matlab Simulink-Beispiele herunterzuladen. Wer eine Lizenz von Matlab Simulink besitzt, kann diese Beispiele selbst anwenden, mit den Modellparametern spielen und so ein besseres Verständnis der zeitkontinuierlichen Simulation entwickeln. x_A x_rel x_S <?page no="103"?> 89 Damit ist das zweite Kapitel über die Simulation zeitkontinuierlicher Systeme abgeschlossen. Das dynamische Verhalten dieser Systeme wird in Differenzialgleichungen beschrieben, die durch schrittweise Integration näherungsweise gelöst werden. Die zwei aus einer Fülle verschiedener Sprachen ausgewählten Beispiele ACSL und Matlab Simulink sollten die Arbeitsweise mit diesen Sprachen zeigen. Der Anwender darf den berechneten Ergebnissen niemals „blind“ glauben. Er muss den Zeitverlauf jeder Größe hinterfragen, die Ergebnisse verstehen. Nur dann, wenn sich die Ergebnisse decken mit der eigenen Vorstellung, steigt das Vertrauen des Anwenders in sein Programm. Wenn nicht, muss er mit dem Simulationsmodell „Experimente“ durchführen, die sein Verständnis vertiefen oder zu Modelländerungen führen. Im folgenden Kapitel wird näher auf die andere Klasse dynamischer Systeme eingegangen, den zeitdiskreten Systemen. <?page no="104"?> 90 <?page no="105"?> 91 3 Simulation zeitdiskreter Systeme 3.1 Einführung 3.1.1 Unterscheidung zeitdiskreter und kontinuierlicher Systeme Die beiden Systemklassen werden in Abhängigkeit des zeitlichen Abstandes unterschieden, in dem die Zustandsänderungen stattfinden. Tendiert der Abstand gegen Null, so spricht man von einem kontinuierlichen System. Hat er dagegen einen bestimmten begrenzten Wert, so handelt es sich um ein zeitdiskretes System. Ein Beispiel für ein kontinuierliches System stellt das einfache 2-Massen-Modell eines Kraftfahrzeuges dar. Die Änderung von Lage und Geschwindigkeit der Fahrzeugmassen verlaufen stetig, auch die Anregungen von der Straße können als quasi stetig angesehen werden. Betrachtet man dagegen den zeitlich veränderlichen Zustand eine Ampelanlage (rot, gelb, grün) oder eines Telefonnetzes (Leitungen belegt bzw. frei), so fällt auf, dass der Systemzustand über einen gewissen Zeitraum konstant ist und nur nach endlichen Abständen geändert wird. Solche Systeme nennt man zeitdiskret. In Abhängigkeit der konkreten Fragestellung kann sein, dass sich ein System unterschiedlich darstellt, zum Beispiel ein Fahrstuhl in einem Hochhaus: • Will man einen Fahrkorb im Schacht so elektrisch ansteuern, dass er zeitoptimal und ruckfrei in eine bestimmte Ebene einfährt, wird man das zugehörige Differenzialgleichungssystem aufstellen unter Berücksichtigung des elastischen Stahlseils und der Dynamik des Elektromotors. Der Fahrstuhl wird also durch ein zeitkontinuierliches System beschrieben. • Will man erreichen, dass der Fahrstuhl optimalen Durchsatz an Fahrgästen hat, muss man eine Ansteuerungslogik entwickeln, die die Wartezeiten der verschiedenen Fahrgäste minimiert. Dann wird der Fahrstuhl als zeitdiskretes System untersucht. Die Untersuchung zeitdiskreter Systeme verlangt bei der Simulation eine völlig andere Vorgehensweise als bei kontinuierlichen Systemen: • Ein kontinuierliches System wird simuliert durch Lösen des dieses System beschreibenden Gleichungsund/ oder Differenzialgleichungssystems. • Ein diskretes System wird simuliert, indem nur die Übergangszeitpunkte betrachtet und registriert werden, bei denen bedingt durch bestimmte Ereignisse Zustandsänderungen der Systemvariablen stattfinden. Dem entsprechend sind diskrete Simulationswerkzeuge prinzipiell anders aufgebaut als die Tools für kontinuierliche Systeme. <?page no="106"?> 92 3.1.2 Konzept der diskreten Simulation Auch zur Simulation diskreter Systeme stehen sehr viele verschiedene Simulationssprachen zur Verfügung. Wenn diese Sprachen auch jeweils auf die unterschiedlichen Belange der Anwendungen zugeschnitten sind, so ist die Software doch immer auf einem ähnlichen Konzept aufgebaut. Sie setzt sich aus vier Komponenten zusammen: • Datenbereiche für die Systemelemente • Simulationsprogramm Das Simulationsprogramm besteht aus verschiedenen Programmstücken, die einzeln aufgerufen werden können und die Zustandsübergänge durch Modifikation der Datenbereiche vornehmen. • Ablaufkontrolle Die Ablaufkontrolle verfügt über eine Liste, die angibt, welches Programmstück zu welcher Zeit aufgerufen werden muss. • Simulationsuhr Die Simulationsuhr gibt die Zeit an. Sie wird in diskreten Schritten weiter geschaltet. Die Ablaufkontrolle überwacht die Simulationsuhr, um entscheiden zu können, welche Zustandsübergänge durchgeführt werden sollen. 3.1.3 Beispiel Fernsprechvermittlung Zur Verdeutlichung der eingeführten Begriffe soll ein Fernsprechnetz gemäß Abbildung 82 zwischen A-Dorf und B-Dorf betrachtet werden, das der Einfachheit halber nur aus drei Leitungen besteht: Abbildung 82: Beispiel Fernsprechvermittlung Verkehrssenken Verkehrsquellen A-Dorf B-Dorf Leitungsbündel (Sprachkanäle) Richtung des Verbindungsaufbaus geschlossener offener Koppelpunkt x = 1 Leitungen belegt Gerufener Teilnehmer Rufender Teilnehmer <?page no="107"?> 93 Allgemein kann man nun nach der Wahrscheinlichkeit dafür fragen, dass ein rufender Teilnehmer in A-Dorf keine freie Leitung mehr erhalten kann und damit den Besetztton hört, weil in diesem Moment bereits x = 3 Ferngespräche von A-Dorf nach B-Dorf geführt werden (Verlustwahrscheinlichkeit). Die im vorigen Kapitel 3.1.2 dargestellten Komponenten der zeitdiskreten Simulation sind in diesem Beispiel: • Systemelemente Leitungen, besetzt/ nicht besetzt • Simulationsprogramm Programmstücke wie: Prozedur ANRUF - Teilnehmer in A-Dorf hat Gesprächswunsch hebt Hörer ab hört Dauerton wählt Nummer für Gespräch nach B-Dorf hört Rufton bei gerufenem Teilnehmer spricht mit gerufenem Teilnehmer und: Prozedur ENDE - Teilnehmer legt Hörer auf • Ablaufkontrolle Definiert „Spielregeln“ Ankunftsprozess: Zwischen 2 Gesprächswünschen in A-Dorf liege die Zeit T a . Bedienungsprozess: Die Dauer eines Gesprächs (1 Leitung besetzt) sei beschrieben durch das Zeitintervall T b . Dabei sind T a und T b Zufallsvariablen, für die jeweils der Mittelwert und die Varianz vorgegeben werden. Betriebsstrategie: Bei blockiertem Leitungsbündel (x=3) erhält ein weiterer rufender Teilnehmer den Besetztton (Verlustsystem). Andere Möglichkeit: Er kann warten, bis eine Leitung frei wird (Wartesystem). Der Ablauf eines typischen Simulationslaufs für das behandelte Beispiel ist in Abbildung 83 auf Seite 94 dargestellt. In den Kästchen oberhalb der Zeitachse sind die Ankunfts-Ereignisse A i eingetragen. Zwischen ihnen liegt jeweils die Zeitspanne T ai . Unterhalb der Zeitachse erkennt man die Ende-Ereignisse E i . Die Gesprächsdauern T bi sind ebenfalls dem Bild zu entnehmen. In der unteren Hälfte des Bildes ist die Anzahl x der momentan im betrachteten Netz geführten Ferngespräche dargestellt. Man sieht, wie sich x im Laufe der Zeit ändert. Zum Zeitpunkt des 5. Anrufwunsches allerdings sind gerade alle 3 Leitungen belegt, und der erste Verlustfall tritt ein. Üblicherweise wird diese Simulation erst nach 10 5 Anrufereignissen abgebrochen, um eine genügend hohe Sicherheit bei der Berechnung der Verlustwahrscheinlichkeit zu erhalten. Sie ist natürlich abhängig von den zufälligen Zeiten T a und T b . Für den begrenzten Stichprobenumfang im Beispiel kann nur die relative Häufigkeit berechnet werden, z.B. die für Abweisen. Diese ergibt sich zu Relative Verlusthäufigkeit = Anzahl der Verlustfälle / Anzahl der Anrufereignisse = 0,2 (1 von 5 Anrufern wurde abgewiesen) <?page no="108"?> 94 Abbildung 83: Typischer Ablauf der Simulation In diesem in die diskrete Simulation einführenden Kapitel sind mehrfach Begriffe aufgetaucht wie Zufallsvariable, Mittelwert, Varianz, Wahrscheinlichkeit usw. Man sieht, dass zur Beschreibung dieser Zusammenhänge die Statistik eine nicht unbedeutende Rolle spielt. Deshalb sollen im nächsten Abschnitt einige Grundlagen vermittelt werden, so weit sie für die Betrachtung eine Rolle spielen. 3.2 Zufallsvariable 3.2.1 Diskrete Zufallsvariable Eine diskrete Zufallsvariable kann innerhalb eines vorgegebenen Zahlenbereichs nur ganzzahlige Werte annehmen. Beispiele hierfür sind: • Augenzahl eines Würfels x = 1,2,3,4,5,6 • Anzahl belegter Leitungen in einem Telefonnetz x = 0,1,2,3 gemäß Abbildung 83 Die statistischen Kenngrößen einer Zufallsvariablen werden durch Messungen ermittelt. Aus diesen Messungen lassen sich Schätzwerte ableiten. A 1 A 2 A 3 A 4 A 5 E 1 E 4 E 2 T b1 T b2 T b3 T b4 T a1 T a2 T a3 Leitungsbündel belegt Verlustfall x 1 2 3 Simulationsuhr <?page no="109"?> 95 Ein Beispiel zeigt Abbildung 84. Aufgetragen über der Zeit ist die Anzahl gleichzeitig belegter Leitungen des Telefonnetzes gemäß Abbildung 82 auf Seite 92. Zu bestimmten Zeitpunkten mit konstantem Abstand werden nun stichprobenartig Messungen der Belegung durchgeführt. Hier sei der Einfachheit halber der Stichprobenumfang u mit nur 10 Messungen klein (üblich sind u = 10 4 ). Man erhält 10 Messungen x 1 bis x 10 . Diese können nun statistisch ausgewertet werden. Abbildung 84: Beispiel einer diskreten Zufallsvariablen 3.2.2 Statistische Kenngrößen diskreter Zufallsvariablen • Absolute Häufigkeit H(x), relative Häufigkeit h(x) Die erfassten Stichproben sind in Tabelle 7 dargestellt: x H(x) h(x) = H(x) / u 0 2 0,2 1 3 0,3 2 1 0,1 3 4 0,4 Tabelle 7: Absolute und relative Häufigkeit Die Summe aller relativen Häufigkeiten ist stets 1,0 . Messergebnisse x 1 2 3 0 1 3 3 3 2 3 1 0 1 X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9 X 10 Messzeitpunkte Stichprobenumfang u = 10 t <?page no="110"?> 96 Auch eine grafische Darstellung in Form eines Histogramms ist möglich, siehe Abbildung 85: Abbildung 85: Histogramm der absoluten Häufigkeiten • Absolute Summenhäufigkeit H(x l), relative Summenhäufigkeit h(x l) Die Summenhäufigkeit gibt an, wie oft bei der Messung genau l oder weniger Leitungen belegt waren. Das Ergebnis in Tabelle 8 zusammengestellt: l H(x l) h(x l) = H(x l) / u 0 2 0,2 1 5 0,5 2 6 0,6 3 10 1,0 Tabelle 8: Absolute und relative Summenhäufigkeit Das Histogramm für diesen Fall zeigt Abbildung 93. Abbildung 86: Histogramm der absoluten Summenhäufigkeit 0 1 2 3 1 2 3 4 x H(x) 5 0 1 2 3 2 4 6 8 l H(x l) 10 <?page no="111"?> 97 • Mittelwert der Stichproben 7 , 1 2 , 1 2 , 0 3 , 0 0 ) x ( h x u ) x ( H x x 3 0 x 3 0 x m = + + + = = = = = • Mittlere quadratische Abweichung der Stichproben (Standardabweichung, Streuung) ( ) 1 1 2 − − = = u x x s u i m i Diese Auswertungen ergeben einen Einblick in die Stichproben. Bei genügend großem Stichprobenumfang lassen sich aus den Messwerten Schätzwerte für die Zufallsvariablen ermitteln, d.h. Werte, die für ∞ ⎯→ ⎯ u beobachtet würden: • Die Wahrscheinlichkeit w(x), dass genau x Leitungen belegt angetroffen werden, ist ) x ( h lim ) x ( w u ∞ ⎯→ ⎯ = . Die Wahrscheinlichkeit ergibt sich aus der relativen Häufigkeit. • Die diskrete Verteilungsfunktion F(l), d.h. die Wahrscheinlichkeit, dass l oder weniger Leitungen belegt angetroffen werden, ist ) l x ( h lim ) l x ( w ) l ( F u ≤ = ≤ = ∞ ⎯→ ⎯ . Die Verteilungsfunktion entsteht aus der relativen Summenhäufigkeit. • Der Erwartungswert x folgt aus dem Mittelwert, wenn man nur die Anzahl der Stichproben genügend erhöht. m 3 0 x u x lim ) x ( w x x = ∞ ⎯→ ⎯ = = <?page no="112"?> 98 3.2.3 Stetige Zufallsvariable Eine stetige Zufallszahl z kann innerhalb eines vorgegebenen Zahlenbereichs unendlich viele verschiedene Werte annehmen. Die Wahrscheinlichkeit für das Auftreten eines bestimmten Wertes z geht daher gegen Null und kann nicht zur Beschreibung der stetigen Zufallsvariablen verwendet werden. Mit der Definition der • Wahrscheinlichkeitsdichte f(z), die angibt, mit welcher Wahrscheinlichkeit z im Bereich zwischen z 1 und z 2 liegt, ≥ 2 1 z z 0 ) z ( f mit , dz ) z ( f lässt sich die stetige Zufallsvariable statistisch erfassen. Außerdem gilt +∞ ∞ − = 1 dz ) z ( f . • Die stetige Verteilungsfunktion F(t) gibt die Wahrscheinlichkeit an, dass z kleiner oder gleich t ist ∞ − = t dz ) z ( f ) t ( F mit f(z) 0. Ein Beispiel der Wahrscheinlichkeitsdichte und der zugehörigen Verteilungsfunktion zeigt Abbildung 87 . Abbildung 87: Wahrscheinlichkeitsdichte und Verteilungsfunktion Zwei weitere Größen seien hier noch angegeben, nämlich der • Erwartungswert z +∞ ∞ − = dz ) z ( f z z <?page no="113"?> 99 und die • Varianz 2 σ +∞ ∞ − − = σ dz ) z ( f ) z z ( 2 2 3.2.4 Beispiel: Messung einer stetigen Zufallsvariablen Bei der Messung einer stetigen Zufallsvariablen werden die Stichproben einzelnen Messbereichen (Klassen) zugeordnet, so z.B. die Gesprächsdauern von Telefongesprächen, siehe Tabelle 9. Klassenhäufigkeit Tb [ sec ] absolut relativ Schätzung für die Verteilungsfunktion 0 - 20 2 0.002 0.002 20 - 40 3 0.003 0.005 40 - 60 36 0.036 0.041 60 - 80 186 0.186 0.227 80 - 100 372 0.372 0.599 100 - 120 296 0.296 0.895 120 - 140 93 0.093 0.988 140 - 160 11 0.011 0.999 160 - 180 1 0.001 1.000 >180 0 0.000 1.000 Tabelle 9: Klassenhäufigkeit von Gesprächsdauern Die Gesprächsdauern werden in Abschnitte von jeweils 20 Sekunden eingeteilt. Damit kann man Klassenhäufigkeiten angeben, wie sie in der zweiten Spalte zu erkennen sind. Neben den absoluten sind auch die relativen Klassenhäufigkeiten von 1000 Stichproben angegeben. In der letzten Spalte wurde daraus eine Schätzung der Verteilungsfunktion dieser Zeiten berechnet. Man erkennt aus der Tabelle, dass sehr wenige Teilnehmer kürzer als 40 Sekunden oder länger als 140 Sekunden telefoniert haben. Der Großteil aller Gespräche dauerte zwischen 1 und 2 Minuten. Dies zeigt sich auch in Abbildung 88 auf Seite 100. Hier wurden die Messergebnisse grafisch dargestellt. Im Bereich zwischen 60 und 120 Sekunden ergibt sich ein steiler Anstieg der Verteilungsfunktion F(t), weil hier auch die meisten Gespräche statt gefunden haben. <?page no="114"?> 100 0 50 100 150 200 250 300 350 400 0 - 20 20 - 40 40 - 60 60 - 80 80 - 100 100 - 120 120 - 140 140 - 160 160 - 180 >180 Absolute Häufigkeit Abbildung 88: Auswertung von Gesprächsdauern 0 100 200 300 400 500 600 700 800 900 1000 0 - 20 20 - 40 40 - 60 60 - 80 80 - 100 100 - 120 120 - 140 140 - 160 160 - 180 >180 Absolute Summenhäufigkeit Absolute Summenhäufigkeit F(t) Absolute Häufigkeit f(t) <?page no="115"?> 101 3.2.5 Typische Verteilungsfunktionen Für jedes Zufallsexperiment, z.B. • Anzahl der Telefongespräche in einem Leitungsbündel, • Anzahl der an einer bestimmten Stelle vorbeifahrenden Autos (dünner Verkehr), • Anzahl der Zerfallserscheinungen eines radioaktiven Präparates, • … gibt es typische, für den entsprechenden Sachverhalt charakteristische Verteilungsfunktionen. Im Laufe der Vergangenheit werteten Statistiker und Mathematiker Experimente aus, die sie mit verschiedenen Systemen der Natur gemacht hatten, und zeichneten bei genügend hohem Stichprobenumfang deren Verteilungsfunktion auf. Danach suchten sie nach einer mathematischen Funktion, die die gemessene Verteilungsfunktion am besten beschrieb. So gibt es unter anderem die • Poisson-Verteilung, • Exponential-Verteilung, • Erlang-k-Verteilung, • Normal-Verteilung (Gauss-Verteilung) • Binomial-Verteilung, • Hypergeometrische Verteilung, • … . Beispielhaft soll hier nur auf einige wenige Verteilungsfunktionen eingegangen werden. Poisson-Verteilung: Die Poisson-Verteilung ist eine diskrete Verteilungsfunktion. Sie beschreibt z.B. die zufällige Anzahl von Ankünften (Telefonate, Kunden, Autos) pro Zeiteinheit. Die Wahrscheinlichkeit, dass in einem Zeitintervall T genau x Ankünfte eintreten, wird hier beschrieben durch ,... 3 , 2 , 1 , 0 x mit ! x e ) T ( ) x ( P T x = λ = λ − Dabei ist λ die Anzahl der im Mittel pro Zeiteinheit eintreffenden Ankünfte (z.B. Kunden). Ist die Anzahl von Ankünften Poisson-verteilt, so ist der zeitliche Abstand dieser Ankünfte negativ exponentiell verteilt. Exponential-Verteilung: Der zeitliche Abstand von Ereignissen stellt eine stetige Zufallsvariable dar. Die Wahrscheinlichkeitsdichte des Ankunftsintervalls wird beschrieben durch z e ) z ( f λ − λ = . <?page no="116"?> 102 Die Verteilungsfunktion F(t) als Integral der Wahrscheinlichkeitsdichte wird somit zu t e 1 ) t ( F λ − − = . Beispielhaft wird hier diese Berechnung dargestellt: t t 0 t 0 z t 0 t 0 z e 1 e e e dz e dz ) z ( f ) t ( F λ − λ − − λ − λ − − = − + = − = = λ = = Die Exponential-Verteilung ist in Abbildung 89 dargestellt für den Fall 1 = λ . 0 0,5 1 0 1 2 3 4 5 Abbildung 89: Exponential-Verteilung Möchte man nun bei der Simulation eines Fernsprechnetzes auf dem Computer eine zufällige Zahlenfolge von Ankunftszeiten realisieren, die exponentiell verteilt ist, so muss man folgendes wissen: Es ist möglich, eine gleichverteilte Zahlenfolge zu berechnen. Hier stehen Unterprogramme zur Verfügung, die rekursiv aufgerufen werden können. Als Ergebnis erhält man eine nicht-zufällige Zahlenfolge, die aussieht, als ob es sich um echte Zufallszahlen handeln würde. Man erhält eine Folge von Zahlen, die nach statistischer Auswertung die gleichen Eigenschaften hat, wie eine gleichverteilte Folge von echten Zufallszahlen. Es handelt sich jedoch um so genannte Pseudo-Zufallszahlen. Für die Simulation ist dies ein Vorteil, weil man bei jedem Simulationsneustart immer dieselbe Folge von „Zufallszahlen“ verwendet. Unterschiedliche Ergebnisse der einzelnen Simulationsläufe ergeben sich also nicht wegen unterschiedlicher Zahlenfolgen, sondern können auf unterschiedliche Modellparameter zurückgeführt werden. Gleichverteilung bedeutet, dass die Zahlenfolge im Bereich von z.B. 0 bis 1 aus Werten zusammengesetzt ist, die gleichhäufig vorkommen. Es ergibt sich damit die in Abbildung 90 auf Seite 103 dargestellte Wahrscheinlichkeitsdichtefunktion: <?page no="117"?> 103 0 0,5 1 0 0,5 1 1,5 z f(z) r=F(t) 0 0,5 1 0 0,5 1 1,5 t Abbildung 90: Gleichverteilung: Wahrscheinlichkeitsdichte und Verteilungsfunktion Zur Berechnung einer exponentiell verteilten Zahlenfolge benötigt man nun deren Umkehrfunktion. Ersetzt man F(t) durch r, so wird aus r = F(t) die neue Funktion t = g(r): t e 1 r λ − − = Logarithmieren der Gleichung führt auf t ) e ln( ) r ln( t λ − = = − λ − 1 und damit auf λ − − = ) r 1 ln( t . Ist nun r eine gleichverteilte Zufallszahl im Bereich 1 r 0 ≤ ≤ , dann ist auch 1-r gleichverteilt. Zur Vereinfachung der Berechnung kann 1-r durch r ersetzt werden. λ = 1 T a ist der Mittelwert der Zufallsvariablen t, d.h., der Mittelwert der Ankunftsabstände. Damit ergibt sich damit folgende Formel zur Ermittlung der Zufallszahl t r log T r log t a − = λ − = . <?page no="118"?> 104 Die Anweisungen z.B. in einem FORTRAN-Programm sehen damit folgendermaßen aus: CALL ZUF( IY, R ) gleichverteilte Zufallszahl T = - TA * ALOG( R ) Ankunftsabstand mit Exponential-Verteilung Erlang-k-Verteilung: Der dänische Mathematiker A.K. Erlang entwickelte für den Fernsprechverkehr die folgende Verteilungsfunktion: − = λ − λ − = 1 k 0 i i t k ! i ) t k ( e 1 ) t ( F Dabei ist λ 1 der Mittelwert von t. Man erkennt, dass für k=1 der Sonderfall der bereits bekannten Exponentialverteilung enthalten ist. Es ergibt sich die in Abbildung 91 dargestellte Verteilungsfunktion für den Fall 1 = λ . 0 0,5 1 0 1 2 3 4 5 t F(t) k=1 k=2 k=5 k= k= Abbildung 91: Erlang-k-Verteilung Ein weiterer Spezialfall ergibt sich, wenn man ∞ → k laufen lässt. Die Verteilungsfunktion ist ebenfalls in Abbildung 91 auf Seite 104 eingetragen. Sie beschreibt getaktete Ereignisse, also nicht-zufällige Ereignisse mit gleichem Abstand voneinander. Für den Fall = 1 ist die Wahrscheinlichkeit null, dass ein Ereignis früher eintritt, als zum Zeitpunkt t = 1. Da es aber auch nie länger dauert, springt die Verteilungsfunktion an dieser Stelle auf den Wert 1. Die Erlang-k-Verteilung stellt somit eine Verallgemeinerung der Exponentialverteilung dar. Systemabhängig ergibt sich durch Messung jeweils ein verschiedener Wert für k. t F(t) <?page no="119"?> 105 Normal-Verteilung (Gauss-Verteilung): Eine Gauss-Verteilung besitzt z.B. die schwankende Bearbeitungszeit in einer Werkzeugmaschine in Abhängigkeit von unterschiedlicher Werkstoffmenge und Materialzähigkeit. Charakteristische Werte für diese Zeiten sind der Mittelwert μ und die Standardabweichung . Normalverteilte Zufallsereignisse haben als Wahrscheinlichkeitsdichtefunktion 2 2 2 ) z ( e 2 1 ) z ( f σ μ − − σ π = (Glockenkurve). Bekanntermaßen wird die Verteilungsfunktion berechnet durch Integration der Wahrscheinlichkeitsdichtefunktion dz e 2 1 dz ) z ( f ) t ( F t 2 ) z ( t 2 2 ∞ − σ μ − − ∞ − σ π = = . Dieses Integral nennt man Gausssches Integral oder Fehlerintegral. Es ist nicht geschlossen lösbar. Für seine Berechnung hilft man sich mit Näherungsverfahren und legt die Funktion in Form einer Wertetabelle ab. Das Gleiche gilt für seine Umkehrfunktion. Auch diese muss zur Berechnung von normalverteilten Zufallszahlen auf dem Rechner abgespeichert werden. Abbildung 92 zeigt den typischen Verlauf der Wahrscheinlichkeitsdichte der Normalverteilung für den Fall = 5. Die Kurve ist symmetrisch zu z = und hat an den Stellen ± Wendepunkte. Abbildung 92: Wahrscheinlichkeitsdichte der Normalverteilung 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 0 1 2 3 4 5 6 7 8 9 10 z f(z) Sigma = 0.5 Sigma = 1.0 Sigma = 2.0 <?page no="120"?> 106 3.3 Betriebsstrategien Wie im Beispiel von Kapitel 3.1.3 dargestellt, spielen im Rahmen der Ablaufkontrolle neben den Ereignistypen und den dazugehörigen Prozessen auch noch die Betriebsstrategien bei der Simulation eine wesentliche Rolle. Bei vollständig belegtem Telefonnetz (3 Leitungen belegt) wurde festgelegt, dass ein weiterer Anrufer den Besetztton hören sollte, also nicht warten durfte, bis wieder eine Leitung frei wurde. Der Anrufer wurde abgewiesen. Es gab keinen Warteplatz für ihn. Ein so genannter Verlustfall trat ein. Je nach Anzahl der Warteplätze w unterscheidet man verschiedene Systeme, siehe Tabelle 10. Anzahl der Warteplätze w System 0 Verlust-System 0 < w < Warte-Verlust-System Warte-System Tabelle 10: Systemunterscheidung als Funktion möglicher Warteplätze Innerhalb einer Warteschlange können unterschiedliche Abfertigungsstrategien realisiert werden: • FIFO first in first out z.B. Tankstelle, Zigarettenautomat • LIFO last in first out z.B. Auswertung des letzten Messwertes • zufällig • Eventuell treten Wartende ohne Bedienung wieder aus der Warteschlange aus (Nachbildung von Ungeduld) Hier sind Regeln für die Simulation festzulegen. Auch bezüglich der Abfertigungsdisziplin von mehreren Warteschlangen untereinander sind verschiedene Möglichkeiten gegeben, siehe Abbildung 93: • Zyklisch • Zufällig • Nach Priorität: z.B. zuerst alle von WS1, dann WS2, dann WS3 Abbildung 93: Abfertigungsstrategien bei Warteschlangen WS1 WS2 WS3 1 Bedienstation <?page no="121"?> 107 3.4 Beispiel Warte-Verlust-System Als Beispiel sei ein System mit einfacher Struktur vorausgesetzt, Abbildung 94. Es soll aus nur einer Wartestation WS und einer Bedienstation BS bestehen und besitzt demnach zwei Systemelemente: Abbildung 94: Beispiel Frisörladen Zum besseren Verständnis kann man sich z.B. vorstellen, dass Kunden einen Frisörladen betreten, um sich die Haare schneiden zu lassen. Es gibt einen Frisör (Bedienstation), der sich jeweils nur einem Kunden widmen kann. Außerdem befindet sich im Laden ein Stuhl, auf dem ein weiterer Kunde Platz nehmen kann um zu warten, bis der Frisör frei wird (Wartestation). Eine mögliche Betriebsstrategie, die von der Ablaufkontrolle verfolgt wird, ist in der folgenden Tabelle 11 zusammengestellt: WS BS Ankommender Kunde frei frei wird sofort bedient frei belegt kann warten belegt frei (Fall tritt nicht auf) belegt belegt wird abgewiesen Tabelle 11: Festlegung der Betriebsstrategie des Beispiels Es handelt sich hier um ein Warte-Verlust-System. Warten ist möglich. Wenn jedoch der eine Warteplatz belegt ist, kann ein ankommender Kunde auch abgewiesen werden. Will man diesen Zusammenhang in ein Programm umsetzen, so wird man für jedes Systemelement eine Speicherzelle vorsehen. Der Inhalt der Zelle legt fest, ob die Station belegt ist oder nicht, wie in Tabelle 12 dargestellt: Inhalt Station 0 frei > 0 belegt Tabelle 12: Kennung der Belegung einer Station des Beispiels WS BS Kunden Stuhl Frisör Frisörladen Eingang Ausgang <?page no="122"?> 108 Weitere Speicherzellen sind zur Beschreibung des Systemverhaltens nötig. Sie sind in der folgenden Tabelle 13 zusammengestellt: Speicher Name Inhalt Bedeutung 1 Zeit Momentane Zeit 2 WS Ta Beim Belegen wird der Ankunftszeitpunkt des Kunden eingetragen 3 BS Tb Beim Belegen wird der Endzeitpunkt des Kunden eingetragen 4 KGES Zähler Alle bisher eingetroffenen Kunden 5 KW Zähler Alle Kunden, die bisher warten mussten 6 KA Zähler Alle Kunden, die bisher abgewiesen wurden Tabelle 13: Darstellung aller Speicherzellen des Beispiels Dabei stellen die Speicher 1 und 4 bis 6 nur Hilfszellen dar, die jedoch für die Auswertung des Systemverhaltens von Bedeutung sind. So sagt, ausreichend großer Stichprobenumfang vorausgesetzt, der Wert von KW, KA und KGES etwas aus über die • Wahrscheinlichkeit für Warten KGES KW W = • Wahrscheinlichkeit für Abweisen KGES KA A = Die Wahrscheinlichkeit B, sofort bedient zu werden, ist also B = 1 - W - A. Zusätzlich wichtig für die Ablaufkontrolle ist die Definition des Ankunftsprozesses und des Bedienprozesses. Anschaulich: Wann kommt der nächste Kunde in den Laden? Wie lange ist der Frisör mit diesem Kunden beschäftigt? Hierfür müssen Ankunftszeiten und Bedienzeiten mit bestimmten statistischen Eigenschaften definiert werden: <?page no="123"?> 109 • Für die Ankunftsabstände t a sei eine Exponentialverteilung vorausgesetzt, mit dem mittleren zeitlichen Abstand 1/ . Dies zeigt Abbildung 95. Abbildung 95: Zeitlicher Abstand der Ankünfte mit Exponentialverteilung • Der Bedienprozess sei beschrieben durch eine einheitliche, d.h. konstante Bediendauer. Sie soll t b = 6 Zeiteinheiten betragen. Die zugehörige Verteilungsfunktion ergibt sich gemäß der folgenden Abbildung 96. Abbildung 96: Verteilungsfunktion für ein getaktetes Ereignis Dabei gilt nt konsta t t b b = = . Dies entspricht einem getakteten Ereignis. Anschaulich: Der Frisör verpasst jedem Kunden seinen Einheitsschnitt. Nach Festlegung aller für die Simulation des Systems nötigen Bedingungen kann ein Simulationslauf erfolgen. Die Ergebnisse für die ersten 15 Zeiteinheiten sind in Abbildung 97 auf Seite 110 dargestellt. Aus Abbildung 98 auf Seite 111 kann man die Struktur der Programme ANK und END entnehmen. Diese Programme steuern den Simulationsablauf. t a 0 0,5 1 0 1 2 3 4 5 t F( t t a ) 0 0,5 1 0 1 2 3 4 5 t F( t t b ) tb __ t b <?page no="124"?> 110 Die Belegung der vorhandenen Speicherzellen über der Zeit kann man aus der folgenden Tabelle 14 erkennen; eingetragen wird nur, wenn sich der Inhalt einer Zelle ändert: Speicher Name Inhalt über der Zeit 1 Zeit 0 1.5 4.5 6.5 7.5 2 WS 0 4.5 0 3 BS 0 7.5 0 ; 13.5 4 KGES 0 1 2 3 5 KW 0 1 6 KA 0 1 Tabelle 14: Simulationsablauf des Beispiels Wertet man den Simulationslauf bereits nach 3 Ankunftsereignissen aus, so ergibt sich die • Wahrscheinlichkeit für Warten 33 . 0 3 1 KGES KW W = = = • Wahrscheinlichkeit für Abweisen 33 . 0 3 1 KGES KA W = = = Üblicherweise erfolgt eine Auswertung erst nach genügend vielen Ereignissen, z.B. nach 10 2 bis 10 3 Kunden des Frisörladens. Es kommt dabei nicht wie im Falle der Lösung einer Differenzialgleichung für ein zeitkontinuierliches System auf den Zeitverlauf an. Es interessiert vielmehr, ob der Durchlauf der Kunden flüssig war, ob vielleicht zu viele Kunden abgewiesen werden mussten und es deshalb besser wäre, weitere Wartestühle aufzustellen oder einen zweiten Frisör einzustellen, … Abbildung 97: Beispiel eines Simulationslaufs t 1.5 4.5 6.5 END 7.5 ANK 6.5 ANK 4.5 ANK 1.5 1.5 t a =3.0 2.0 8.5 END 13.5 ANK 15.0 Beginn der Simulation Ende der Simulation 15.0 7.5 13.5 WS BS Kunde Ankunftsabstände 0.0 <?page no="125"?> 111 Abbildung 98: Programmstruktur Es wird deutlich, dass zur Beschreibung des Ablaufs in den Unterprogrammen einiger Aufwand getrieben werden muss, um alle Änderungen, Ereigniszähler usw. penibel genau zu erfassen. Um es den Anwendern einfacher zu machen, wurden verschiedenste benutzerfreundliche Programmiersprachen entwickelt, die die rein „buchhalterischen“ Aufgaben übernehmen. Ein Beispiel aus der Menge dieser Sprachen wird im Folgenden exemplarisch vorgestellt. Ziel ist, Wirkung und Zusammenspiel der in diesen Sprachen verwendeten Elemente verständlich zu machen. Das Wissen hierüber ist entscheidend bei der Anwendung anderer, neuerer Sprachen der zeitdiskreten Simulation. START BS frei ? BS belegen Belegzeitpunkt: TE = Zelle 1 + TB WS frei ? ja ja nein nein RETURN Zur Kalenderkontrolle Unterprogramm ANK START RETURN nein Unterprogramm END TA Zelle 1 TE Zelle 1 Ankunft zählen KGES + 1 KGES BS freimachen 0 Zelle 3 TE Zelle 3 Wartefall zählen KW + 1 KW Abweisen zählen KA + 1 KA WS belegen TA Zelle 2 TE in Kalender eintragen BS belegen Belegzeitpunkt: TE = Zelle 1 + TB TE Zelle 3 TE in Kalender eintragen WS frei ? ja WS freimachen 0 Zelle 2 <?page no="126"?> 112 3.5 Simulation auf dem Digitalrechner (Beispiel GPSS) 3.5.1 Allgemeines GPSS General Purpose System Simulator (G.Gordon, USA, 1961) • Modellbeschreibung durch Blockdiagramme • 1 Block = 1 Aktivität (engl. activity) z.B. Bedienen • Verbindungslinie zwischen den Blöcken gibt Reihenfolge der Ausführung der Aktivitäten an • 37 Blocktypen, gekennzeichnet jeweils durch Name, Symbol, Datenfeld (Parameter a, b, c, d, usw.) • Transaktionen wandern von einem Block zum nächsten In GPSS werden die Elemente, die durch das Modell laufen, Transaktionen genannt. Beispiele: Kunden in einem Laden, Werkstücke in einer Fabrik • GENERATE-Block: erzeugt Transaktionen • TERMINATE-Block: löscht Transaktionen • 1 Block kann gleichzeitig mehrere Transaktionen aufnehmen Beispiel: Wartespeicher mit mehreren Warteplätzen • Transport einer Transaktion von einem Block zum nächsten geschieht zu einem vorbestimmten Zeitpunkt oder abhängig vom jeweiligen Systemzustand • Blockdiagramm besteht aus Blöcken, die durch eine Kennung markiert sind (Kennung wird vom Übersetzerprogramm automatisch vergeben), Symbolischer Name für die Kennung ist zulässig, wichtig z.B. bei Verzweigungen, siehe Kapitel 3.5.3 3.5.2 Erzeugung und Transport der Transaktionen durch das Modell Während der Simulation wird die „innere Uhr“ (Integer-Variable) fortlaufend inkrementiert, bis das nächste Ereignis stattfindet (Kalendersteuerung). Zu festgelegten Zeiten werden Transaktionen erzeugt im GENERATE-Block: Bildet den zeitlichen Abstand T a zwischen zwei nacheinander in das Modell kommenden Transaktionen nach. Beispiel: Ankunftsabstand zwischen Kunden Symbol: a,b <?page no="127"?> 113 Parameter: a Mittelwert a T b Werteintervall a T Δ Die durch einen Zufallszahlen-Generator berechneten Werte T a sind gleichverteilt im Intervall a a a a a T T ..... T ..... T T Δ + Δ − Ein wichtiger Sonderfall ist dabei: 0 T b a = Δ = . Daraus folgt: . const a T T a a = = = Eingabezeile: GENERATE a, b Beispiel: GENERATE 2, 0 d.h., T a hat den konstanten Wert 2. Bemerkung: • T a kann auch nach einer vorgebbaren Verteilungsfunktion vom Programm ermittelt werden (wird hier nicht behandelt). • Neben den Parametern a und b gibt es weitere, die die Anzahl der Transaktionen begrenzt, ihren frühesten Erscheinungszeitpunkt festlegt, …, (wird hier nicht behandelt). Diese Transktionen bewegen sich im Laufe der Zeit durch das Modell, um es am Ende wieder zu verlassen im TERMINATE-Block: löscht die Transaktionen aus dem Modell Symbol: Parameter: a Anzahl der Transaktionen, die beim Verlassen berücksichtigt werden für den Simulationsende-Zähler: a = 1: Transaktion wird gezählt a = 0: Transaktion wird nicht gezählt Der Fluss der Transaktionen durch das Modell bis zum TERMINATE-Block kann für eine wählbare Zeit gestoppt werden im ADVANCE-Block: Bildet den Zeitbedarf T b für eine Aktivität nach. Beispiel: Bearbeitungsdauer für ein Werkstück in einer Bedienstation a <?page no="128"?> 114 Symbol: Jede Transaktion bleibt im Mittel T b Zeiteinheiten in diesem Block, bis sie zum nächsten Block geht. Parameter: a Mittelwert b T b Werteintervall b T Δ Die auftretenden Werte von T b liegen gleichverteilt im Intervall b b b b b T T ..... T ..... T T Δ + Δ − Ein wichtiger Sonderfall ist dabei: 0 T b b = Δ = . Daraus folgt: . const a T T b b = = = Eingabezeile: ADVANCE a, b Beispiel: ADVANCE 4, 3 D.h., T b nimmt die Werte (4-3) = 1 bis (4+3) = 7 gleichwahrscheinlich an. 3.5.3 Reihenfolge der Ereignisse Die Kalendersteuerung bewirkt, dass die Transaktionen zu bestimmten Zeitpunkten durch das Modell bewegt werden. Bei Gleichzeitigkeit: Die Transaktionen werden nach Prioritätsklassen weitergeleitet (innerhalb jeder Prioritätsklasse: FIFO). In den Blöcken wird keine Verweildauer nachgebildet (Ausnahme: ADVANCE). Daher erfolgt die Bewegung einer Transaktion so lange durch das Blockdiagramm, bis Fall 1: ein ADVANCE-Block mit Parameter a > 0 erreicht wird. Fall 2: ein nächster Block momentan nicht zugänglich ist. Die bewegte Transaktion ist somit blockiert und wird erst automatisch weiterbewegt, wenn die Blockierung beendet ist. Fall 3: ein TERMINATE-Block erreicht ist. Die Transaktion wird gelöscht. Das Programm führt zunächst alle Bewegungen von Transaktionen in einem bestimmten Zeitpunkt durch und schreitet dann zum Zeitpunkt des nächsten Ereignisses in Inkrementalschritten fort. a,b <?page no="129"?> 115 3.5.4 Verzweigungen Um die Transaktionen während der Simulation in verschiedene Richtungen zu lenken, gibt es in GPSS den TRANSFER-Block Symbol: Parameter: a Auswahlmodus b Kennung des nächsten Blockes im Hauptzweig c Kennung des nächsten Blockes im Nebenzweig Es stehen 9 verschiedene Auswahlmodi zur Verfügung. Die Wichtigsten sind: 1. Zufällige Auswahl gemäß vorgegebener Wahrscheinlichkeit r für die Weiterleitung an den Nebenzweig, a = .r Eingabezeile: TRANSFER .r, b, c Beispiel: TRANSFER .100, WEG1, WEG2 100 Promille = 10% der Transaktionen werden zum Block mit der Kennung WEG2 weiter geleitet (Nebenzweig). 90% der Transaktionen wandern zum Block mit der Kennung WEG1 (Hauptzweig). 2. Auswahl gemäß Eintrittsmöglichkeit: a = BOTH Eingabezeile: TRANSFER BOTH, b, c Beispiel: TRANSFER BOTH, WEG1, WEG2 Weiterleitung nach Block b, falls dieser nicht blockiert ist; sonst Weiterleitung nach Block c. Falls Block b und Block c momentan beide blockiert sind: Warten, bis einer der beiden Blöcke b oder c wieder frei sind. a b c Hauptzweig Nebenzweig <?page no="130"?> 116 3.5.5 Einfaches Beispiel Fertigung und Prüfung von Teilen (Modell 1) Gegeben: Eine Werkzeugmaschine stößt nach jeweils T a = 5 Minuten (getaktet) ein Werkstück aus, welches einer Prüfung unterzogen wird. Die Prüfdauer soll unterschiedlich lang sein und T b = 4 ± 3 Minuten betragen (Zeitdauern mit Gleichverteilung). Prüfungsergebnis: Es sollen ca. 10% Ausschuss erzeugt werden. Simulation mit GPSS Für die Simulation muss festgelegt werden, was einer Transktion entspricht. Das sind im Beispiel die Werkstücke. Sie werden in der Maschine erzeugt, laufen durch das Modell und verlassen es nach erfolgter Prüfung wieder. Der Zeittakt für die voranschreitende Zeit ist 1 Minute. Alle Zeiten im Beispiel sind in Minuten angegeben. Also: 1 Werkstück 1 Transaktion 1 Minute 1 Zeiteinheit Aufgabe ist nun, durch geeignete Auswahl an Blöcken den gegebenen Sachverhalt des Textes umzusetzen in ein GPSS-Programm. In diesem Fall ist das zugehörige Blockdiagramm (Modell 1) dargestellt in Abbildung 99: Abbildung 99: GPSS-Blockschaltbild Modell 1 .100 GUT 5,0 4,3 FEHL 1 1 GENERATE ADVANCE TRANSFER TERMINATE Werkzeugmaschine Prüfung (Dauer 1 … 7 Minuten) Prüfungsergebnis (10% fehlerhaft 90% gut) Verlassen des Modells <?page no="131"?> 117 Es ist nicht angegeben, was mit den guten bzw. schlechten Teilen nach der Prüfung geschehen soll. Deshalb sind keine weiteren Blöcke notwendig, etwa für notwendige Nacharbeit, … Die Transaktionen werden deshalb nach erfolgter Prüfung direkt über die TERMINATE-Blöcke gelöscht. Bei jeder in den Blöcken GUT bzw. FEHL gelöschten Transaktionen wird der allgemeine Endezähler um 1 erhöht. Dies geschieht so lange, bis ein Grenzwert erreicht ist, der damit auch das Endekriterium der ganzen Simulation darstellt. Nach Auswahl der Blöcke und der Darstellung im Blockdiagramm erfolgt nun die programmtechnische Umsetzung. Codierung: Sie ist nicht spaltengebunden. Die allgemeine Zeilensyntax ist: Name Block Par1, Par2, Par3,… Kommentar Für Modell 1 folgt damit: * Modell 1 Kommentarzeile beginnt mit * in Spalte 1 * SIMULATE GENERATE 5 statt 5,0 zulässig ADVANCE 4,3 Prüfstation TRANSFER .100, GUT, FEHL Ergebnis GUT TERMINATE 1 FEHL TERMINATE 1 START 1000 Grenzwert des Endezählers END Bemerkungen zum Programm: SIMULATE und bilden den Rahmen für das Programm. END GENERATE Das erste Werkstück wird bei TIME = 5 generiert. START Die Simulation soll beendet werden, wenn 1000 Transkationen das Modell verlassen haben. Im Beispiel spielt es keine Rolle, ob es sich um gute oder schlechte Teile handelt. Deshalb werden die Teile in beiden TERMINATE-Blöcken gezählt. Wollte man simulieren, bis 1000 gute Teile erzeugt sind, dann müsste man den Parameter a des TERMINATE-Blocks auf null setzen, die schlechten Teile also nicht berücksichtigen im Ende- Zähler. <?page no="132"?> 118 Simulationsergebnis: Anzahl der Werkstücke, die im TERMINATE-Block bei GUT bzw. im TERMINATE- Block bei FEHL gezählt wurden. Dies ist abhängig von den Zufallszahlen, die der TRANSFER-Block automatisch und reproduzierbar erzeugt (Pseudozufallszahlen). 3.5.6 Eine oder mehrere gleichartige Bedienstationen (Modelle 2 und 3) Bisher wurde davon ausgegangen, dass jedes Werkstück, welches die Maschine verlässt, auch sofort geprüft werden kann. Wenn wie im letzten Beispiel in Modell 1 gezeigt die Prüfzeit für die Werkstücke im Mittel kürzer (ADVANCE 4) als die Erzeugerzeit (GENERATE 5) ist, wird klar, dass zur Prüfung der Teile ein einzelner Prüfer genügt. Der Prüfvorgang ist bereits beendet, bevor ein neues Teil aus der Maschine kommt. Der Prüfer kann warten auf die nächste Prüfaufgabe. Ist die Erzeugerrate allerdings größer als die Prüfrate, wird sich mit der Zeit ein wachsender Berg von ungeprüften Teilen bilden. Jetzt bedarf es zusätzlicher „Spielregeln“. Es muss festgelegt werden, ob dann der Prüfer nicht jedes Teil, sondern nur stichprobenartig prüfen soll, - oder sich die neuen Werkstücke besser ansammeln sollen, um nicht ungeprüft einen anderen Weg im Fertigungsablauf zu nehmen. Ohne weitere Maßnahme - und so war es im Modell 1 würden die Transaktionen alle direkt in den ADVANCE-Block laufen und dort die Prüfung sofort gestartet werden, als ob beliebig viele Prüfer zur Verfügung ständen. Diese Situation einer begrenzten Anzahl von Prüfern kann mit den bisher verfügbaren Mitteln nicht simuliert werden. Das „Belegen“ nur eines Prüfers mit nur einem Werkstück muss im Modell geeignet nachgebildet werden. Die Belegung oder Freigabe einer oder mehrerer gleichartiger Bedienstationen (hier: Prüfer) geschieht dabei entweder mit den beiden Blocktypen SEIZE und RELEASE oder ENTER und LEAVE, siehe Tabelle 15. Bedeutung Bezeichnung in GPSS Zugehörige Blocktypen 1 Bedienstation FACILITY SEIZE ergreifen RELEASE loslassen Mehrere gleichartige Bedienstationen STORAGE ENTER betreten LEAVE verlassen Tabelle 15: Typen von Bedienstationen Je ein Beispiel für eine FACILITY (Modell 2) und einen STORAGE (Modell 3) zeigen Abbildung 100 auf Seite 119 und Abbildung 101 auf Seite 121. Die Symbole der vier neuen Blöcke können den Abbildungen entnommen werden. Modell 2 und Modell 3 bauen jeweils auf Modell 1 auf; die dazu gekommenen Blöcke für die FACILITY und den STORAGE steuern den Zufluss in den ADVANCE-Block. <?page no="133"?> 119 Über SEIZE wird die Transaktion in die FACILITY zugelassen bzw. über ENTER in den STORAGE (die Prüfumgebung). Läuft eine Transaktion nach Verlassen des ADVANCE-Blocks über den RELEASEbzw. LEAVE-Block, wird der Kalendersteuerung klar, dass ein Prüfplatz in der Prüfstation frei wurde. Eventuell davor Wartende dürfen jetzt nachrücken. 1. Beispiel FACILITY (Modell 2): Abbildung 100: GPSS-Blockschaltbild Modell 2 .100 GUT 5,0 4,3 FEHL 1 1 GENERATE ADVANCE TRANSFER TERMINATE SEIZE RELEASE Prüf Prüf Prüf : Wählbarer Name der Facility <?page no="134"?> 120 Das zugehörige Programm lautet: * Modell 2 PRUEF EQU 1, F Symboldefinition für PRUEF SIMULATE GENERATE 5 SEIZE PRUEF Belegung der FACILITY ADVANCE 4, 3 RELEASE PRUEF Freigabe der FACILITY TRANSFER .100, GUT, FEHL GUT TERMINATE 1 FEHL TERMINATE 1 START 1000 END Bemerkungen zum Programm: EQUIVALENCE-Anweisung: Es gibt 1 Facility im Programm, die automatisch die Nummer 1 erhält. SEIZE 1 bzw. RELEASE 1 wären damit zulässige Anweisungen zur Belegung und Freigabe der FACILITY. Möchte man die FACILITY jedoch mit einem speziellen Namen im Programm ansprechen, so muss dieser frei wählbare Name der Nummer 1 zugeordnet werden. Dies geschieht über die EQUIVALENCE-Anweisung PRUEF EQU 1,F. Der Name PRUEF ist damit der FACILITY 1 zugeordnet. Die beiden Blöcke SEIZE und RELEASE, die eine FACILITY definieren, sorgen im Gegensatz zu Modell1 dafür, dass ein Werkstück im Block SEIZE hängen bleibt, wenn die Prüfstation belegt ist, d.h. im ADVANCE-Block bereits eine Transaktion ist. Erst wenn der Prüfvorgang abgeschlossen ist und die Transaktion den RELEASE- Block durchläuft, wird eine weitere Transaktion aus dem SEIZE-Block zum AD- VANCE-Block vorgelassen. Simulationsergebnis: Ergebnisse wie bei Modell 1; Zusätzlich alle Informationen, die durch Erweiterung des Modells von Bedeutung sind, hier also Ergebnisse für die FACILITY PRUEF: • Gesamtzahl der eingetroffenen Transaktionen in die FACILITY • Auslastung als Prozentsatz der Zeit, in der die FACILITY belegt war • Mittlere Verweildauer einer Transaktion in der FACILITY PRUEF 2. Beispiel STORAGE (Modell 3): Die Produktion von Teilen läuft ab wie bei den Modellen 1 und 2. Zur Prüfung werden aber in Modell 3 drei gleichartige Prüfer und eine geänderte Prüfzeit je Werkstück von 12 ± 5 Minuten. Das zugehörige Blockdiagramm zeigt Abbildung 101 auf Seite 121. <?page no="135"?> 121 Blockdiagramm: Abbildung 101: GPSS-Blockschaltbild Modell 3 Das zugehörige Programm lautet: * Modell 3 PGRUP EQU 1, S Symboldefinition für PGRP SIMULATE GENERATE 5 ENTER PGRUP Belegung eines STORAGE-Platzes ADVANCE 12,5 LEAVE PGRUP Freigabe eines STORAGE-Platzes TRANSFER .100, GUT, FEHL GUT TERMINATE 1 FEHL TERMINATE 1 START 1000 PGRUP STORAGE 3 Größe des STORAGE : 3 END 1 .100 GUT 2,0 12,5 FEHL 1 GENERATE ADVANCE TRANSFER TERMINATE ENTER LEAVE PGRUP PGRUP PGRUP: Name des Storage <?page no="136"?> 122 Bemerkungen zum Programm: EQUIVALENCE-Anweisung: Es gibt 1 STORAGE im Programm, der automatisch die Nummer 1 erhält: ENTER 1 bzw. LEAVE 1 wären damit zulässige Anweisungen zur Belegung und Freigabe des STORAGE. Möchte man den STORAGE jedoch mit einem speziellen Namen im Programm ansprechen, so muss dieser frei wählbare Name der Nummer 1 zugeordnet werden. Dies geschieht über die EQUIVALENCE-Anweisung PGRUP EQU 1,S. Der Name PGRUP ist damit dem STORAGE 1 zugeordnet. Die beiden Blöcke ENTER und LEAVE, die einen STORAGE definieren, sorgen im Gegensatz zu Modell 2 dafür, dass ein Werkstück im Block ENTER hängen bleibt, wenn die Prüfstation belegt ist, d.h. im ADVANCE-Block bereits so viele Transaktionen sind, wie es die Größe des STORAGE zulässt. Erst, wenn der Prüfvorgang einer Transaktion abgeschlossen ist und die Transaktion den LEAVE-Block durchläuft, wird eine Transaktion aus dem ENTER-Block zum ADVANCE-Block vorgelassen. Die Größe des STORAGE mit dem Namen PGRUP ist 3, maximal 3 Transaktionen dürfen sich gleichzeitig im STORAGE befinden. Simulationsergebnis: Ergebnisse wie bei Modell 1 ; Zusätzlich für den Storage PGRP: • Minimale Anzahl gleichzeitig belegter Bedienstationen • Maximale Anzahl gleichzeitig belegter Bedienstationen • Anzahl der momentan belegten Bedienstationen bei Simulationsende • Gesamtzahl der eingetroffenen Transaktionen • Auslastung als Prozentsatz der Zeit, in der der STORAGE belegt war • Mittlere Verweildauer einer Transaktion im STORAGE PGRUP Durch Einführung weiterer Blöcke ist es damit gelungen, die Transaktionen so zu steuern, wie es einem realistischen Ablauf bei einer begrenzten Anzahl von Prüfern entspricht. Aus der Fülle weiterer Blöcke sollen nun einige weitere eingeführt werden, die nicht der Steuerung des Ablaufs dienen, aber zusätzliche informative Messgrößen erzeugen. 3.5.7 Block-Typen für weitere Messungen (Modell 4) • Warteschlangen-Element, falls Blockierung im Modell auftreten kann: - Blocktyp QUEUE a erhöht - Blocktyp DEPART a erniedrigt die Anzahl der Wartenden in Warteschlange a um 1. <?page no="137"?> 123 Zusätzliches Simulationsergebnis: maximale Länge der Warteschlange mittlere Länge der Warteschlange - Verteilungsfunktion der Wartezeiten • Statistik über Durchlaufszeiten der Transaktionen: - Blocktyp MARK registriert Ankunftszeit im MARK-Block - Blocktyp TABULATE registriert Ankunftszeit im TABULATE-Block und wertet die Durchlaufszeit vom MARK-Block bis zum TABULATE-Block aus. Simulationsergebnis: Zur Ermittlung der Klassenhäufigkeit wird der Messwert in eine Tabelle a eingetragen (TABLE) Modell 4, Abbildung 102 auf Seite 124, entspricht Modell 3, nur werden die Werkstücke bei Blockierung eine Warteschlange bilden, und die Daten der Warteschlange werden erfasst. Zudem wird die Zeit zwischen MARK und TABULATE jeder Transaktion gestoppt und statistisch ausgewertet. Damit kommen in Modell 4 alle bisher dargestellten Blöcke zum Einsatz. Das Programm von Modell 4 lautet: * Modell 4 Symboldefinition WART EQU 1, Q für WART PGRUP EQU 1, S für PGRUP TAB EQU 1, T für TAB SIMULATE GENERATE 5,0 Werkstück erzeugen QUEUE WART in Warteschlange einordnen ENTER PGRUP Storage betreten DEPART WART Warteschlange verlassen MARK Zeitpunkt registrieren ADVANCE 12,5 Werkstück prüfen LEAVE PGRUP Storage verlassen TABULATE TAB Durchlaufszeit messen TRANSFER .100, GUT, FEHL Trennung gut/ schlecht GUT TERMINATE 1 Modell verlassen FEHL TERMINATE 1 Modell verlassen START 400 Für 400 Teile simulieren PGRUP STORAGE 3 Storage-Größe TAB TABLE M$1, 5, 5, 7 Tabellen-Format END <?page no="138"?> 124 Blockdiagramm: Abbildung 102: GPSS-Blockschaltbild Modell 4 Bemerkungen zum Programm: Die allgemeine Form zur Programmierung einer Tabelle lautet name TABLE a, b, c, d .100 GUT 5,0 FEHL 1 1 GENERATE ADVANCE TRANSFER TERMINATE ENTER LEAVE PGRUP PGRUP QUEUE DEPART WART WART MARK TABULATE 12,5 TAB <?page no="139"?> 125 Die Parameter sindParameter: a Art des Messwerts; im Beispiel M$1 = die seit dem letzten Durchgang durch einen parameterlosen MARK-Block verflossene Zeit b Obergrenze des unteren Restintervalls c Intervallbreite d Anzahl der Intervalle, einschließlich unterem Restintervall Für das Beispiel gilt gemäß Abbildung 103: Abbildung 103: Klasseneinteilung in Tabelle TABLE 3.5.8 Bedingte Verzweigung (Modell 5) Beispiel: Werkstücke kommen aus der Werkzeugmaschine auf ein Fließband, welches an zwei Prüfern vorbeiläuft. Das Band ist in mehrere Abschnitte aufgeteilt. Die Bandlaufzeit je Abschnitt beträgt 2 Minuten. Der Prüfer, der gerade Zeit hat, nimmt das Werkstück vom Band, um es zu prüfen, Prüfdauer: 8 ± 6 Minuten. Falls beide Prüfer beschäftigt sind, verlässt das Teil die Fertigung ungeprüft. Weitere Aktivitäten bei den geprüften Teilen gibt es nicht. 5 10 15 20 25 30 35 Unteres Restintervall Oberes Restintervall Parameter b c Intervalle 1 2 3 4 5 6 7 d Häufigkeitsklassen der Messgröße M$1 (Parameter a) Werkzeugmaschine Prüfer 1 Prüfer 2 ungeprüft geprüft geprüft <?page no="140"?> 126 Das Blockschaltbild zeigt Abbildung 104: Blockschaltbild: Abbildung 104: GPSS-Blockschaltbild Modell 5 GENERATE 5,0 1 ADVANCE TERMINATE MARK TABULATE TAB BOTH ADVANCE SEIZE RELEASE PRF1 PRF1 8,6 2,0 BOTH TRANSFER ADVANCE SEIZE RELEASE PRF2 PRF2 2,0 8,6 MARKE WEITR VERL Sprung zu MARKE 0 TRANSFER TRANSFER <?page no="141"?> 127 Das Programm von Modell 5 lautet: * Modell 5 PRF1 EQU 1, F Symboldefinition für PRF1 PRF2 EQU 2, F Symboldefinition für PRF2 TAB EQU 1, T Symboldefinition für TAB * SIMULATE GENERATE 5 MARK ADVANCE 2 * TRANSFER BOTH, , WEITR SEIZE PRF1 Belegung der Facility PRF1 ADVANCE 8, 6 RELEASE PRF1 Freigabe der Facility PRF1 MARKE TABULATE TAB TERMINATE 1 * WEITR ADVANCE 2 TRANSFER BOTH, , VERL SEIZE PRF2 Belegung der Facility PRF2 ADVANCE 8, 6 RELEASE PRF2 Freigabe der Facility PRF2 TRANSFER , MARKE unbedingter Sprung * VERL TERMINATE 0 * TAB TABLE M$1,5,5,10 * START 1000 * END Bemerkungen zum Programm: ADVANCE 2: Mit dieser Vereinbarung wird erreicht, dass die Werkstücke erst 2 Zeittakte später an der nächsten Station ankommen, im Beispiel den TRANSFER- Blöcken. SEIZE PRF1: folgt unmittelbar im Hauptzweig des TRANSFER-Blocks, ohne besondere Sprungmarke. Eine Angabe des Parameters b im TRANSFER-Block ist in diesem Fall nicht nötig. TRANSFER , MARKE: Dies ist ein unbedingter Sprung nach MARKE. Eine Bedingung vor dem Komma muss in diesem Fall nicht angegeben werden. Alle Werkstücke, die in diesem Zeig des Programms laufen, sollen an der Stelle MARKE weiterlaufen. <?page no="142"?> 128 Ein letzter Block soll noch erwähnt werden, mit dem man ähnlich dem TRANSFER- Block den Fluss der Transaktionen steuern kann, der Test-Block Die Umlenkung der Transaktionen in verschiedene Programmzweige erfolgt durch logischen Vergleich zweier Zahlenwerte, den so genannten Standard Numeric Attributes SNA, die den Systemzustand charakterisieren. Ist die Bedingung erfüllt, läuft die Transaktion weiter in die nächste Programmzeile, in Richtung Hauptzweig. Ist sie nicht erfüllt, springt die Transaktion an eine durch eine Sprungmarke gekennzeichnete andere Stelle im Programm, in Richtung Nebenzweig. Symbol: Parameter: aux Vergleichsoperator; Der Vergleich kann wahr oder nicht wahr sein. E Equal = NE Not Equal L Less < LE Less Equal G Greater > GE Greater Equal a 1. Standard Numeric Attribute b 2. Standard Numeric Attribute c Kennung des nächsten Blockes, in den die Transaktion läuft, wenn die Bedingung nicht erfüllt ist (Nebenzweig) Eingabezeile: TEST_aux SNA1, SNA2, WEG Beispiel: TEST_GE Q$1, 10, WEG3 Im Folgenden sollen die Standard Numeric Attributes SNA näher erläutert werden. nicht wahr Vergleich a mit b c Hauptzweig Nebenzweig wahr <?page no="143"?> 129 SNA Standard Numeric Attributes: Standard Numeric Attributes sind Zahlenwerte, die den mit der Zeit veränderlichen Systemzustand beschreiben. Diese Werte werden vom Programm automatisch aktualisiert. Damit dienen sie auch als Entscheidungshilfe, um z.B. Transaktionen zustandsabhängig durch das Programm zu leiten. Beispiele für SNA-Werte: Zeit AC$1 absolute Zeit C$1 relative Zeit Block-Count von Block 1 N$1 total entries W$1 momentaner Inhalt Zufallszahlengenerator 1 RN$1 Random Number (0 - 999) STORAGE 3 S$3 momentaner Inhalt Average SA$3 mittlerer Inhalt Maximum SM$3 maximaler Inhalt QUEUE 2 Q$2 momentaner Inhalt Average QA$2 mittlerer Inhalt Maximum QM$2 maximaler Inhalt Zeitbedarf bzw. M$1 Zeitdauer seit GENERATE bzw. „Alter“ der Transaktion bzw. seit dem letzten MARK-Block Anwendungsbeispiele im Block TEST: 1. TEST_GE RN$1, 100, WEG2 Erklärung: Wenn die Zufallszahl von Zufallszahlengenerator 1 kleiner ist als 100, dann leite die Transaktion nach WEG2. Bei der erzeugten Gleichverteilung der Zufallszahlen im Bereich von 0 bis 999 bedeutet dies, dass 10% der Transaktionen nach WEG2 laufen. 2. TEST_LE Q$2, 10, WEG1 Erklärung: Wenn die Länge der Warteschlange 2 größer als 10 ist, dann leitet der Block die Transaktion nach WEG1 (und eben nicht in die nächste Zeile des Programms, wo z.B. die Transaktionen in Richtung zur Warteschlange 2 laufen). Hiermit könnte nachgebildet werden, dass sich Kunden nur dann an einer Schlange anstellen, wenn diese nicht zu lang ist. <?page no="144"?> 130 3.5.9 Simulation Modell 4 Als Beispiel soll Modell 4 aus Kapitel 3.5.7 simuliert werden. Allerdings sollen nur 2 Prüfer eingesetzt werden, der STORAGE hat also die Größe 2. Nach Aufruf des Programms Modell4.gps durch das Kommando GPSS Modell4.gps wird die Simulation durchgeführt und die Ergebnisse in einer Datei MODELL4.LST aufgelistet, die hier in Ausschnitten wiedergegeben wird: 1. Zuerst wird das Programm dargestellt, inklusive der von GPSS den einzelnen Blöcken zugeteilten Nummern: Line Block 1 * MODELL 4 2 WART EQU 1, Q 3 PGRUP EQU 1, S 4 TAB TABLE 1, T 5 * 6 SIMULATE 7 1 GENERATE 5 8 2 QUEUE WART 9 3 ENTER PGRUP 10 4 DEPART WART 11 5 MARK 12 6 ADVANCE 12,5 13 7 LEAVE PGRUP 14 8 TABULATE TAB 15 9 TRANSFER .100, GUT, FEHL 16 10 GUT TERMINATE 1 17 11 FEHL TERMINATE 1 18 * 19 PGRUP STORAGE 2 20 TAB TABLE M$1,5,5,19 21 * 22 START 400 23 * 24 END • Die Simulation lief also so lange, bis das 400. Teil geprüft war. Dabei kam es nicht darauf an, wie das Ergebnis ausfiel, die Werkstücke beider TERMINATE-Blöcke 10 und 11 wurden berücksichtigt. • Die Prüfzeit war für jedes Teil unterschiedlich und konnte (gleichverteilt) Werte zwischen 7 und 17 Zeiteinheiten annehmen. 2. Als erstes Ergebnis wird die Endezeit angegeben: ABSOLUTE CLOCK 2419 <?page no="145"?> 131 Interpretation des Ergebnisses: • Nach 2419 Zeitschritten war die Simulation beendet, das 400. Teil war durch einen der beiden TERMINAL-Blöcke gelaufen. Die Größenordnung des Ende-Zeitpunktes erscheint vernünftig: 400 Teile * 12 Zeiteinheiten im Mittel zur Prüfung / 2 Prüfer = 2400 • Dazu kommen noch die Zeiten, bis ein erstes Teil erzeugt war. Das erste Teil wird zum Zeitpunkt 5 erzeugt. Und die Prüfzeit streut außerdem um ±5 Zeiteinheiten. Es dauert immer eine gewisse Zeit, bis sich nach Start der Simulation der Ablauf mit dem Anfangszustand „alle Blöcke sind leer“ eingespielt hat. 3. Danach folgt eine Auflistung der so genannten BLOCK COUNTS. Es wird die Anzahl der in den Blöcken 1 bis 11 eingetroffenen Transaktionen, hier Werkstücken, dargestellt (TOTAL: insgesamt; CURRENT: die bei Simulationsende im Block anwesenden): BLOCK COUNTS BLOCK CURRENT TOTAL 1 1 484 2 82 483 3 0 401 4 0 401 5 0 401 6 1 401 7 0 400 8 0 400 9 0 400 10 0 362 11 0 38 Interpretation des Ergebnisses: • Die Maschine hat bei Simulationsende 483 Werkstücke gefertigt, das 484. Teil ist noch in Arbeit (GENERATE-Block 1) • 483 Teile sind in eine Warteschlange gelaufen (QUEUE-Block 2); 401 Teile davon haben die Warteschlange verlassen (DEPART-Block 4). Die Warteschlange ist am Ende 82 Teile lang. • Diese 401 Teile kamen in der Prüfstation (STORAGE) an, im ENTER- Block 3, und wurden von den beiden Prüfern nacheinander geprüft (AD- VANCE-Block 6). • Nach 399 geprüften Teilen waren kurz vor Simulationsende die Teile 400 und 401 bei den Prüfern in Bearbeitung. Eines der beiden war zuerst fertig geprüft, verließ den STORAGE (LEAVE-Block 7) und lief als 400. Element durch einen der beiden TERMINATE-Blöcke 10 oder 11. Das andere Teil verblieb im STORAGE, dort im ADVANCE-Block 6. • Es sollen 10% Ausschuss erzeugt werden. Dies wird über den Transfer- Block 9 mit 38 Werkstücken von 400 auch nahezu erreicht. <?page no="146"?> 132 4. Im vierten Teil der Ergebnisliste stehen die Daten des STORAGE PGRUP: STORAGE PGRUP CAPACITY 2 AVERAGE CONTENT 1.99 AVERAGE UTILIZATION 1.00 ENTRIES 401 AVERAGE TIME / TR 12.01 CURRENT CONTENT 1 MAXIMUM CONTENT 2 Interpretation des Ergebnisses: • Im Mittel befanden sich 1.99 Teile in der Station. Nicht genau 2, denn zu Anfang mussten die Prüfer auf die ersten Teile warten und hatten noch nichts zu tun. • Im Mittel war die Aufenthaltsdauer der Teile etwas größer als die vorgesehene Dauer von 12 Zeiteinheiten. Dies liegt am Zufallszahlengenerator für die Prüfzeiten: 400 Zufallszeiten sind zu wenig, um den Mittelwert genauer zu erreichen. • Maximal 2 der insgesamt vorhandenen 2 Prüfer waren gleichzeitig beschäftigt. Dies ist nicht verwunderlich; die Prüfer müssen jedes Werkstück lange prüfen, und neue Teile kommen schneller dazu, alle 5 Zeiteinheiten. • Es sind einfach zu wenige Prüfer. • Deshalb baut sich auch sukzessive eine Warteschlange auf. 5. Im fünften Teil der Ergebnisliste sind Daten der Warteschlange QUEUE zusammengestellt: QUEUE WART´ MAXIMUM CONTENT 82 AVERAGE CONTENT 41.63 TOTAL ENTRIES 483 ZERO ENTRIES 6 PERCENT ZERO 1.24 AVERAGE TIME/ TR 208.51 $AVERAGE TIME/ TR 211.13 CURRENT CONTENT 82 Interpretation des Ergebnisses: • Die Zahlenwerte deuten auf ein stetiges Wachstum der Warteschlange hin. Zu Anfang war die Länge Null, am Ende mit 82 Teilen am längsten, der Mittelwert beträgt 41. • 6 ZERO ENTRIES bedeutet, dass nur 6 mal der Fall eintrat, dass ein Teil aus der Maschine kam und eine leere Warteschlange vorfand. • Entsprechend gibt es zwei verschiedene mittlere Aufenthaltsdauern, einmal mit Berücksichtigung nur der Teile, die wirklich warten mussten (477), zum anderen mit Berücksichtigung aller durchgelaufenen Teile (483). <?page no="147"?> 133 6. Am Ende wird die Tabelle TAB ausgedruckt mit der Häufigkeitsverteilung entsprechend der programmierten Klasseneinteilung für die benötigten Durchlaufzeit jeder Transaktionen vom MARK-Block bis zum TABULATE-Block. Alle den Benutzer eventuell interessierenden statistischen Daten werden ausgegeben, siehe Tabelle 16. TABLE NO. 1 TAB ENTRIES IN TABLE MEAN ARGUMENT STANDARD DEVIATION SUM OF ARGUMENTS 400 12.02 3.07 4810.0 UPPER LIMIT FRE- QUENCY PERCENT OFTOTAL CUMMU- LATIVE PERCENTAGE CUMMU- LATIVE REMAINDER MULTIPLE OF MEAN DEVIATION FROM MEAN 5 0 0.00 0.00 100.00 0.42 -2.29% 10 136 34.00 34.00 66.00 0.83 -0.66% 15 194 48.50 82.50 17.50 1.25 0.97% 20 70 17.50 100.00 0.00 1.66 2.60% REMAINING VALUES ARE ZERO Tabelle 16: Klassenhäufigkeit der Durchlaufzeiten in Modell 4 Interpretation des Ergebnisses: • Die meiste Zeit verstreicht während der Prüfung im entsprechenden ADVANCE-Block 6 (ADVANCE 12, 5). Die Aufenthaltsdauer im STORAGE beträgt zwischen 7 und 17 Zeiteinheiten. • Die Klasseneinteilung ist also nicht besonders gut gewählt: Beginn bei 5 und Klassenbreite 5. Die Teile fallen in die Klassen (5 - 10), (10 - 15) und (15 - 20). Eine Gleichverteilung ist so nicht zu erkennen. Und die Klassen ab Klasse 20 - 25 bleiben leer und werden nicht in die Tabelle aufgenommen (Remaining Values are Zero). • Für eine bessere Klasseneinteilung sollten deshalb bei weiteren Simulationsläufen 10 Klassen ab 7 mit einer Breite von 1 gewählt werden. Im Gegensatz zu den Simulationsbeispielen der zeitkontinuierlichen Simulation, wo hauptsächlich der Verlauf der Größen über der Zeit interessiert, sind bei der zeitdiskreten Simulation andere Größen von Bedeutung. Man möchte wissen, ob der Fluss der Werkstücke (Transaktionen) durch das Modell problemlos verläuft, ob es Engpässe gibt, Stauraum für Teile benötigt wird oder ob Mitarbeiter überlastet sind. Am Produktionsablauf des Beispiels erkennt man, dass die gewählten Parameter nicht ideal sind: • Die Maschine produziert alle 5 Zeiteinheiten ein Werkstück. • 2 Prüfer müssen jedes Werkstück aufwendig (im Mittel 12 Zeiteinheiten) prüfen und sind damit überfordert (Stress; Prüfqualität fraglich). • Deshalb bildet sich eine stetig wachsende Warteschlange. Damit einher geht ein stetig wachsender Bedarf an Stauraum (Logistik, evtl. Gebäudekosten). <?page no="148"?> 134 In weiteren Simulationsläufen wäre nun zu prüfen, welche Maßnahmen am besten wirken, um diese Probleme zu beseitigen. Lösungsansätze könnten z.B. sein: • Maschine etwas langsamer laufen lassen (GENERATE 6) • Einsatz von weiteren Prüfern (Vergrößerung des STORAGE auf 3 oder 4) • Reduktion des Prüfaufwands (ADVANCE 10,7) • Reduktion des Prüfaufwands durch Beschränkung auf stichprobenartige Prüfung der Teile. (Dann müsste ein weiterer TRANSFER-Block dafür sorgen, dass nur ein bestimmter Prozentsatz der Teile im STORAGE ankommt, die anderen Teile müssten ungeprüft an der Prüfstation vorbei laufen.) Dadurch könnten die Systemparameter besser aufeinander abgestimmt werden und die Staubildung wäre beseitigt. Mit diesem Beispiel zur zeitdiskreten Simulation ist auch das letzte Kapitel über die Simulationstechnik abgeschlossen. In den Anhängen sind ein paar weitere Informationen zum bisherigen Inhalt dargestellt. Für den interessierten Leser gibt es zum Abschluss Angaben zu weiterführender Literatur. <?page no="149"?> 135 4 Anhang: ACSL 4.1 ACSL Model Definition Statements (Examples) Statement Explanation ABS (x) Absolute value ACOS (x) Arc-cosine; result in radians AINT (x) Integer part of real expression x ALOG (x) Natural logarithm ARRAY v(1,2,3) Specify up to three dimensions BOUND (bb,tb,x) Limit expression x to be between bottom and top bounds CINTERVAL CINT=0.1 Define name and value for communication interval COS (x) Cosine of expression x in radians CONSTANT d=a Set constant d to value a DEAD (bb,tb,x) Dead zone between bottom and top bounds DELAY (x,ic,tdl, nmx) Delay expression through fixed time tdl DERIVT (ic,x) Differentiate expression x DYNAMIC Begins section entered every communication interval EXP (x) Natural exponent of expression x INITIAL Begins section executed at beginning of every run (i.e., after START) INTEG (xd, xic) Integrates derivative xd; given initial value xic LIMINT (yd,ic,bb,tb) Limited integrator LSW (p,j1,j2) Integer switch; returns j1 when p is TRUE, else j2 OUTPUT (v1,v2,…) Record values of v i each communication interval PROGRAM string First card of model definition deck; dollar sign not allowed in string SQRT (x) Square root of expression x; error if x negative STEP (tz) Step function; output changes from zero to one at time tz TERMT (lexpr) Identifies run termination condition; forces transfer to TERMINAL-section when logical expression lexpr is TRUE TERMINAL Begins section entered at termination of a run VARIABLE T [TIC=0.0] Defines name and initial condition of independent variable Tabelle 17: ACSL Model Definition Statements (Examples) <?page no="150"?> 136 4.2 ACSL Runtime Executive Commands (Examples) Command Subcommands Explanation DISPLY Display values of named variables END Complete PROCED definition OUTPUT ‘NCIOUT’ = ‘CLEAR’ List variables to be output to DIS unit during run Number of communication intervals between OUTPUT Clear OUTPUT list PLOT ‘ALL’ ‘CHAR’ = ‘HI’ = ‘LO’ = ‘TAG’ = ‘LOG’ ‘OVER’ ‘XAXIS’= ‘XHI’ = ‘XLO’ = ‘XTAG’ = ‘XLOG’ Printer and/ or line and/ or strip plots Plot all variables on PREPAR list Use given character for plot symbol Specify upper Y axis value Specify lower Y axis value Specify Y axis character tag string Use logarithmic scale for Y axis Overplot; no axis drawn Specify X axis variable Specify right-most X axis value Specify left-most X axis value Specify X axis character tag string Use logarithmic scale for X axis PREPAR ‘CLEAR’ List variables to be saved on intermediate file during run Clear prepared list PRINT ‘NCIPRN’ = ‘ALL’ Print variables listed in column form; variables must have been listed previously in PREPAR statement Number of communication intervals between print lines Print all variables on PREPAR list PROCED Define beginning of runtime PRECEDure blocks; Block terminated by END RANGE ‘ALL’ ‘IHI’ = ‘ILO’ = ‘IVAR’ = Print maximum and minimum values of named variables; Variables must be on PREPAR list Determine range of all variables on PREPAR list High value for independent variable Low value for independent variable Define independent variable SET Set values of model constants START Start simulation run STOP Terminate ACSL runtime session; return to system XERROR Establish absolute errors for named state variables Tabelle 18: ACSL Runtime Executive Commands (Examples) <?page no="151"?> 137 4.3 ACSL Simulationen mit dem nichtlinearen Pendel IM FOLGENDEN WIRD EIN GEDAEMPFTES, NICHTLINEARES PENDEL UNTERSUCHT. DAS ACSL-PROGRAMM STEHT IM FILE PENDLM.CSL ZUR VERFUEGUNG. DAS PENDEL WIRD DURCH DIE REIBUNG AM DREHPUNKT GEDAEMPFT. BEI ENTSPRECHENDER STARTGESCHWINDIGKEIT KANN SICH DAS PENDEL BELIEBIG OFT UEBERSCHLAGEN. ES SOLLEN KLEINE SCHWINGUNGEN, GROSSE SCHWINGUNGEN UND UEBERSCHLAEGE BETRACHTET WERDEN. ___________________________________________________________________________ TITLE VEREINBART EINE UEBERSCHRIFT FUER ALLE ZEICHNUNGEN S TITLE='VERHALTEN DES NICHT-LINEAREN PENDELS' MIT S (FUER SETZEN) UND VARIABLENNAMEN (Z.B. TITLE) KOENNEN DIESEN VARIABLEN NEUE WERTE ZUGEWIESEN WERDEN. ___________________________________________________________________________ VOR DEM ERSTEN LAUF MUESSEN MIT PREPAR ALLE VARIABLEN VEREINBART WERDEN, WELCHE VON INTERESSE SIND. DIE WERTE DER VARIABLEN WERDEN DANN BEI ALLEN FOLGENDEN LAEUFEN AUF DAS PREPAR-FILE GESCHRIEBEN. DAS PREPAR-FILE HAT DIE ENDUNG .RRR . VON INTERESSE SIND: T ZEIT (SEKUNDEN) TH WINKEL THETA (BOGENMASS) THD WINKELGESCHWINDIGKEIT (BOGENMASS/ SEKUNDE) THDD WINKELBESCHLEUNIGUNG (BOGENMASS/ SEKUNDE/ SEKUNDE) PREPAR T, TH, THD, THDD ___________________________________________________________________________ MIT OUTPUT WERDEN NUN VARIABLEN VEREINBART, WELCHE WAEHREND ALLEN SIMULATIONSLAEUFEN AM BILDSCHIRM ERSCHEINEN. ES KOENNEN NUR VARIABLEN DER PREPAR-LISTE AUSGEGEBEN WERDEN MIT NCIOUT=X WIRD VEREINBART, DASS NUR JEDER X-TE WERT, WELCHER INS PREPAR-FILE GESCHRIEBEN WIRD, AUCH AM BILDSCHIRM ERSCHEINT. OUTPUT T, TH, THD, 'NCIOUT'=20 WENN ERFORDERLICH KOENNEN DIE PREPAR- UND OUTPUT-VEREINBARUNGEN SPAETER GEAENDERT WERDEN. ___________________________________________________________________________ BEZEICHNUNGEN DER STARTBEDINGUNGEN: THIC STARTWINKEL THDIC STARTWINKELGESCHWINDIGKEIT MIT D KOENNEN AUGENBLICKLICHE WERTE VON VARIABLEN AM BILDSCHIRM ANGEZEIGT WERDEN. SET THIC=1.0 SET THDIC=0.0 D THIC, THDIC THIC 1.00000000 THDIC 0. <?page no="152"?> 138 MIT START BEGINNT NUN DIE EIGENTLICHE SIMULATION MIT DEN VORGEGEBENEN BEDINGUNGEN. WAEHREND DER SIMULATION WERDEN DIE WERTE DER OUTPUT-VARIABLEN AUF DEN BILDSCHIRM UND DER PREPAR-VARIABLEN AUF DAS PREPAR-FILE GESCHRIEBEN. START T 0.0 TH 1.00000000 THD 0.0 T 0.50000000 TH-0.27461500 THD-3.47618000 T 1.00000000 TH-0.74189600 THD 2.02740000 T 1.50000000 TH 0.70504800 THD 1.79906000 T 2.00000000 TH 0.21992100 THD-3.08541000 T 2.50000000 TH-0.77442400 THD 0.33405600 T 3.00000000 TH 0.32556300 THD 2.55948000 T 3.50000000 TH 0.48188900 THD-2.04242000 T 4.00000000 TH-0.62459100 THD-0.87337500 T 4.50000000 TH 0.00619890 THD 2.51286000 T 5.00000000 TH 0.56011700 THD-0.96356500 __________________________________________________________________________ MIT TITLE(11)= WIRD EINE WEITERE UEBERSCHRIFT FUER DIE ZEICHNUNGEN VEREINBART. S TITLE(11)='STARTWINKEL=1.0' ___________________________________________________________________________ MIT PRINT KOENNEN TABELLEN DER PREPAR-VARIABLEN AUF DAS FILE PENDEL.OUT GESCHRIEBEN WERDEN. MIT NCIOUT WIRD VEREINBART, NUR JEDEN 25-STEN WERT DER PREPAR-LISTE ZU SCHREIBEN. MIT ALL WIRD VEREINBART, DASS ALLE VARIABLEN DER PREPAR-LISTE GESCHRIEBEN WERDEN. AUF DAS FILE PENDEL.OUT KANN ERST NACH DEM ENDE DIESER RUNTIME-SITZUNG ZUGEGRIFFEN WERDEN. MIT HVDPRN= (T FUER TRUE, F FUER FALSE) WIRD VEREINBART, OB DIE TABELLE AUCH AUF DEM BILDSCHIRM ERSCHEINEN SOLL ODER NUR INS FILE PENDLM.OUT GESCHRIEBEN WIRD. S HVDPRN=.T. $ PRINT AUCH AUF BILDSCHIRM PRINT 'NCIPRN'=25, 'ALL' LINE T TH THD THDD 0 0.0 1.0000000 0.0 -13.547700 25 0.6250000 -0.6543550 -2.4878700 10.279900 50 1.2500000 -0.0014729 3.4217200 -0.6373630 75 1.8750000 0.5679320 -2.3699800 -8.2021400 100 2.5000000 -0.7744240 0.3340560 11.194300 LINE T TH THD THDD 125 3.1250000 0.5905300 1.6026500 -9.2741300 150 3.7500000 -0.1538310 -2.6134400 2.9718300 175 4.3750000 -0.2995140 2.2739000 4.3110800 200 5.0000000 0.5601170 -0.9635650 -8.3675300 S HVDPRN=.F. $ 'PRINT NICHT AUF BILDSCHIRM' ___________________________________________________________________________ MIT RANGE WERDEN DIE EXTREMWERTE (MINIMUM UND MAXIMUM) DER VARIABLEN AUF DEN BILDSCHIRM GESCHRIEBEN. MIT TIMES WIRD ANGEZEIGT, ZU WELCHEM ZEITPUNKT DER EXTREMWERT AUFTRITT <?page no="153"?> 139 RANGE 'TIMES', 'ALL' T 0. AT 0.0 5.00000000 AT 5.000000 TH-0.91880100 AT 0.825000 1.00000000 AT 0.0 THD-3.69781000 AT 0.400000 3.42172000 AT 1.250000 THDD-13.5477000 AT 0.0 12.8163000 AT 0.800000 SOLLEN DIE EXTREMWERTE NUR IN EINEM BESTIMMTEN BEREICH EINER BESTIMMTEN VARIABLEN GESUCHT WERDEN, SO MUSS DIE VARIABLE MIT IVAR UND DER BEREICH MIT ILO UND IHI VEREINBART WERDEN RANGE 'IVAR'=T, 'ILO'=2., 'IHI'=3., 'TIMES', 'ALL' T 2.02500000 AT 2.025000 3.00000000 AT 3.000000 TH-0.77926500 AT 2.475000 0.32556300 AT 3.000000 THD-3.16878000 AT 2.050000 2.93478000 AT 2.875000 THDD-5.64395000 AT 3.000000 11.3334000 AT 2.450000 ___________________________________________________________________________ MIT PLOT KOENNEN NUN ERGEBNISZEICHNUNGEN AM BILDSCHIRM ERZEUGT WERDEN. ZUVOR SOLLTE ABER DIE ART DER ZEICHNUNG VEREINBART WERDEN. MIT CALPLT = WIRD VEREINBART OB LINIENPLOTS GEZEICHNET WERDEN MIT GRDPLT= WIRD VEREINBART OB GITTERLINIEN GEZEICHNET WERDEN S CALPLT=.T. $ ZEICHNE LINIENPLOTS (NUR EIN ACHSENKREUZ) S GRDCPL=.T. $ ZEICHNE GITTERLINIEN ___________________________________________________________________________ DAS ERSTE PLOT-KOMMANDO PLOT 'XTAG'='(SEC)',TH, THD ' BEDEUTET: ZEICHNE DEN WINKEL TH UND DIE WINKELGESCHWINDIGKEIT THD. UEBER DIE ERSTE VARIABLE DER PREPAR-LISTE (UEBER T). BESCHRIFTE X-ACHSE MIT (SEC) ('XTAG'='(SEC)'). DIE SKALIERUNG WIRD VON ACSL UEBERNOMMEN. DIE UEBERSCHRIFT WURDE SCHON (MIT TITLE=) VEREINBART. PLOT 'XTAG'='(SEC)', TH, THD, $ PAUSE <?page no="154"?> 140 DAS ZWEITE PLOT-KOMMANDO PLOT 'XAXIS'=TH, 'XTAG'='(RAD)', THD, 'TAG'='(RAD/ SEC)' BEDEUTET: NEUE X-ACHSE=TH ('XAXIS'=TH). BESCHRIFTE X-ACHSE MIT (RAD). ZEICHNE THD UEBER X-ACHSE (TH). BESCHRIFTE Y-ACHSE MIT (RAD/ SEC) ( TAG'='(RAD/ SEC'). PLOT 'XAXIS'=TH, 'XTAG'='(RAD)', THD, 'TAG'='(RAD/ SEC)' $ PAUSE ___________________________________________________________________________ MIT STRPLT= WIRD VEREINBART, OB STRIP-PLOTS GEZEICHNET WERDEN STRIP-PLOTS ZEICHNEN DIE VARIABLEN NICHT IN EIN DIAGRAMM,' SONDERN IN MEHRERE DIAGRAMME UEBEREINANDER. MIT XINSPL= WIRD VEREINBART WIE LANG DIE X-ACHSE DER STRIP-PLOTS WIRD (IN INCH, 1 INCH = 2.54 CM ). S CALPLT=.F. $ ZEICHNE KEINE LINIENPLOTS S STRPLT=.T. $ ZEICHNE STRIP-PLOTS S XINSPL=5.0 $ STRIP-PLOT-X-ACHSE= 5 INCH LANG <?page no="155"?> 141 DAS DRITTE PLOT-KOMMANDO PLOT 'XAXIS'=T, 'XTAG'='(SEC)', TH, 'TYPE'=6, THD, 'TYPE'=2 BEDEUTET: NEUE X-ACHSE=T. BESCHRIFTE X-ACHSE MIT (SEC). ZEICHNE TH MIT FARBE 6 UND THD MIT FARBE 2 PLOT 'XAXIS'=T, 'XTAG'='(SEC)', TH, 'TYPE'=6, THD, 'TYPE'=2 $ PAUSE ___________________________________________________________________________ DAS VIERTE PLOT-KOMMANDO PLOT TH, 'TYPE'=6, THDD, 'TYPE'=4 BEDEUTET: ZEICHNE TH MIT FARBE 6 UND THDD MIT FARBE 4. X-ACHSE BLEIBT T, DA NICHTS NEUES VEREIBART WURDE. BESCHRIFTUNG DER X-ACHSE BLEIBT (SEC). <?page no="156"?> 142 PLOT TH, 'TYPE'=6, THDD, 'TYPE'=4 $ PAUSE ___________________________________________________________________________ DIE AUSGABEEINHEIT FUER DAS PLOT-KOMMANDO IST IM PARAMETER PLT ADDRESSIERT. DIESER PARAMETER KANN EBENFALLS MIT D AM BILDSCHIRM BETRACHTET UND MIT S GEAENDERT WERDEN. D PLT PLT 0 DIE NACH PLT ANGEGEBENE ZAHL IST DIE CODEADRESSE DER AKTUELLEN PLOT-AUSGABEEINHEIT. WEITERE CODEADRESSEN SIND ZUM BEISPIEL: 3 FUER SCHNITTSTELLE COM1 4 FUER SCHNITTSTELLE COM2 AN COM1 UND COM2 KOENNEN PLOTTER ANGESCHLOSSEN WERDEN. ___________________________________________________________________________ DAS MODELL WIRD NUN UNTER VERSCHIEDENEN BEDINGUNGEN SIMULIERT. UM DIE AUSGABE ZU VEREINFACHEN WERDEN ZUNAECHST ZWEI PROZEDUREN DEFINIERT. SIE ERHALTEN DIE NAMEN STRIP UND PHASE . WENN DIE PROZEDUR MIT IHREM NAMEN AUFGERUFEN WIRD, ERZEUGT SIE JEDESMAL DIE SELBEN ZEICHNUNGEN MIT DEN AKTUELLEN WERTEN DES PREPAR-FILES. PROCED STRIP SET STRPLT=.T., CALPLT=.F. PLOT 'XAXIS'=T, 'XTAG'='(SEC)', TH, 'TYPE'=6, THD, 'TYPE'=2 $ PAUSE PLOT TH, 'TYPE'=6, THDD, 'TYPE'=2 $ PAUSE END <?page no="157"?> 143 PROCED PHASE SET STRPLT=.F., CALPLT=.T. PLOT 'XAXIS'=TH, 'XTAG'='(RAD)', THD,'TYPE'=6, 'TAG'='(RAD/ SEC)' $ PAUSE END ___________________________________________________________________________ NUN SOLL DAS MODELL MIT DEM NEUEN STARTWINKEL THIC=3.0 SIMULIERT WERDEN. FUER DIE PLOTS WIRD DESHALB AUCH EINE NEUE UEBERSCHRIFT VEREINBART. S THIC=3.0, TITLE(11)='STARTWINKEL=3.0' ___________________________________________________________________________ DER NEUE LAUF UEBERSCHREIBT DAS PREPAR-FILE MIT DEN NEU BERECHNETEN WERTEN. START T 0.0 TH 3.00000000 THD 0.0 T 0.50000000 TH 2.62318000 THD-1.93472000 T 1.00000000 TH 0.21616000 THD-7.61792000 T 1.50000000 TH-2.22702000 THD-1.21545000 T 2.00000000 TH-1.22168000 THD 5.53437000 T 2.50000000 TH 1.68442000 THD 3.03100000 T 3.00000000 TH 1.27391000 THD-4.60479000 T 3.50000000 TH-1.44390000 THD-3.16180000 T 4.00000000 TH-1.01807000 THD 4.65704000 T 4.50000000 TH 1.41089000 THD 2.26764000 T 5.00000000 TH 0.58004200 THD-5.04921000 ___________________________________________________________________________ DER AUFRUF VON STRIP UND PHASE ERZEUGT NUN DIE IN DEN PROZEDUREN VEREINBARTEN ZEICHNUNGEN. DIESE ZEIGEN UNS DEN VERLAUF VON GROSSEN SCHWINGUNGEN. ZUERST WIRD STRIP AUFGERUFEN. STRIP <?page no="158"?> 144 ___________________________________________________________________________ MIT PHASE WIRD DANN DIE ZWEITE PROZEDUR AUFGERUFEN. PHASE <?page no="159"?> 145 DAS VERHALTEN DES MODELLS BEI UEBERSCHLAEGEN SOLL NUN BETRACHTET WERDEN. DAZU WIRD GESETZT DIE STARTWINKELGESCHWINDIGKEIT THDIC=10.0, DER STARTWINKEL THIC =0.0 . S THIC=0, THDIC=10, TITLE(11)='STARTGESCHWINDIGKEIT=10.0' ___________________________________________________________________________ DAS SEITHERIGE ENDE DER SIMULATION WIRD MIT D BETRACHTET UND MIT S VERLAENGERT, UM DEM PENDEL ZEIT ZUM EINSCHWINGEN ZU GEBEN. MIT NCIOUT=40 WIRD DIE OUTPUT-TABELLE VERKUERZT, ES WIRD NUR. JEDER 40STE WERT DER PREPAR-LISTE AUF DEM BILDSCHIRM ANGEZEIGT. D TSTP TSTP 4.99000000 S TSTP=9.9 S NCIOUT=40 ___________________________________________________________________________ START T 0.0 TH 0.0 THD 10.0000000 T 1.00000000 TH 7.32159000 THD 7.90605000 T 2.00000000 TH 12.3935000 THD 8.17762000 T 3.00000000 TH 14.9515000 THD-2.28093000 T 4.00000000 TH 10.3747000 THD-0.06594160 T 5.00000000 TH 14.4386000 THD 1.06887000 T 6.00000000 TH 10.8904000 THD-0.84339200 T 7.00000000 TH 14.0923000 THD-0.20574300 T 8.00000000 TH 11.2787000 THD 1.71899000 T 9.00000000 TH 13.4307000 THD-3.27041000 T 9.90000000 TH 11.9206000 THD 3.47950000 __________________________________________________________________________ FUER DIE ZEICHNUNGEN WERDEN NUN NOCH EINIGE VEREINBARUNGEN NEU GETROFFEN. S XINSPL=8.0 STRIP-PLOT-X-ACHSE 5 INCH LANG S XTISPL=2.0 MARKIERE X-ACHSE ALLE 2 INCH S PSFSPL=0.66 VERKLEINERE Y-ACHSENMASSSTAB AUF 66% DIE VERKLEINERUNG ERLAUBT ES DREI DIAGRAMME AUF DEN BILDSCHIRM ZU ZEICHNEN. S STRPLT=.T., CALPLT=.F. $ ZEICHNE STRIPPLOTS, KEINE LINIENPLOTS PLOT 'XAXIS'=T, 'XTAG'='(SEC)', TH, THD, THDD <?page no="160"?> 146 ___________________________________________________________________________ ABSCHLIESSEND SOLL NOCH EINE PARAMETERSTUDIE DURCHGEFUEHRT WERDEN. DER EINFLUSS DER DAEMPFUNG SOLL GEZEIGT WERDEN. ALS ERSTES STARTET DAS PENDEL BEI THIC=1.0 OHNE DAEMPFUNG. MIT OUTPUT 'CLEAR' WIRD VEREINBART, DASS KEINE WERTE WAEHREND DES LAUFS AUF DEM BILDSCHIRM ERSCHEINEN. S TITLE(11)='PARAMETER SWEEP - KDAMP=0.00, 0.02, 0.20' S THIC=1.0, THDIC=0.0, S TSTP=4.99 $ SIMULATIONSENDE BEI 4.99 S KDAMP=0.0 $ DAEMPFUNG=0.0 OUTPUT 'CLEAR' START ___________________________________________________________________________ ES SOLLEN NUN ZWEI WEITERE LAEUFE MIT ANDERER DAEMPFUNG DURCHGEFUEHRT UND IN EINEM DIAGRAMM VERGLICHEN WERDEN. MIT NRWITG=.T. WIRD VEREINBART DAS PREPAR-FILE NICHT ZU UEBERSCHREIBEN, SONDERN DIE NEUEN SIMULATIONSERGEBNISSE ANZUFUEGEN. <?page no="161"?> 147 S NRWITG=.T. S KDAMP=0.02 START S KDAMP=0.20 START NEUE VEREINBARUNGEN FUER DIE ZEICHNUNGEN S CALPLT=.T., STRPLT=.F. $ ZEICHNE LINIENPLOTS, KEINE STRIP-PLOTS S SYMCPL=.T. $ ZEICHNE SYMBOLE AN JEDE KURVE S FTSPLT=.T. $ UNTERDRUECKT DIE LINIE VOM ENDE DER EINEN KURVE ZUM ANFANG DER NAECHSTEN KURVE S XINCPL=5.0 $ LINIENPLOT-X-ACHSE 5 INCH LANG MIT EINEM PLOTKOMMANDO KOENNEN NUN ALLE DREI KURVEN IN EIN DIAGRAMMM GEZEICHNET WERDEN. PLOT TH, 'CHAR'='1' MIT CHAR=1 WIRD VEREINBART, DASS DIE ERSTE KURVE DAS SYMBOL 1 ERHAELT. PLOT TH, 'CHAR'='1' ___________________________________________________________________________ MIT STOP WIRD DIE RUNTIME-SITZUNG BEENDET UND DOS MELDET SICH WIEDER AM BILDSCHIRM STOP <?page no="162"?> 148 <?page no="163"?> 149 5 Anhang: Aktive Federung mit Matlab Simulink In Kapitel 2.4 wurde beschrieben, wie man mit Matlab Simulink das Verhalten dynamischer Systeme untersuchen kann, die sich durch ein Differenzialgleichungssystem beschreiben lassen. Drei Beispiele wurden behandelt und deren Ergebnisse diskutiert: ein Spannungsmesser, ein nichtlineares Pendel und ein Zweimassenschwinger (vereinfachtes Fahrzeugmodell) Das Modell des Zweimassenschwingers soll nun nochmals erweitert werden im Sinne eines mechatronischen Gesamtsystems, bestehend aus mechanischen Bauteilen, Aktuatoren, Sensoren und Steuergerät. Es entsteht ein Zeimassenschwinger mit integrierter aktiver Federung. Grundlagen: Im ursprünglichen Modell gemäß Abbildung 58 auf Seite 71 ist neben dem Dämpfer nur eine konventionelle Stahlfeder zwischen Aufbau und Rad vorgesehen (passives Federungssystem). Die Kraft in der Feder ist eine alleinige Funktion des Relativweges zwischen Aufbau und Rad (Ein-/ Ausfederung). So ist es bei allen passiv gefederten Fahrzeugen: sie nicken deshalb beim Bremsen, wanken in der Kurve nach außen und verlieren an Höhe bei Beladung, und zwar jeweils so weit, bis die Kräfte in den Federn im Gleichgewicht stehen zu den physikalisch bedingten äußeren Kräften und Momenten, die bei diesen Fahrmanövern am Aufbau wirken. Die Erfahrung zeigt, dass alle Insassen eines Fahrzeugs es als sehr angenehm empfinden, wenn derartige dynamischen Bewegungen des Aufbaus unterdrückt werden. Mit dem passiven Federungssystem ist dies jedoch prinzipiell nicht möglich. In Sportwagen versucht man diese Effekte mit harten Federn zu minimieren; dies bedeutet jedoch Verlust an Komfort bei Fahrten über schlechte Straßen. Selbst ein verstellbarer Dämpfer ist dazu nur bedingt in der Lage. Die Luftfeder erlaubt zumindest eine langsame Niveauregelung bei Beladung. Richtig erfolgreich ist die Bekämpfung dieser Bewegungen erst, wenn sie durch ein aktives Federungssystem statisch und dynamisch kompensiert werden. Im Vergleich zur passiven Federung bedeutet dies jedoch viele zusätzliche Bauteile, erhöhte Kosten und zusätzliche Energie vom Verbrennungsmotor. Denn typisch für alle aktiven Systeme ist, dass sie Energie benötigen, um aktiv (nicht reaktiv) Kräfte zwischen den einzelnen Massen des Fahrzeugs aufzubauen. Abhängig vom Typ und Einbauort der aktiven Stellelemente kann entweder nur der Wankfreiheitsgrad des Aufbaus beeinflusst werden (z.B. bei BMW und Porsche) durch einen aktiven Stabilisator an jeder Achse, oder zusätzlich auch der Nickfreiheitsgrad und der vertikale Hub des Aufbaus (z.B. Mercedes S-Klasse), wenn an jedem Rad die passive Feder durch ein aktives Federbein ersetzt wird, System ABC, Active Body Control. Der Stabilisator kann dann entfallen. Aktives Federbein: Im Gegensatz zur passiven Feder ist das aktive Federbein eine Kombination aus einer aktiven Feder und einem konventionellen Dämpfer in einer <?page no="164"?> 150 Baueinheit. Dabei besteht die aktive Feder aus einer konventionellen Stahlfeder in Reihe zu einem hydraulischen Stellzylinder, wie in Abbildung 105 dargestellt: Abbildung 105: Ersatzschaltbild des Zweimassenschwingers mit aktiver Federung Wird die Feder wegen erhöhter Kraft F_ZL zusammengedrückt (rechts), dann fährt der Kolben des Zylinders um das entsprechende Maß x_K aus. Dadurch bleibt die Länge des Federbeins konstant und damit auch der Abstand zwischen Rad und Aufbau. (Beispiel: Beim Bremsen müssen die beiden Zylinder an der Vorderachse ausfahren, die an der Hinterachse einfahren.) Modell in Matlab Simulink Das Modell des passiven Zeimassenschwingers, Abbildung 61 auf Seite 74, muss um diesen hydraulischen Stellzylinder erweitert werden. Dieser wird im Modell des aktiven Zweimassenschwingers ZMS_aktiv gemäß Abbildung 106 auf Seite 151 in einem Subsystem Federbein integriert. (Die hydraulische Versorgung ist nicht nachgebildet, z.B. die Pumpe und Druckregelventile.) Proportional zur gemessenen Relativbewegung x_rel zwischen Rad und Aufbau (Subsystem x_Rel_Sensor) fließt ein Ölstrom (Subsystem Regelalgorithmus) in den Stellzylinder des Federbeins und führt zu einer Kolbenbewegung x_Aktuator. Dadurch bleibt der Abstand zwischen Aufbau und Rad konstant (Mechanik). Das ist das Ziel der Regelung in diesem geschlossenen Regelkreis. Rad m_R Aufbau m_A Feder c_A Dämpfer d_A x_S x_R x_A Reifen c_R F_ZL passiv Rad m_R Aufbau m_A Feder c_A Dämpfer d_A x_S x_R x_A Reifen c_R F_ZL aktiv Stellzylinder x_K <?page no="165"?> 151 Abbildung 106: Modell des Zweimassenschwinger mit aktiver Federung Die notwendige Erweiterung im Subsystem Federbein ist in Abbildung 107 auf Seite 152 dargestellt. Die Federkraft F_F ist proportional zur Längenänderung x_F der Feder. Um den hydraulischen Druck p in bar im Aktuator und die zur Verstellung notwendige Leistung P = Q_Soll * p_hyd in Watt mit Q_Soll > 0, realisiert über einen Saturation-Block -, richtig berechnen zu können, müssen • die achskinematischen Verhältnisse berücksichtigt werden, nämlich ein Übersetzungsverhältnis i_Achse=2 zwischen Rad und Einbauort des Federbeins • und die Gesamtkraft force in der Feder berechnet werden, unter Einbeziehung des Gewicht des Aufbaus static load. Aus der Gesamtkraft force in der Feder berechnet sich der im Aktuator herrschende Hydraulikdruck p_hyd durch Division durch die Kolbenfläche. Zur Berechnung der benötigten hydraulischen Leistung dürfen nur positive Ölströme berücksichtigt werden, also Öl, das von der Pumpe in den Zylinder fließt. (Negative Ölströme fließen dann, wenn der Kolben einfährt. Das dabei verdrängte Öl fließt in Richtung Ölbehälter und belastet die Pumpe nicht. ) <?page no="166"?> 152 Abbildung 107: Federbein des aktiven Systems Zur Berechnung der Längenänderung der Feder wird das Subsystem Länge der Feder eingeführt, Abbildung 108. Die Änderung ist nicht (x_R - x_A), wie zuvor im Falle der passiven Federung, sondern (x_R - x_A + x_Aktuator). Der Weg des Aktuators in Metern berechnet sich durch Integration des Ölstroms Q_Soll in l/ min. Abbildung 108: Subsystem zur Berechnung der Federlänge <?page no="167"?> 153 Simulation: Die Simulation soll das Zusammenspiel aller Komponenten zeigen, wenn im Manöver 5 die Zusatzlast F_ZL über der Zeit linear wächst und nach t = 1s einen konstanten maximalen Wert von F_ZL=1000 N erreicht, siehe Abbildung 109. (In einem realen Fahrzeug könnte das z.B. am linken Vorderrad beim Einlenken des Lenkrades wegen einer Rechskurve sein. Zur Simulation der Verhältnisse am kurveninneren Vorderrad rechts müsste die Zusatzlast entsprechend negativ gewählt werden.) Die Ergebnisse sind ab Abbildung 110 auf Seite 153 dargestellt. Abbildung 109: Manöver 5: Mit der Zeit wachsende Zusatzlast F_ZL Das Fahrzeug wird so weit angehoben, dass trotz dieser Zusatzlast F_ZL der Relativweg x_Rel zwischen Rad und Aufbau nach erfolgter Regelung wieder den Wert Null erreicht. Bei einer Federsteifigkeit von c_A = 160000N/ m drückt sich die Feder zwar um einen Betrag von (i_Achse * F_ZL) / c_A = (2*1000N) / (160000N/ m) = 0,0125m zusammen. Der Kolben des Aktuators fährt aber genau um dieses Maß aus und hält den Abstand zwischen Rad und Aufbau konstant, Abbildung 110: Abbildung 110: Manöver 5: Kompensation der Federverkürzung Dies bedeutet aber nicht, dass der Aufbau die Höhe absolut einhält. Er bewegt sich wegen der erhöhten Last gemeinsam mit dem Rad nach unten, und zwar um ein Stück x_A = x_Rad = - F_ZL / c_R = 1000N / (200000N/ m) = 5mm. Die Reifen- <?page no="168"?> 154 einfederung zwischen Rad und Straße wird nicht berücksichtigt im Regelalgorithmus und deshalb nicht kompensiert, Abbildung 111. (Im realen Fahrzeug ist deshalb der absolute Wankwinkel nicht Null, sondern nur der Wankwinkel relativ´zu den Rädern.) Abbildung 111: Manöver 5: Absolutbewegung von Rad und Aufbau Der hdraulische Druck im Aktuator steigt von anfangs ca. 68,9 bar (dem statischen Druck ohne Zusatzlast) auf ca. 86,4 bar (mit Zusatzlast) an, Abbildung 112. Der Ölstrom Q_Soll in den Kolben ist proportional zur Regelabweichung, dem Relativweg x_Rel. Sein Maximalwert beträgt ca. 0,69 l/ min. Abbildung 112: Manöver 5: Druck im und Ölstrom in den Hydraulikzylinder Die hydraulische Leistung P berechnet sich aus dem Produkt von Hydraulikdruck p und Ölstrom Q_Soll. Der Spitzenwert beträgt P = 101 Watt, Abbildung 113 auf Seite 155. Dies ist auf den ersten Blick nicht viel, es handelt sich dabei allerdings nur um den Leistungsanteil, den das Federbein zum Verstellen des Kolbens benötigt. (Der Versorgungsdruck im realen System ist viel größer, p_System = 180 bar, um ausreichend Stelldynamik für alle Regelaufgaben zu erreichen. Der überschüssige Leistungsanteil ist reine Verlustleistung, das Öl wird auf dem Weg zum Tank abgedrosselt auf p_Tank = 1 bar. Die in das System gesteckte Energie wird dabei in diesen hydraulischen Drosseln in Wärme umgewandelt.) <?page no="169"?> 155 , Abbildung 113: Manöver 5: Vom System benötigte hydraulische Leistung Damit ist die Simulation des aktiven Zweimassenschwingers beendet. Zur Vertiefung der Kenntnisse in Simulationstechnik und, um weitere Erfahrungen zu sammeln im Umgang mit Simulationstools, ist es empfehlenswert, die im Buch gezeigten Beispiele nochmals selbst zu rechnen, unterschiedliche Parametervarianten durchzuspielen und jeweils die Auswirkung der verschiedenen Werte auf den Zeitverlauf der Zustandsgrößen zu prüfen. Hierzu gibt es im Download-Bereich des Expert Verlags unter http: / / www.expertverlag.de/ im Unterkapitel „Verschiedenes“ die im Buch behandelten Matlab Simulink-Beispiele als Download. Voraussetzung für die Arbeit mit diesen Modellen ist allerdings eine Matlab Simulink-Lizenz. Dabei kann man üben, wie man vorgeht bei der Modellerstellung mit Simulink-Blöcken verschiedene Integrationsverfahren einsetzen und vergleichen die unterschiedlichen Ausgabemöglichkeiten testen - Systemparameter variieren und vor allem trainieren, das dynamische Verhalten zu interpretieren <?page no="170"?> 156 Darüber hinaus steht auch ein weiteres einfaches Beispiele zum Download zur Verfügung, das Modell eines eines hüpfenden Balls. Es sind aber auch Demonstrationsbeispiele auf den Internet-Seiten von Matlab direkt wählbar unter http: / / www.mathworks.com/ industries/ auto/ demos.html <?page no="171"?> 157 6 Weiterführende Literatur Die folgende Liste zeigt eine kleine Auswahl von Büchern, die geeignet sind, das Grundlagenwissen über Simulationstechnik zu ergänzen und zu vertiefen: Ammon, D.: Modellbildung und Systementwicklung in der Fahrzeugdynamik, Teubner, 1997. Bossel, H.: Systeme, Dynamik, Simulation, Books on Demand GmbH, 2004. Deger, Y.: Die Methode der Finiten Elemente, Expert, 5. Aufl., 2008. Feindt, B. E.-G.: Computersimulation von Regelungen, Oldenbourg, 1999. Gipser, M.: Systemdynamik und Simulation, Teubner, 1999. Karnopp, D., Margolis, D., Rosenberg, R.: System Dynamics, Wiley, 3. Aufl., 2000. Scherf, H.E.: Modellbildung und Simulation dynamischer Systeme, Oldenbourg, 2003. Stelzmann, U., Groth, C., Müller, G.: FEM für Praktiker - Band 2: Strukturdynamik, Expert, 5. neu bearb. Auflage, 2008. The MathWorks Inc: Matlab Simulink Handbücher (für Lizenzinhaber) <?page no="172"?> 158 <?page no="173"?> 159 7 Sachregister Ablaufkontrolle 92, 95, 108 GPSS Seize Block 118 ACSL 48, 49 GPSS Standard Numeric 129 ACSL Modell-Statements 50, 52, 135 Attribute ACSL Run-Time-Befehle 49, 50, 52-56, GPSS Storage 118 136 GPSS Storage 118 ACSL Schrittweitensteuerung 57 GPSS Table 124 ACSL Systemvariable 56, 57 GPSS Tabulate Block 123 ACSL-Programm-Sections 50, 51, 52, 54 GPSS Terminate Block 112, 113 Analogiemodelle 5, 6 GPSS Test Block 128 Analogrechner 6, 7, 9 GPSS Transaktionen 112 Attribute 2 GPSS Transfer Block 115 Betriebsstrategien 106 Häufigkeit 95, 100 Humanoid-Dummies 23, 27 Crash-Simulation 17, 19 Hydropuls-Anlage 29 Datenbereiche 92 Insassenklima-Simulation 26, 27 Dauerfestigkeitssimulation 20, 22 Integration numerisch 44 Detaillierung 12 Integrationsverfahren 40, 45, 46, 67 Detaillierungsgrad 2 Differenzialgleichung 11, 41-43 Klimakammer 17 Komponenten 2 Elemente 2 Ergebnis-Animation 20, 21, 22 MatlabSimulink 48, 62 Ergebnis-Interpretation 60 MatlabSimulink Compare To 86 Ergonomie 23-25 Zero Block Erlang-k-Verteilung 104 MatlabSimulink From-Block 75 Erwartungswert 97 MatlabSimulink Goto-Block 74 Exponentialverteilung 101, 102, 109 MatlabSimulink Manual Switch 75 MatlabSimulink Saturation 84, 151 Fahrdynamiksimulation 20, 21 Block Fahrsimulator 32-36 MatlabSimulink Scope 73 MatlabSimulink Signal Routing 74, 75 Gauss-Verteilung 105 MatlabSimulink Sine Wave 78 Geräuschprüfstand 17,19 Block Getriebesimulator 30, 31 MatlabSimulink Step-Block 76 Gleichverteilung 102, 103 MatlabSimulink Subsystem 73 Glockenkurve 105 MatlabSimulink To Workspace 81 GPSS 112 Block GPSS Advance Block 113 Modell 1, 4 GPSS Depart Block 122 Modell Validierung 38 GPSS Enter Block 118 Modellbildung 12, 15 GPSS Facility 118 Modelle mathematisch 10, 11, 18 GPSS Generate Block 112 Modelle physikalisch 5, 17 GPSS Leave Block 118 GPSS Mark Block 123 Normalverteilung 105 GPSS Queue Block 122 GPSS Release Block 118 Pendel, nichtlinear 58-61, 69 <?page no="174"?> 160 Poisson-Verteilung 101 Systementwurf 37 Pseudo-Zufallszahlen 102 Systemoptimierung 37 Systempostulierung 37 Regen-Simulation 30 Systemstudie 1, 39 Schwerelosigkeit-Simulation 30 Varianz 94, 99 Simulation 16 Verteilungsfunktion 82, 83, 97, Simulation ko- 20 98,101 Simulation offline 18, 20 Simulation Sprachen 47 Wahrscheinlichkeit 94, 97 Simulationstechnik 1 Wahrscheinlichkeitsdichte 98 Simulink 62 Warteschlange 106 Simulink Programmerstellung 65-67 Warte-Verlust-System 106, 107 Spannungsmessinstrument 63 Weltraumsimulator 30, 31 Summenhäufigkeit 96, 97, 100 Windkanal 17, 18 System 1, 2 Wirkungsplan 63, 64 System kontinuierlich 38, 41, 91 System zeitdiskret 38, 91 Zufallsvariable 74 Systemanalyse 37 Zufallsvariable diskrete 74 Systeme gemischte 28 Zufallsvariable stetige 78