DE10023377C2 - Verfahren zur Erhöhung der Leistungsfähigkeit einer Computereinrichtung bei Finite-Elemente-Simulationen und eine solche Computereinrichtung - Google Patents

Description

Die Erfindung betrifft ein Verfahren zur Erhöhung der Leistungsfähigkeit einer Computereinrichtung bei Finite-Elemente-Simulationen durch effiziente automatische Konstruktion geeigneter Basisfunktionen für die Berechnung von Näherungslösungen und eine solche Computereinrichtung.

Eine Vielzahl technisch-physikalischer Phänomene lässt sich durch partielle Differenzialgleichungen beschreiben. Hierunter fallen u. a. Probleme aus der Strömungsmechanik (z. B. Strömung um eine Tragfläche), der elektromagnetischen Feldtheorie (z. B. elektrischer Feldverlauf in einem Transistor) oder der Elastizitätstheorie (z. B. Deformation einer Autokarosserie). Die genaue Kenntnis und Beschreibung derartiger Vorgänge ist ein zentraler Bestandteil bei der Konstruktion und Optimierung technischer Objekte. Um zeit- und kostenintensive Experimente einzusparen ist das Interesse an computergestützten Simulationen groß. Hierbei haben sich Finite-Elemente-Verfahren (FE-Verfahren) etabliert und sind seit langem Gegenstand intensiver Forschung. Dies betrifft insbesondere auch automatische Netzgenerierungsverfahren als Grundlage für die Konstruktion geeigneter Basisfunktionen.

Fig. 1 veranschaulicht in der linken Bildhälfte den Stand der Technik beim Ablauf einer FE-Simulation für ein lineares Randwertproblem als typisches Modellbeispiel. Ausgehend von den Daten, die die Geometrie des zu simulierenden technischen Objektes beschreiben, wird zunächst ein System von Basisfunktionen konstruiert, das einerseits die Einhaltung der Randbedingungen ermöglicht und andererseits zur Approximation der gesuchten Lösung geeignet ist. Dann wird mit Hilfe dieser Basisfunktionen ein lineares Gleichungssystem mittels numerischer Integrationsverfahren aufgestellt. Schließlich werden die Koeffizienten der gesuchten Approximation als Lösung dieses linearen Gleichungssystems ermittelt.

FE-Verfahren oder deren Anwendung sind Gegenstand einer Reihe von Patenten. Beispielhaft sei US 4,819,161 genannt, wo die FE-Approximationen einer großen Klasse von Differenzialgleichungen automatisiert werden, sowie US 5,731,817, wo ein Verfahren zur Generierung von Hexaeder-Netzen die Grundlage für ein Verfahren zur FE-Simulation bildet.

Bei den meisten praxisrelevanten FE-Verfahren werden Basisfunktionen verwendet, die auf einer durch Vernetzung des Simulationsgebietes erzeugten Zerlegung definiert sind. Fig. 2a zeigt eine Auswahl gebräuchlicher Elemente deren Dimension, Grad, Glattheit und Parameter in Fig. 2b aufgelistet sind. Eine Übersicht über Vernetzungsmethoden ebener Gebiete findet sich beispielsweise in K. Ho-Le, Finite element mesh generation methods: a review and classification, Comp. Aided Design 20 (1988), 27-38. Die Vernetzung komplizierterer dreidimensionaler Gebiete gestaltet sich nach derzeitigem Kenntnisstand überaus schwierig, wie S. Owen, A survey of unstructured mesh generation technology, Proceedings, 7th International Meshing Round Table, Sandia National Lab (1998), 239-257, zeigt. Die Verfahren benötigen lange Rechenzeiten und sind teilweise auch noch nicht vollkommen automatisiert. Es gibt jedoch in jüngerer Zeit eine Reihe sehr innovativer neuer Ansätze. Beispielsweise wird in A. Fuchs, Optimierte Delaunay-Triangulierungen zur Vernetzung getrimmter NURBS-Körper, Universität Stuttgart, 1999 eine Kräfteverteilung simuliert, um eine optimale Verteilung von Triangulierungspunkten zu erzielen. In US 5,729,670 werden zwei- und dreidimensionale Netze durch Lösen von Strömungsproblemen erzeugt, was eine interessante Umkehrung des üblichen FE-Mechanismusses ist. Daneben wurden viele Algorithmen entwickelt, um einzelne Teilaspekte der Netzgenerierung zu verbessern. Beispielsweise wird in DE 196 21 434 A1 und US 5,774,696 ein Verfahren zur Elimination von Überschneidungen mit vorgegebenen Kanten oder Grenzflächen bei Delaunay-Triangulierungen beschrieben.

Vernetzungsfreie FE-Verfahren haben bisher keine Bedeutung für die Anwendungen erlangt. Sowohl bei dem Lagrange-Multiplier-Verfahren, s. z. B. J. H. Bramble, The Lagrange multiplier method for Dirichlet's problem, Math. Comp. 37 (1981), 1-11, als auch bei der Penalty-Methode, s. z. B. P. Bochev und M. Gunzburger, Finite element methods of least squares type, SIAM Review 40 (1998), 789-837, stellt bei Verwendung einfacher, stabiler Basisfunktionen die Behandlung von Randbedingungen ein Hauptproblem dar.

Bei vielen technischen Simulationen ist die automatische Netzgenerierung sehr aufwendig und nimmt den weitaus größten Teil der Rechenzeit in Anspruch. Des Weiteren ist die Approximationskraft der üblicherweise verwendeten linearen und multilinearen Basisfunktionen gering. Um genaue Resultate zu erzielen, muss also eine große Anzahl von Basisfunktionen verwendet und damit ein entsprechend großes Gleichungssystem gelöst werden. Ansatzfunktionen höheren Grades über Triangulierungen besitzen in der Regel ebenfalls ein ungünstiges Verhältnis zwischen der erreichbaren Genauigkeit und der Anzahl der verwendeten Basisfunktionen. Schließlich lassen sich auf unstrukturierten Netzen keine glatten Basisfunktionen auf einfache Weise definieren. Bereits für stetig differenzierbare Elemente sind sehr spezielle Konstruktionen notwendig (s. Fig. 2a).

HÄUSSLER-COMBE, U., u. a.: "Die elementfreie Galerkin Methode - Überblick und Anwendungsbeispiele", Bauingenieur, Band 74, Nummer 3, März 1999, Seite 130 bis 138 beschreibt ein Verfahren, bei dem eine Familie von Gewichtsfunktionen dazu verwendet wird, Polynome zu lokalisieren und damit eine Basis zu konstruieren. Randbedingungen werden mittels des sogenannten Penalty-Ansatzes berücksichtigt, wobei die Wahl der Konstanten erhebliche Probleme aufwirft. Die Wahl der Gewichtsfunktion gemäß Bild 3 dieses Dokuments dient zum Modellieren eines Risses, und nicht - wie gemäß der vorliegenden Erfindung - zur Einhaltung von Randbedingungen. Gemäß der vorliegenden Erfindung hat die Gewichtsfunktion eine grundsätzlich andere Bedeutung, nämlich die, dass homogene Randbedingungen per Konstruktion exakt eingehalten werden und somit gerade der in diesem Dokument angegebene Penalty- Ansatz überflüssig wird und die damit verbundenen Probleme beseitigt werden. Die Lokalität der Basis ist gemäß der vorliegenden Erfindung durch der Verwendung der B-Splines gesichert.

BARTELS, R. H., u. a.: "An Introduction to Splines for Use in Computer Graphics & Geometric Modeling", MORGAN KAUFMANN PUB., 1987, Seite 46 bis 63, stellt einige Grundzüge der B-Spline Theorie dar. Ein konkreter Bezug zur Finite-Elemente-Methode ist aber nicht offenbart oder nahegelegt.

V. L. Rvachev und T. I. Sheiko: "R-Functions in boundary value problems in mechanics", Applied Mechanics Reviews, 48 (1995), Seiten 151 bis 188 und V. Shapiro: "Theory of R-Functions and Applications: A Primer", Cornell Programmable Automation, Sibley School of Mechanical Engineering, Ithaca NY (1988), CPA 883, geben Verfahren an, wie Gewichtsfunktionen für Gebiete konstruiert werden können, die sich speziell durch Schnitt, Vereinigung und Komplementbildung aus elementaren Teilgebieten erzeugen lassen. Diese Gewichtsfunktionen können zur Einhaltung von wesentlichen Randbedingungen verwendet werden.

V. L. Rvachev und T. I. Sheiko, V. Shapiro und I. Tsukanov:
"On Completeness of RFM-Solution Structures", Computational Mechanics, 25 (2000), Seiten 305 bis 316 und V. Shapiro und I. Tsukanov: "Meshfree Simualation of Deforming Domains", Computer-Aided Design, 31 (1999), Seiten 459 bis 471 setzen zwar B-Splines als Basisfunktionen ein, das Stabilitätsproblem bleibt aber ungelöst.

P. Oswald: "Multilevel Solvers for elliptic problems on domains", Multiscale Wavelet Methods for PDEs, W. Dahmen et al. (Herausgeber), Academic Press (1996), Seiten 3 bis 58, beschreibt ein Verfahren, wie B-Splines durch eine Erweiterungskonstruktion stabilisiert werden können. Die dabei verwendete Konstruktion unterscheidet sich aber von der erfindungsgemäßen Konstruktion. Darüber hinaus werden wesentliche Randbedingungen bei der Konstruktion nicht berücksichtigt.

Der Erfindung liegt das Problem zugrunde, die Leistungsfähigkeit bekannter FE-Verfahren und die FE-Verfahren ausführender Computereinrichtungen durch die effiziente Konstruktion von Basisfunktionen mit günstigen Eigenschaften zu erhöhen. Insbesondere soll die Vernetzung des Simulationsgebietes vollständig entfallen, eventuell gegebene Randbedingungen eingehalten werden, genaue Lösungen mit relativ wenigen Koeffizienten erreichbar und das entstehende Gleichungssystem effizient lösbar sein. Dadurch sollen die Nachteile des Standes der Technik überwunden und somit die Genauigkeit und Geschwindigkeit der Simulation physikalischer Eigenschaften bei der Konstruktion und Optimierung technischer Objekte verbessert werden.

Einige zentrale Begriffe und Bezeichnungen, die in der nachfolgenden Beschreibung des erfindungsgemäßen Verfahrens verwendet werden, sollen vorab erläutert werden.

Das Simulationsgebiet Ω ist eine beschränkte Menge der Dimension d = 2 oder d = 3, auf der die zu untersuchenden physikalischen Größen mittels der FE-Methode approximiert werden sollen. Der Rand des Simulationsgebietes wird mit T bezeichnet. Unter einem Gitter mit Gitterweite h versteht man eine Zerlegung einer Teilmenge der Ebene oder des Raumes in Gitterzellen Z_k. Dabei ist jede Gitterzelle je nach Dimension d ein Quadrat oder ein Würfel mit Kantenlänge h. Genauer gilt Z_k = kh + [0, h]^d, wobei k eine Menge von ganzzahligen d-Vektoren durchläuft. Die uniformen Tensorprodukt B-Splines in d Variablen vom Grad n zur Gitterweite h werden mit b_k bezeichnet, s. z. B. C. de Boor, A Practical Guide to Splines, Springer, 1978. Dies sind (n-1)-mal stetig differenzierbare Funktionen, die auf den Gitterzellen mit Polynomen vom Grad n übereinstimmen, wie in Fig. 4 dargestellt. Die obere Bildhälfte zeigt einen Träger Q_k des B-Splines vom Grad n = 2, der Dimension d = 2 und der Glattheit m = 1. In der unteren Bildhälfte ist der resultierende Tensorprodukt B-Spline b_k abgebildet. Der Träger Q_k, also die Vereinigung aller Gitterzellen, auf denen der B-Spline b_k nicht identisch Null ist, besteht aus (n+1)^d Gitterzellen; es gilt genauer Q_k = kh + [0, (n+1)h]^d. In allen Figuren ist der B-Spline b_k an der Stelle kh markiert, also beispielsweise im Fall d = 2 in der linken unteren Ecke des Trägers. Für die FE-Simulation sind nur solche B-Splines von Belang, deren Träger das Simulationsgebiet Ω schneiden; sie werden als relevante B-Splines bezeichnet. Die relevanten B-Splines werden nochmals in zwei Gruppen eingeteilt: diejenigen B-Splines, bei denen der Teil des Trägers im Inneren des Simulationsgebietes größer als eine vorgegebene Schranke s ist, werden als innere B-Splines bezeichnet. Alle übrigen relevanten B-Splines werden als äußere B-Splines bezeichnet.

Die Aufgabe der vorliegenden Erfindung wird durch das im Anspruch 1 bestimmte Verfahren gelöst. Besondere Ausführungsarten der Erfindung sind in den Unteransprüchen bestimmt. Der Anspruch 11 bestimmt eine erfindungsgemäße Computereinrichtung.

In der rechten Bildhälfte der Fig. 1 ist die Einbindung des erfindungsgemäßen Verfahrens in den Ablauf einer FE-Simulation nach dem Stand der Technik bzw. die Substitution einiger Verfahrensschritte einer FE-Simulation nach dem Stand der Technik durch das erfindungsgemäße Verfahren dargestellt.

Die Eingabe 1 des Simulationsgebietes Ω kann über ein Eingabemittel erfolgen, insbesondere auch durch Einspeichern von Daten, die aus einer computerunterstützten Konstruktion (CAD/CAM) abgeleitet sind. Beispielsweise können so die bei der Konstruktion eines Kraftfahrzeugs verwendeten Daten unmittelbar in die erfindungsgemäße FE-Simulation einfließen.

Bei der Eingabe 2 und Abspeicherung des Typs der Randbedingungen wird zwischen natürlichen und wesentlichen Randbedingungen unterschieden. Die erfindungsgemäße Basis wird für homogene Randbedingungen gleichen Typs konstruiert. Insbesondere verschwinden bei wesentlichen Randbedingungen die Basisfunktionen auf dem Rand Γ. Inhomogene Randbedingungen lassen sich dann beim Aufstellen des FE-Systems mittels Methoden, die dem Stand der Technik entsprechen, behandeln.

Schließlich werden Steuerparameter eingelesen 3. Diese betreffen den Grad n und die Gitterweite h der zu verwendenden B-Splines sowie die Schranke s zur Klassifikation der inneren und äußeren B-Splines. Bei fehlender Spezifikation können alle diese Eingabeparameter mit Hilfe von empirisch oder analytisch konstruierten Bewertungsfunktionen automatisch bestimmt werden.

Die folgende erfindungsgemäße Konstruktion der Basisfunktionen gliedert sich in die in Fig. 3 schematisch dargestellten Schritte, die nun beschrieben werden.

Nach dem Einlesen des Simulationsgebietes Ω wird im ersten Verfahrensschritt ein das Simulationsgebiet Ω überdeckendes Gitter erstellt. Dann wird geprüft, welche der Gitterzellen ganz, teilweise oder nicht im Simulationsgebiet Ω liegen, die Zelltypen bestimmt 4 und diese Information über die Zelltypen gespeichert. Hierzu sind im Wesentlichen Innen/Außen-Tests sowie Schnittbestimmungen zwischen dem Rand Γ des Simulationsgebietes Ω und den Strecken bzw. Quadraten, die die Gitterzellen beranden, durchzuführen. Fig. 7 zeigt die Ein- und Ausgabedaten für diesen Verfahrensschritt.

Im zweiten Verfahrensschritt werden zunächst anhand der Information über die Zelltypen die relevanten B-Splines bestimmt. Dann wird die Klassifikation 5 in innere und äußere B-Splines durchgeführt; die entsprechenden Listen der Indizes werden mit I und J bezeichnet. Hierzu ist mit Hilfe der im ersten Verfahrensschritt gewonnenen Daten die Größe derjenigen Teile der Träger der B-Splines zu bestimmen, die im Inneren des Simulationsgebietes liegen und mit der vorgegebenen Schranke s zu vergleichen. Fig. 9 zeigt die Ein- und Ausgabedaten für diesen Verfahrensschritt.

Im dritten Verfahrensschritt werden Kopplungskoeffizienten e_i,j berechnet 6, die innere und äußere B-Splines gemäß der Vorschrift

verbinden. Jedem inneren B-Spline b_i wird damit ein erweiterter B-Spline B_i zugeordnet. Die Konstruktion und die Eigenschaften der Indexmengen J(i) und der Kopplungskoeffizienten e_i,j sind wie folgt gegeben. Die Indexmengen J(i) bestehen aus Indizes äußerer B-Splines. Ihnen entsprechen komplementäre Indexmengen I(j) von Indizes innerer B-Splines; es gehört also i genau dann zu I(j), wenn j zu J(i) gehört. Für einen gegebenen äußeren Index j ist die Indexmenge I(j) ein Feld, also ein quadratischer bzw. würfelförmiger Bereich von (n+1)^d inneren Indizes, der durch einen minimalen Abstand zum Index j gekennzeichnet ist. Für einen gegebenen äußeren Index j und einen inneren Index i in der Indexmenge I(j) sei p_i dasjenige d-variate Polynom vom Grad n in jeder Variablen, das an der Stelle i den Wert 1 und an allen anderen Punkten des Feldes I(j) den Wert 0 hat. Dann ist der Kopplungskoeffizient e_i,j als Wert von p_i an der Stelle j gegeben, es gilt also e_i,j = p_i(j). Die konkreten Werte der Kopplungskoeffizienten können entweder für verschiedene Grade und relative Positionen von j und I(j) tabelliert werden oder aber mit Hilfe von Lagrangepolynomen einfach berechnet werden. Fig. 11 zeigt die Ein- und Ausgabedaten für diesen Verfahrensschritt.

Sind natürliche Randbedingungen gegeben, so werden die in Gleichung (1) definierten erweiterten Splines ohne Modifikation zur weiteren Durchführung des FE-Verfahrens verwendet. Sind dagegen wesentliche Randbedingungen gegeben, so ist noch eine Gewichtung gemäß der Vorschrift

durchzuführen. Die zugehörige Abfrage erfolgt in einem optionalen Verfahrensschritt 6a. Die so definierten Funktionen werden als gewichtete erweiterte B-Splines (WEB-Splines) bezeichnet. Formal gesehen entsprechen die bei natürlichen Randbedingungen verwendeten erweiterten B-Splines dem Sonderfall w(x) = 1. Sie werden deshalb ebenfalls als WEB-Splines bezeichnet. Für den Fall wesentlicher Randbedingungen ist die Gewichtsfunktion w wie folgt gekennzeichnet: Für alle Punkte x des Simulationsgebietes lässt sich w(x) nach oben und unten durch von x unabhängige positive Konstanten mal den Abstand dist(x) des Punktes x vom Rand Γ abschätzen. Mit anderen Worten ist w im Inneren von Ω positiv und verschwindet in der Nähe des Randes Γ so schnell wie die Abstandsfunktion dist. Für Simulationsgebiete, die durch elementare geometrische Objekte (Kreise, Ebenen, Ellipsen, etc.) berandet sind, kann eine geeignete Gewichtsfunktion gegebenenfalls in expliziter analytischer Form angegeben werden. Anderenfalls sind Berechnungsvorschriften heranzuziehen, die typischerweise eine Glättung der Abstandsfunktion darstellen. Der Skalierungsfaktor 1/w(x_i) entsteht durch Auswertung der Gewichtsfunktion im Gewichtspunkt x_i. Dies kann ein beliebiger Punkt im Träger des B-Splines b_i sein, der vom Rand mindestens den halben Schrankenwert s/2 entfernt ist.

Als Resultat des erfindungsgemäßen Verfahrens erhält man eine Berechnungsvorschrift für die WEB-Splines B_i (vgl. Definitionen (1) und (2)), die alle günstigen Eigenschaften gemäß der Aufgabenstellung aufweisen. Damit kann nun das FE-Verfahren gemäß dem Stand der Technik weitergeführt werden. Es ist dabei jedoch möglich, nach Anspruch 11 die reguläre Gitterstruktur der Basisfunktionen vorteilhaft auszunutzen.

Zusammenfassend ist zu sagen, dass für das erfindungsgemäße Verfahren unter anderem die Kopplung der äußeren an die inneren B-Splines von Bedeutung ist. Diese bewirkt, dass die konstruierte Basis die für FE-Berechnungen wesentlichen Eigenschaften besitzt. Insbesondere ist eine erfindungsgemäße Basis B_i (i aus der Indexmenge I) gleichmäßig in Bezug auf die Gitterweite h stabil und besitzt die gleiche Fehlerordnung wie die B-Splines b_k bei Approximation glatter Funktionen, die die gleichen Randbedingungen erfüllen. Zum anderen ist durch die Verwendung der Gewichtsfunktion w die Einhaltung wesentlicher Randbedingungen gewährleistet.

Die zuvor beschriebenen Verfahrensschritte und deren Eigenschaften sowie die bevorzugte Ausführungsart und die zugehörigen Ausführungsbeispiele sollen mit Hilfe der nachfolgend aufgelisteten Figuren weiter erläutert bzw. illustriert werden. Dabei können die in den Ansprüchen und in der Beschreibung erwähnten Merkmale jeweils einzeln für sich oder in beliebiger Kombination erfindungswesentlich sein.

Fig. 1 zeigt zum Stand der Technik die einzelnen Schritte beim Ablauf einer Finite-Elemente-Simulation und ordnet die erfindungsgemäße Ermittlung der WEB-Basis in diesen Ablauf ein,

Fig. 2a stellt einige Finite Elemente des Standes der Technik dem in Fig. 2b gezeigten WEB-Element gegenüber und listet die für Finite-Elemente-Approximationen relevanten Parameter auf,

Fig. 3 zeigt in einem Flussdiagramm die Verfahrensschritte zur Ermittlung der WEB-Basis,

Fig. 4 zeigt einen Träger und den zugehörigen Tensorprodukt-B-Spline vom Grad 2,

Fig. 5 illustriert die Problemstellung des ersten Ausführungsbeispiels (Auslenkung einer Membran unter konstantem Druck) und zeigt die zugehörige Lösung,

Fig. 6 zeigt die Zelltypen anhand des ersten Ausführungsbeispiels,

Fig. 7 gibt eine Übersicht der Ein- und Ausgabedaten des Verfahrens zur Bestimmung der Zelltypen,

Fig. 8 illustriert exemplarisch die Klassifikation der B-Splines anhand des ersten Ausführungsbeispiels,

Fig. 9 gibt eine Übersicht der Ein- und Ausgabedaten des Verfahrens zur Klassifikation der B-Splines,

Fig. 10 zeigt die Kopplungskoeffizienten eines äußeren B-Splines und die zugehörigen inneren B-Splines anhand des ersten Ausführungsbeispiels,

Fig. 11 gibt eine Übersicht der Ein- und Ausgabedaten des Verfahrens zur Berechnung der Kopplungskoeffizienten,

Fig. 12 illustriert die Konstruktion der Gewichtsfunktion der bevorzugten Ausführungsart,

Fig. 13 zeigt den Träger eines WEB-Splines und die zugehörigen Kopplungskoeffizienten anhand des ersten Ausführungsbeispiels,

Fig. 14 erläutert die Problemstellung eines zweiten Ausführungsbeispiels (inkompressible Strömung) und deren Lösung anhand der Strömungslinien und einer Abbildung zur Verteilung der Strömungsgeschwindigkeit,

Fig. 15a bis 15c zeigen für das zweite Ausführungsbeispiel zu den B-Spline-Graden n = 1, 2, 3 und gleicher Gitterweite die sich ändernde Klassifikation der B-Splines,

Fig. 16a bis 16c geben Auskunft über die Fehlerentwicklung bei der Finite-Elemente-Approximation mit Hilfe von WEB-Splines sowie über das Laufzeitverhalten der WEB-Approximation für das zweite Ausführungsbeispiel,

Fig. 17 stellt die WEB-Basis einem auf linearen Ansatzfunktionen über einer Triangulierung (Stand der Technik) basierenden Verfahren gegenüber und

Fig. 18 zeigt eine erfindungsgemäße Computereinrichtung.

Eine besonders günstige Ausführungsart des erfindungsgemäßen Verfahrens, die als WEB-Verfahren bezeichnet wird, ist durch folgende Festlegungen gegeben:
Die Schranke s wird so gewählt, dass die inneren B-Splines dadurch gekennzeichnet sind, dass mindestens eine der Gitterzellen ihres Trägers vollständig im Simulationsgebiet Ω liegt. Da zur Bestimmung der relevanten B-Splines ohnehin der Schnitt zwischen den Gitterzellen und dem Rand Γ zu berechnen ist, ist somit zur Klassifizierung kein wesentlicher zusätzlicher Rechenaufwand notwendig. Der Gewichtspunkt x_i wird als Mittelpunkt einer vollständig im Simulationsgebiet Ω liegenden Gitterzelle im Träger des B-Splines b_i gewählt. Auch dies ist effizient möglich, da die Bestimmung einer solchen Zelle bereits Bestandteil der Klassifizierungsroutine ist.

Falls keine explizite analytische Form der Gewichtsfunktion bekannt ist, wird diese wie folgt angesetzt:

Fig. 12 illustriert die Konstruktion der Gewichtsfunktion. Dabei gibt der Parameter δ die Breite eines Streifens Ω_δ an, innerhalb dessen die Gewichtsfunktion zwischen dem Wert 0 auf dem Rand des Simulationsgebietes und dem Wert 1 auf einem Plateau über Ω/Ω_δ variiert. Der Parameter δ wird so bestimmt, dass die Glattheit der Gewichtsfunktion gewährleistet ist.

Ein wesentlicher Vorteil des Verfahrens ist, dass keine Vernetzung des Simulationsgebietes erforderlich ist. Dies bewirkt in technischen Anwendungen eine deutliche Ersparnis an Rechenzeit und Speicherkapazität und vereinfacht den Simulationsablauf. Dabei ist die Verfahrensstruktur für zwei- und dreidimensionale Probleme formal und auch in der technischen Durchführung weitgehend identisch. Dies ermöglicht die zeit- und kostensparende Implementierung von Lösern für verschiedenste Anwendungen basierend auf einheitlichen Programmstrukturen. Die Verwendung von B-Splines entspricht dem industriellen Standard bei der Modellierung geometrischer Objekte und schafft damit eine natürliche Verbindung zwischen FE- und CAD/CAM-Anwendungen. Bestehende umfangreiche Programmbibliotheken aus beiden Bereichen können zur Implementierung einer auf dem erfindungsgemäßen Verfahren beruhender FE-Simulation eingesetzt werden. Die nach dem WEB-Verfahren konstruierten Basisfunktionen besitzen alle Standardeigenschaften Finiter Elemente. Hierzu gehört insbesondere die Stabilität der Basis. Diese bewirkt, dass beispielsweise bei linearen elliptischen Randwertproblemen die Konditionszahl des entstehenden Gleichungssystems bei kleiner werdender Gitterweite nicht schneller anwächst als bei optimalen Triangulierungen. Für die Anwendungen bedeutet dies beispielsweise, dass sich lineare Gleichungssysteme, wie sie bei der FE-Methode typischerweise entstehen, mittels iterativer Verfahren effizient lösen lassen. Weiterhin ist bei gegebenem Grad die Approximationsordnung maximal und die Anzahl der benötigten Parameter minimal. Damit sind sehr genaue Approximationen schon mit relativ geringer Parameterzahl möglich. Konkret kann dies bedeuten, dass Genauigkeiten, die bislang den Einsatz von Großrechnern erforderten, auf Arbeitsplatzrechnern erzielt werden können. Die reguläre Gitterstruktur der erfindungsgemäßen Basis ermöglicht eine sehr effiziente Implementierung, insbesondere bei der Aufstellung und Lösung von FE-Systemen. Ferner gelingt es durch Verwendung der Gewichtsfunktion Randbedingungen bei der Simulation einzuhalten, ohne die reguläre Gitterstruktur der Basisfunktionen zu beeinträchtigen. Schließlich kann bei der Lösung linearer elliptischer Randwertprobleme durch den Einsatz von Mehrgitter-Verfahren zur Lösung des entstehenden linearen Gleichungssystems erreicht werden, dass die Gesamtlösungszeit proportional zur Anzahl der gesuchten Koeffizienten und damit optimal ist.

Das erfindungsgemäße Verfahren in der speziellen bevorzugten Ausführungsart (WEB-Verfahren) wird anhand des in Fig. 5 gezeigten ersten Ausführungsbeispiels veranschaulicht. Die Differenzialgleichung und Randbedingungen sind sehr elementar gewählt, so dass ohne großen Mehraufwand neben der erfindungsgemäßen Konstruktion der WEB-Basis der gesamte Ablauf der FE-Simulation verfolgt werden kann.

In der oberen Bildhälfte ist eine elastische Membran gegeben, die entlang dem Rand Γ eines ebenen Simulationsgebietes Ω fixiert ist und auf die im Inneren des Gebietes ein konstanter Druck f = 1 wirkt. Die Verschiebung u erfüllt bei geeigneter Normierung die Poisson-Gleichung mit homogenen Randbedingungen,

-Δu = 1 in Ω
u = 0 auf Γ.

Die Verschiebung u bzw. Auslenkung der Membran ist in der unteren Bildhälfte der Fig. 5 ersichtlich. Wie oben beschrieben, gliedert sich das WEB-Verfahren in die folgenden Schritte.

Eingabe 1 des Simulationsgebietes Ω: Der Rand Γ ist eine periodische Spline-Kurve vom Grad 6, die durch ihre Kontrollpunkte 20 (in Fig. 5 mit schwarzen Punkten gekennzeichnet) gespeichert wird.

Eingabe 2 der Randbedingungen: Die homogene Randbedingung ist wesentlich, so dass die Konstruktion einer Gewichtsfunktion erforderlich ist.

Eingabe 3 der Steuerparameter: Es wird der Grad n = 2 und, um die Figuren übersichtlich zu gestalten, eine relativ große Gitterweite h = 1/3 verwendet.

Bestimmung 4 der Zelltypen: Wie in Fig. 6 illustriert, wird das Simulationsgebiet durch ein Gitter 21 überdeckt, das die Gitterzellen der Träger aller potenziell für die Basiskonstruktion relevanten B-Splines enthält. Die Typbestimmung ergibt in dem Beispiel 69 äußere Gitterzellen 22 und 11 innere Gitterzellen 24, sowie 20 Gitterzellen 23 auf dem Rand.

Klassifikation 5 der B-Splines: Der Träger des B-Splines b_k ist hier das Quadrat Q_k mit den Ecken

(k₁, k₂)h, (k₁ + 3, k₂)h, (k₁ + 3, k₂ + 3)h, (k₁, k₂ + 3)h

Q_(-4,0) und Q_(2,1) sind in Fig. 8 gezeigt. Die Gitterpunkte kh der relevanten B-Splines, für die Q_k das Innere des Simulationsgebietes schneidet, sind in Fig. 8 durch einen Punkt oder einen Kreis markiert. Alle Gitterpunkte ih für innere B-Splines (i aus der Indexliste I), für die mindestens eine Zelle des Trägers Q_i vollständig im Inneren von Ω liegt, sind durch einen Punkt markiert. Beispielsweise liegt für i = (-4, 0) die Gitterzelle (-2, 0)h + [0, h]² ganz in Ω. Alle Gitterpunkte jh für äußere B-Splines (j aus der Indexliste J), für die keine Zelle des Trägers Q_j ganz in Ω liegt, sind durch einen Kreis markiert.

Berechnung 6 der Kopplungskoeffizienten: Zur Bestimmung der Kopplungskoeffizienten e_i,j wird für jedes feste j der Indexliste J der nächstgelegene 3 × 3-Bereich

I(j) = {l₁, l₁ + 1, l₁ + 2} × {l₂, l₂ + 1, l₂ + 2}

von Indizes in I gesucht. In Fig. 10 ist für den mit einem Kreis markierten äußeren Gitterpunkt j = (-1, 2) der Bereich I(j) durch Punkte gekennzeichnet. Die Fig. 10 zeigt ebenfalls die zugehörigen Kopplungskoeffizienten in Matrixdarstellung. Sie werden mit bivariater Interpolation berechnet. Beispielsweise ist das Interpolationspolynom zu i = (-1, -1)

p_i(x) = -x₁(x₁ + 2)x₂(x₂ - 1)/2.

Sein Wert an der Stelle x = j = (-1, 2) ist p_i(j) = 1 = e_i,j. Es fällt auf, dass viele der Kopplungskoeffizienten 0 sind. Dies ist ein typisches Phänomen. Die Kopplungskoeffizienten e_i,j sind nur dann für alle i der Indexliste I(j) ungleich 0, wenn die Indizes i in jeder Komponente von dem Index j verschieden sind.

Berechnungsvorschrift 7 für die Gewichtsfunktion: Die Gewichtsfunktion ist durch die Gleichung (3) gegeben mit n = 2 und δ = 0.2. Der Parameter δ wird numerisch berechnet. Er muss klein genug sein, so dass Singularitäten der Abstandsfunktion vermieden werden. Zur Berechnung der Abstandsfunktion erstellt das Verfahren ein Programm, das mit Hilfe des Newton-Verfahrens arbeitet. Da die Gewichtsfunktion nur in einem Randstreifen ungleich 1 ist, ist der Aufwand bei der späteren Auswertung gering.

Ausgabe

Fig. 13 zeigt den Träger eines WEB-Splines B_i und die zu seiner Beschreibung benötigten Daten. Dies sind die Indexliste J(i) der mit b_i gekoppelten äußeren B-Splines b_j, die Kopplungskoeffizienten e_i,j und der Gewichtspunkt x_i. Diese Daten werden unter Einbeziehung der Gewichtsfunktion zur Erstellung einer Berechnungsvorschrift für die WEB-Splines benutzt.

Der weitere Ablauf der FE-Simulation folgt dem Stand der Technik.

Aufstellen 9 des FE-Systems: Die Einträge der Systemmatrix und der rechten Seite sind

G_k,i = ∫_Ωgrad B_k grad B_i, F_k = ∫_Ωf B_k, k, i ∈ I.

Das Gleichungssystem GC = F für die Basiskoeffizienten C_i hat im betrachteten Beispiel die Dimension 31. Die Matrixeinträge G_k,i werden mit Hilfe von numerischer Integration berechnet, ebenso die Integrale F_k.

Lösen 10 des FE-Systems: Das Galerkin-System wird iterativ mit dem Verfahren der konjugierten Gradienten gelöst, wobei die SSOR-Vorkonditionierung zur Konvergenzbeschleunigung Verwendung findet. Nach 24 Iterationsschritten wird die Lösung im Rahmen der Rechnergenauigkeit gefunden (Toleranz ≦ 1e - 14).

Berechnung 11 und Ausgabe 12 der Approximation: Die mit dem erfindungsgemäßen Verfahren berechnete Approximation ist u = Σ_iC_iB_i und wird in in der unteren Bildhälfte von Fig. 5 grafisch dargestellt. Der relative Fehler in der L₂-Norm beträgt 0.028.

Die Leistungsfähigkeit des erfindungsgemäßen Verfahrens in der speziellen bevorzugten Ausführungsart (WEB-Verfahren) wird in einem zweiten Ausführungsbeispiel anhand der Simulation einer inkompressiblen Strömung veranschaulicht. Die in Fig. 14 gezeigte Anordnung zweier kreisförmiger Hindernisse in einem parallel berandeten Kanal dient der Illustration des prinzipiellen Vorgehens. Bei komplizierteren Geometrien, wie sie in konkreten Anwendungen typischerweise vorliegen, arbeitet das Verfahren vollkommen analog und effizient. In Fig. 14 sind die Strömungslinien 25 innerhalb des durch Γ₁ bis Γ₄ sowie durch Γ₅ und Γ₆ berandeten Gebiets dargestellt. Die Differenzialgleichung ist

Δu = 0 in Ω

mit den Randbedingungen

Die Flussgeschwindigkeit v = -gradu ist in der unteren Bildhälfte dargestellt.

Die erfindungsgemäße Konstruktion der WEB-Basis läuft vollkommen analog zu dem ersten Ausführungsbeispiel ab. Der einzige Unterschied ist, dass aufgrund der natürlichen Randbedingungen keine Gewichtsfunktion erforderlich ist.

Fig. 15a bis Fig. 15c zeigen die Klassifikation der relevanten B-Splines für verschiedene Grade n (s. auch Fig. 8). In der Abbildung sind die inneren B-Splines b_i, die in die WEB-Basis ohne Erweiterung übernommen werden, durch schwarze Dreiecke markiert. Für kleines h nimmt die Anzahl dieser B-Splines zu, d. h., es gilt B_i = b_i für den überwiegenden Teil der WEB-Basis. Im betrachteten Beispiel trifft dies bei Grad n = 3 für 236 von 252 Indizes i zu.

Fig. 16a zeigt in zwei Diagrammen den numerisch bestimmten relativen L₂-Fehler des Potenzials (linke Bildhälfte) in Abhängigkeit von der Gitterweite h = 2^-k mit k = 1, . . ., 5 und die numerisch geschätzte Konvergenzordnung m (rechte Bildhälfte). Dabei werden für die verschiedenen Grade der WEB-Splines als Kennzeichner * (n = 1), ○ (n = 2), ∆ (n = 3), (n = 4) und (n = 5) verwendet. Erwartungsgemäß ergibt sich m ≈ n + 1, d. h. eine ungefähre Fehlerreduktion um den Faktor 2ⁿ⁺¹ bei Halbierung der Gitterweite. Analog ergibt sich für den in der Fig. 16b abgebildeten relativen Approximationsfehler der Flussgeschwindigkeit (H¹-Norm der Lösung, linke Bildhälfte) eine Konvergenzordnung von m ≈ n (rechte Bildhälfte) und die damit verbundene Fehlerreduktion um etwa den Faktor 2ⁿ bei Halbierung der Gitterweite.

Fig. 16c (rechte Bildhälfte) zeigt die Rechenzeit in Sekunden für die Konstruktion der WEB-Basis in Abhängigkeit von der Zahl der resultierenden Basisfunktionen gemessen auf einem Pentium II-Prozessor mit 400 MHz. Beispielsweise werden für die Konstruktion einer WEB-Basis vom Grad 3 zur Gitterweite h = 0.125 mit 2726 WEB-Splines 1.32 Sekunden benötigt. Es fällt auf, dass der Aufwand zum Erstellen der WEB-Basis weitgehend unabhängig vom Grad n der Basis ist. In der linken Bildhälfte der Fig. 16c ist die Anzahl der CG-Iterationen bezogen auf die Anzahl der Basisfunktionen gezeigt. So werden für das zugehörige System mit 2726 Unbekannten 65 PCG-Iterationen benötigt. Die Gesamtrechenzeit inklusive Aufstellen und Lösen des Galerkin-System beträgt etwa 2.48 Sekunden.

In Fig. 17 wird die WEB-Methode mit einem Standardlösungsverfahren, das das Simulationsgebiet vernetzt bzw. trianguliert (obere Bildhälfte) und Hutfunktionen verwendet, verglichen. Die Grafik zeigt in der unteren Bildhälfte den L₂-Fehler bezogen auf die Anzahl der Parameter. Die Ergebnisse des Standardlösers sind mit fetten Karo-Kennzeichnern markiert und werden den mit Hilfe von WEB-Basen der Grade 1 bis 5 erzielten Ergebnissen gegenübergestellt. Beispielsweise wird eine Genauigkeit von 10^-2 mit der WEB-Methode durch Verwendung von 213 Basisfunktionen mit Grad 2 erreicht bei einer Gesamtrechenzeit von 0.6 Sekunden. Zum Erzielen der gleichen Genauigkeit benötigt das Standardverfahren mit linearen Hutfunktionen 6657 Basisfunktionen.

Bei der Beurteilung des Standardlösungsverfahrens sind noch zwei weitere Aspekte zu beachten. Zum einen verdeutlicht Fig. 17, dass selbst eine moderate Genauigkeit von 10^-3 mittels Hutfunktionen nur erreichbar ist, wenn weit mehr als eine Million Koeffizienten verwendet werden. Dies zeigt, dass genaue Resultate bei der Verwendung von Hutfunktionen in der Regel eine enorme Rechen- und Speicherkapazität erfordern oder aber nach dem Stand der Technik gar nicht erreichbar sind. Zum anderen steigt der zur Vernetzung erforderliche Aufwand mit der Komplexität des Simulationsgebietes. Im Gegensatz zu realistischen Anwendungen ist das hier untersuchte Gebiet aufgrund seiner sehr einfachen Struktur noch vergleichsweise einfach zu triangulieren.

Das zweidimensionale Beispiel zeigt den Performancegewinn durch das WEB-Verfahren. Eine noch größere Leistungssteigerung ist bei dreidimensionalen Problemen möglich. Zum einen ist dort der Aufwand für die Vernetzung, die bei dem WEB-Verfahren entfällt, wesentlich größer. Zum anderen macht sich hier auch die Reduktion der Anzahl der benötigten Basisfunktionen noch wesentlich deutlicher bemerkbar als im zweidimensionalen Fall.

Die Fig. 18 zeigt eine erfindungsgemäße Vorrichtung, insbesondere eine Computereinrichtung 30, mit Eingabemitteln 31, 32, 33, Ausgabemitteln 34 sowie eine den Verfahrensablauf steuernde Steuereinrichtung 35. Zur Abarbeitung des erfindungsgemäßen Verfahrens und insbesondere zum Zwecke der Parallelisierung des zugehörigen Rechenablaufs bedient sich die zentrale Steuereinrichtung 35 vorzugsweise mehrerer arithmetischer logischer Einheiten (Arithmetic Logical Unit, ALU) oder sogar mehrerer zentraler Verarbeitungseinheiten (Central Processing Unit, CPU) 36. Diese erlauben insbesondere eine parallele Abarbeitung der Verfahrensschritte Klassifizieren 5 der B-Splines, insbesondere auch Schneiden des regulären Gitters mit dem Simulationsgebiet Ω, Ermitteln 6 der Kopplungskoeffizienten e_i,j und/oder Auswerten der Gewichtsfunktion w(x) an Punkten x des Simulationsgebietes Ω.

Die Recheneinheiten 36 greifen dabei auf den gemeinsamen Datenbestand von Speichermitteln 37 zu. Die Eingabe der Daten kann beispielsweise mittels einer Tastatur 31, einem maschinenlesbaren Datenträger 38 über eine zugehörige Lesestation 32 und/oder über ein drahtgebundenes oder drahtloses Datennetz mit einer Empfängerstation 33 erfolgen. über die Lesestation 32 bzw. einen zugehörigen Datenträger 38 kann auch das den Verfahrensablauf steuernde Steuerprogramm eingegeben werden und beispielsweise auf dem Speichermittel 37 dauerhaft abgelegt werden. Entsprechend kann das Ausgabemittel 34 einen Drucker, einen Monitor, eine Schreibstation für einen maschinenlesbaren Datenträger und/oder eine Sendestation eines drahtgebundenen oder drahtlosen Datennetzes aufweisen.

Claims (12)

1. Verfahren zur Erhöhung der Leistungsfähigkeit einer Computereinrichtung zur Finiten Elemente-Simulation durch automatische Generierung geeigneter Basisfunktionen mit Hilfe von uniformen Tensorprodukt B-Splines mit den folgenden Schritten:

• - Definieren (1) eines Simulationsgebietes (Ω) und Abspeichern von Gebietsdaten des Simulationsgebietes (Ω),
• - Eingabe (2) und Abspeichern von Randbedingungen,
• - Festlegen (3) einer vorgebbaren Gitterweite (h) und eines vorgebbaren Grades (n) der B-Splines,
• - Ermitteln (4) eines sich über das Simulationsgebiet (Ω) erstreckenden Gitters und eines Typs der Gitterzellen,
• - Klassifizieren (5) der B-Splines, deren Träger das Simulationsgebiet (Ω) schneiden, durch Ermitteln von inneren und äußeren B-Splines, wobei für innere B-Splines der Schnitt des Trägers mit dem Simulationsgebiet (Ω) ein Quadrat bzw. einen Würfel einer vorgebbaren Kantenlänge (s), welche mit der Weite des Gitters (h) übereinstimmen kann, enthält,
• - Ermitteln (6) von Kopplungskoeffizienten (e_i,j) zur Bildung von Linearkombinationen von inneren und äußeren B-Splines durch Lagrange-Interpolation auf zu den äußeren Indizes nächstliegenden Feldern innerer Indizes,
• - Speicherung und Ausgabe der die Basisfunktionen bestimmenden Parameter.

2. Verfahren nach Anspruch 1, dadurch gekennzeichnet, dass vor dem Speichern und Ausgeben der Parameter der Schritt durchgeführt wird:
Festlegen (7) einer vorgebbaren Gewichtsfunktion (w) zur Einhaltung wesentlicher Randbedingungen sowie Bestimmen (8) von Gewichtspunkten und Skalierungsfaktoren.

3. Verfahren nach Anspruch 2, dadurch gekennzeichnet, dass die Gewichtsfunktion w durch einen glatten Übergang von einem konstanten Plateau im Inneren des Simulationsgebietes (Ω) zum Wert Null auf dem Rand (Γ) festgelegt wird.

4. Verfahren nach einem der Ansprüche 1 bis 3, dadurch gekennzeichnet, dass B-Splines mit mindestens einer ganz im Simulationsgebiet (Ω) enthaltenen Gitterzelle des Trägers als innere B-Splines klassifiziert werden.

5. Verfahren nach einem der Ansprüche 1 bis 4, dadurch gekennzeichnet, dass der Gewichtspunkt als Mittelpunkt einer ganz im Simulationsgebiet (Ω) enthaltenen Gitterzelle des Trägers des entsprechenden inneren B-Splines gewählt wird.

6. Verfahren nach einem der Ansprüche 1 bis 5, dadurch gekennzeichnet, dass das Simulationsgebiet (Ω) durch Einspeichern von aus einer computergestützten Konstruktion (CAD/CAM) ableitbaren Daten definiert wird.

7. Verfahren nach einem der Ansprüche 1 bis 6, dadurch gekennzeichnet, dass die Gitterweite h unter Verwendung von abgespeicherten empirisch und/oder analytisch gewonnenen Werten mittels einer zugehörigen ersten Bewertungsfunktion automatisch festgelegt wird.

8. Verfahren nach einem der Ansprüche 1 bis 7, dadurch gekennzeichnet, dass der Grad n unter Verwendung von abgespeicherten empirisch und/oder analytisch gewonnenen Werten mittels einer zugehörigen zweiten Bewertungsfunktion automatisch festgelegt wird.

9. Verfahren nach einem der Ansprüche 1 bis 8, gekennzeichnet durch die Schritte:

• - Aufstellen (9) eines bei der FE-Simulation zu lösenden Gleichungssystems,
• - Lösen (10) des Gleichungssystems,
• - Berechnen (11) einer Näherungslösung, und
• - Ausgeben (12) der Näherungslösung.

10. Verfahren nach Anspruch 9, dadurch gekennzeichnet, dass beim Lösen (10) des Gleichungssystems ein Mehrgitter-Verfahren eingesetzt wird.

11. Vorrichtung zur Durchführung eines Verfahrens nach einem der Ansprüche 1 bis 10, insbesondere Computereinrichtung, mit Eingabemitteln (31, 32, 33) und Ausgabemitteln (34), Speichermitteln (37) und einer Recheneinheit (35, 36), dadurch gekennzeichnet, dass die reguläre Gitterstruktur zur Optimierung des Rechenablaufs, insbesondere durch Parallelisierung, ausgenutzt wird.

12. Maschinenlesbarer Datenträger (38), insbesondere Magnetband, Magnetscheibe, Compact Disc (CD) oder Digital Versatile Disc (DVD), dadurch gekennzeichnet, dass auf dem Datenträger (38) ein Steuerprogramm für eine Computereinrichtung (30) gespeichert ist, gemäß dem durch die Computereinrichtung (30) ein Verfahren nach einem der Ansprüche 1 bis 10 abgearbeitet werden kann.