Bestimmung des effektiven Elastizitätsmoduls von Mikrostrukturen Ingenieurwissenschaften

RVE_stochastische Strukturen_3

Die meisten Werkstoffe haben eine Mikrostruktur, die man erst unter dem Mikroskop erkennt. Die Verformbarkeit eines Bauteils wird maßgeblich durch diese Struktur bestimmt. In dem Fachgebiet der Mikromechanik geht man der Frage nach, wie das makroskopische Verhalten über die Kenntnis der Mikrostruktur vorhergesagt werden kann.
Vereinfacht gesagt geht es um die Frage, wie man das E-Modul des Hookschen Gesetzes, das die Proportion von Spannung zu Dehnung beschreibt, auf Basis der Werkstoffstruktur berechnen kann. Dazu muss man wissen, dass Bauteile ingenieurstechnisch mit der kontinuumsmechanischen Theorie berechnet werden und dort das Materialverhalten vollkommen glatt (in der Fachsprache homogenisiert) angenommen wird. Das entspricht einer Mittelwertbildung der eigentlich streuenden Verformungen im Werkstoff.

Ich habe für einen Würfel einmal die erste Hauptnormaldehnung visualisiert und einen Schnitt animiert. Gedanklich muss man sich vorstellen, dass die homogenisierte Dehnung  im ganzen Würfel konstant einen mittleren Wert annimmt, statt wie dargestellt von lokalen Gegebenheiten abzuhängen.

Russ_LE

Schnitt durch einen Würfel und Darstellung der ersten Hauptnormaldehnung. Schwarze Bereiche entsprechen den Agglomeraten.

Im einfachsten Fall verhält sich der Werkstoff in allen Verformungsrichtungen des Kontinuums gleich, was man als isotropes Verhalten bezeichnet. Diese Verhalten kann durch zwei Module, bzw. Materialparameter beschrieben werden. Das ist allerdings eher der Ausnahmefall. Z.B. ist bei gewalztem Stahl oder eingebetteten Kohlenstofffasern das Verformungsverhalten in die entsprechenden Richtungen steifer. Hier spricht man von einem anisotropen Verhalten, das durch eine unabhängige Steifigkeitsmatrix statt durch isotrope Module modelliert werden kann, dessen Einträge entsprechend den Raumrichtungen und den Wechselwirkungen zwischen den Raumrichtungen bestimmt werden müssen. Konsequent müsste man die Steifigkeitsmatrix eigentlich Modulmatrix nennen und den Begriff Steifigkeit nur auf makroskopischer Ebene verwenden. In den Ingenieurswissenschaften herrscht hier ein sprachlicher Mischmasch. Es gibt die Steifigkeit als Werkstoffeigenschaft und als Bauteileigenschaft. Da man die Werkstoffeigenschaft nur an Bauteilen messen kann (denn jeder Probekörper ist schließlich ein Bauteil), bleibt die Werkstoffeigenschaft des materiellen Punktes immer nur ein Modell, das nur über Geometrie und spezielle Lastfälle validiert werden kann. Um das sprachlich nicht durcheinanderzubringen, sollte man die Eigenschaft des Werkstoffs anders benennen, z.b. Modul, Modulmatrix, Parameter des Materialgesetzes, etc. Die Steifigkeit wäre dann das qualitative Pendant auf der Bauteilebene, die man wiederum über Moden quantifizieren kann. Zu allem Überfluss brennt sich während der Ingenieursausbildung der Begriff Steifigkeit als Synonym für die Verformbarkeit einer Schraubfeder ein. Damit hat man allerdings Schwierigkeiten die Steifigkeit zwischen zwei materiellen Punkten im Kontinuum zu verstehen, denn zwischen diesen denkt man sich dann eben auch intuitiv eine Schraubfeder, was falsch ist. Aber ich schweife ab. Also zurück zum Thema…

Es ist alles andere als trivial, beliebige Mikrostrukturen ingenieurstechnisch für Bauteile zu berücksichtigen. Die erforderliche Homogenisierung, bzw. Mittelwertbildung, funktioniert zunächst einmal generell nicht für alle Mikrostrukturen, da eine klassische Formulierung des Kontinuums kein Momenten-Gleichgewicht am infinitesimalen Volumenelement berücksichtigt, von Nichtlinearitäten ganz zu schweigen.

Nun sollte man allerdings differenzieren zwischen dem Anspruch an eine möglichst genaue Modellierung innerhalb der Materialwissenschaft und dem ingenieurstechnischen Bedarf an vernünftigen Materialparametern für praktische Berechnungen.
Für eine Bauteilauslegung brauche ich nicht immer ein perfektes Materialmodell. Z.B. interessiert mich häufig nur die Verformungscharakteristik einer komplexen Geometrie, für die ich kein empirisches Gefühl habe, um Schwachstellen im Spannungsfluss zu identifizieren. Oder es interessiert mich, ob die ersten Eigenschwingmoden im 1-10 Hz oder 10-100 Hz Bereich liegen. Es ist ein Irrtum zu glauben, dass Ingenieure Berechnungen nutzen, um die Realität detailgetreu abzubilden. Ingenieure wollen am Ende des Tages ein funktionierendes Bauteil in der Hand haben, egal wie sie es geschafft haben, um das einmal salopp zu sagen. Bei kritischen und sicherheitsrelevanten Bauteilen werden eh experimentelle Tests am Bauteil durchgeführt. Wenn mal ehrlich ist, werden aufwändige Berechnungen häufig nur durchgeführt, wenn so ein Test nicht bestanden wird. Durch die Berechnung identifiziert man dann die Schwachstelle und kann sie nun mit einfachen Mitteln am Rechner optimieren. Die Methoden, mit denen Bauteile ausgelegt werden, sind wesentlich vielfältiger als durch alleinige Berechnung mit einem FEM-Programm. Und das ist natürlich absolut gut und richtig so. Auf der anderen Seite gibt es die standardisierten Sicherheitsnachweise, doch die führen in der Regel zu vollkommen überdimensionierten Bauteilen. Die Auslegung scheitert meist nicht an der Berechnung, sondern in der Regel an den Lastkollektiven, die man für die Rechnung verwendet. Doch auch hier mache ich ein anderes Fass auf. Zurück zum Thema…

Ich frage mich als Ingenieur, wie ich eine optimale Näherung in Form von isotropen Modulen aus einer beliebigen Mikrostruktur ermitteln kann. Es muss doch möglich sein, die Mikrostruktur selbst in einem FEM-Programm zu berechnen, daraus die homogenisierten Module zu gewinnen, die ich dann ich einer Bauteilberechnung verwenden kann. Warum sollte ich in der heutigen Zeit auf analytische Näherungslösungen der Mikromechanik zurückgreifen? Die habe ich nachgerechnet. Es sind so grobe Näherungen, zumindest in dem folgenden Beispiel, das sie in der Praxis nicht zu gebrauchen sind.

Vergleich mit analytischen Lösungen

Vergleich der FEM-Lösung mit analytischen Verfahren der Mikromechanik.

Die Idee ist also recht simpel. Ich modelliere die Mikrostruktur in einem Würfel und mache einen virtuellen Zugversuch in einem FEM-Programm. Über Kraft, Oberfläche und Verformung berechne ich dann das gemittelte, effektive Modul. Das kann man zwar genau so machen, nur kontinuumsmechanisch korrekt ist das nicht. Man gewichtet auf diese Weise einen speziellen Lastfall, den Uniaxialen. Wie gesagt, ein Modul entspricht nicht dem Steifigkeitswert einer Schraubfeder. Mein Anspruch ist deshalb, die Module, bzw. die Steifigkeitsmatrix für alle Verformungen im Kontinuum korrekt zu ermitteln. Ich möchte außerdem nicht nur den effektiven Wert, sondern auch den potentiellen Fehler bestimmen können, den ich durch die isotrope Homogenisierung eines eigentlich anisotropen Verhaltens vornehme.

Dieser Fehler ist am einfachsten zu kontrollieren, wenn sich eine anisotrope Mikrostruktur auf makroskopischer Ebene durch stochastische Verteilung quasi isotrop verhält, wie z.b. ein partikelverstärkter Kunststoff.

Zunächst habe ich zusammen mit einer Kollegin, die Berechnungsstrategien für Knochengewebe entwickelt, um Implantverankerungen zu optimieren, die also auch an der Homogenisierungsmethode interessiert war, einen Algorithmus entwickelt, um dreidimensionale Strukturen auf stochastischer Basis zu erstellen. Dazu werden in dem Würfel zufällige Startzellen ausgewählt. Anschließend werden weitere Zellen zufällig ausgewählt, allerdings verworfen, wenn sie nicht direkt zu bestehenden benachbart sind.

Algorithmus zur Zellenerzeugung

Per Zufall werden Materialeigenschaften zugeordnet, aber nur an benachbarten Zellen, bis der gewünschte Volumenanteil erreicht ist

So entstehen langsam abhängig von Anfangszahl zufällige Agglomeratstrukturen, bis der gewünschte Füllgrad erreicht ist. Rußstrukturen im Elastomer kann man mit dieser einfachen Regel erstaunlich realistisch nachbilden. Knochengewebe funktioniert auf den ersten Blick eher nicht so gut, die Ergebnisse erweisen sich aber am Ende ebenfalls als absolut brauchbar.
An solchen Würfeln berechneten wir nun jeweils sechs verschiedene Lastfälle mit der FEM. Die Art der Randbedingungen und die Würfelgröße haben einen großen Einfluss auf die Güte des Ergebnisses, doch ich will die vielfältigen numerischen Komplikationen hier nicht vertiefen. Bei Interesse können Sie die wissenschaftlichen Details gerne den verlinkten Veröffentlichungen am Ende des Artikels entnehmen.

Für alle Lastfälle werden die 6 Spannungen und 6 Dehnungen aller Elemente ausgewertet und jeweils gemittelt. Für jeden Lastfall lässt sich dann eine Spalte der 6×6 Steifigkeitsmatrix berechnen.

Bestimmungs der Einträge in die Steifigkeitsmatrix

Nach sechs Lastfällen ist sie vollständig ermittelt, aber mit 36 Werten überbestimmt. Da die Steifigkeitsmatrix symmetrisch sein muss, mitteln wir die symmetrischen Einträge und es verbleiben 21 Module. Das inhomogene Gefüge ist nun aus kontinuumsmechanischer Sicht korrekt in ein anisotropes Gefüge homogenisiert. Ein gleichmäßiger, bzw. homogener, Würfel mit der ermittelten, anisotropen Steifigkeitsmatrix verformt sich im Mittel genau so, wie die inhomogene Mikrostruktur.
Da die Struktur stochastisch variiert, erzeugten wir über eine Monte-Carlo-Simulation so lange verschiedene Strukturen, bis der Mittelwert für die jeweiligen statistischen Parameter auskonvergiert ist. So weit nichts Neues aus wissenschaftlicher Sicht. Unser Großrechner glühte voll automatisiert einige Wochen lang und erzeugte einige Terrabyte an Daten.

Im letzten Schritt bestand nun die Herausforderung, die beiden effektiven, isotropen Module aus der anisotropen Steifigkeitsmatrix zu ermitteln.

anisotrop2isotrop

Bestimmte Nebendiagonalen müssen dazu verschwinden. Da sich die verbleibenden Einträge gegenseitig bedingen, gibt es unendlich viele Möglichkeiten, um isotrope Module zu bestimmen. Hier haben wir dann Neuland betreten. Wissenschaftlich ist es meist wesentlich einfacher, den Grenzbereich, die sogenannten Schranken, eines gesuchten Wertes zu bestimmen, als den Wert an sich. Diese Schranken haben wir, analog zu den analytischen Schranken, mit mikromechanischen Methoden bestimmt. Im Idealfall haben beiden Schranken den gleichen Wert. Es zeigte sich, dass sie tatsächlich mit zunehmender Würfelgröße zusammenfallen.

Konvergenzverhalten_

Monte-Carlo-Simulation zur Bestimmung von konvergierten Lösungen.

Schließlich konnten wir eine Würfelgröße mit 50x50x50 Elementen identifizieren, bei der wir guten Gewissens den Mittelwert aus beiden Schranken bilden konnten. Das waren nun die gewünschten effektiven Module. Da die makroskopische Struktur vom Prinzip isotrop sein müsste, konnten wir auch noch den Anisotropiefehler unserer Homogenisierung bestimmen. Er lag unter 5%. Für makroskopisch anisotrope Werkstoffe würde dieser Wert natürlich nicht gegen null konvergieren. Als Maß, wie fehlerhaft die isotrope Homogenisierung ist, funktioniert es jedoch sehr gut.

Da wir nur die Anfangszellen und damit die Strukturbildung und den Füllgrad variiert haben, lassen sich die Ergebnisse in einfachen Diagrammen veranschaulichen.

Einfluss des Volumenanteils_

E-Modul in Abhängigkeit des Füllgrads und der Struktur mit empirischer Ausgleichsformel

Einfluss der Struktur

E-Modul-Optimum für bestimmte, stochastische Strukturen

Hinter jeden Datenpunkt stecken 48 Rechnungen, die jeweils 4 Stunden dauerten. Mit ein bisschen Knobelei konnten wir eine empirische Formel finden, die eine Ausgleichkurve durch alle FEM-Einzelergebnisse legt. Mit dieser Formel kann man nun z.b. für Elastomere vorherbestimmen, welcher Rußtyp und welcher Volumenanteil zu welcher Modulsteigerung führt.

Die Ergebnisse decken sich übrigens verblüffend gut mit experimentellen Versuchen, erstaunlicherweise auch für Knochen.

Vergleich mit Messungen_

Vergleich mit experimentellen Messwerten

Nähere Infos hier: EnglischDeutsch

Bestimmung des effektiven Elastizitätsmoduls von Mikrostrukturen
5 (100%) 1 vote

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht. Erforderliche Felder sind mit * markiert.