Klassische Plastizität und ihre thermische Erweiterung ... · 1 Einführung 1 Einführung Die...

Technische Universität BerlinFakultät V – Verkehrs- und MaschinensystemeInstitut für MechanikLehrstuhl für Kontinuumsmechanik und Materialtheorie

Klassische Plastizität und ihre thermischeErweiterung mittels finiter Elementmethode

Zur Erlangung des akademischen Grades Bachelor of Science

Vorgelegt von: Andre Klunker, 325834

Berlin, der 21. Mai 2013

Gutachter: Prof. Dr. rer. nat. W. H. MüllerB. E. Abali, M. Sc.

Eidesstattliche Versicherung

Die selbständige und eigenhändige Anfertigung versichere ich an Eides statt.

Andre Klunker

Inhaltsverzeichnis

Abbildungsverzeichnis I

1 Einführung 1

2 Kontinuumstheoretische Grundlagen 2

2.1 Kinematik des Kontinuums . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2.1.1 Deformations- und Verzerrungsmaße . . . . . . . . . . . . . . . . . . . . . 3

2.1.2 Lagrangesche und Eulersche Darstellung der Bewegung . . . . . . . . 4

2.2 Kinematik kleiner Deformationen . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.3 Bilanzgleichungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.3.1 Massenbilanz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.3.2 Impulsbilanz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.3.3 Bilanz der inneren Energie . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.3.4 Entropiebilanz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3 Klassische Plastizität im mathematischen Modell 13

3.1 Die Prandtl-Reuss-Gleichungen . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.1.1 Vorbetrachtungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.1.2 Assoziierte Fließregel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3.1.3 Fließkriterium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3.1.4 Verfestigungsmodell – isotrope Verfestigung . . . . . . . . . . . . . . . . . 15

3.1.5 Grenzen des Modells – der Bauschinger-Effekt . . . . . . . . . . . . . . 19

3.2 Der Bauschinger-Effekt im Modell . . . . . . . . . . . . . . . . . . . . . . . . . 20

4 Analytische Berechnung von Spannungs-Dehnungs-Beziehungen 24

4.1 Modellproblem Scherklotz – Lösung der Impulsbilanz unter semi-inversem Ansatz 24

4.2 Modellproblem Scherklotz – Direkte Berechnung . . . . . . . . . . . . . . . . . . 29

4.3 Modellproblem Zugstab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

5 Thermische Erweiterung 33

5.1 Kopplung von innerer Energie und Entropie . . . . . . . . . . . . . . . . . . . . . 33

5.1.1 Grundgleichungen und Annahmen . . . . . . . . . . . . . . . . . . . . . . 33

5.1.2 Zustandsraumwechsel und Materialgleichungen . . . . . . . . . . . . . . . 34

5.2 Aufstellen der Feldgleichung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

6 Numerische Behandlung 39

6.1 Schwache Formulierung der Impulsbilanz . . . . . . . . . . . . . . . . . . . . . . . 39

6.2 Schwache Formulierung der Energiebilanz . . . . . . . . . . . . . . . . . . . . . . 40

6.3 Zeitdiskretisierung und Numerik der Materialgleichungen . . . . . . . . . . . . . 40

6.4 Gekoppelte variationelle Formulierung und Randbedingungen . . . . . . . . . . . 42

6.5 Prinzipielle Struktur der Berechnung . . . . . . . . . . . . . . . . . . . . . . . . . 43

7 Konkrete numerische Berechnungen 45

7.1 Vergleich zu den analytischen Lösungen – isotherme Berechnung mit isotropemVerfestigungsmodell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

7.2 Problematik bei kinematischer Verfestigung und thermischer Kopplung . . . . . . 46

7.3 Zugprobe unter wechselseitiger Belastung in zwei Dimensionen – Verfestigungs-modelle im Vergleich . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

7.4 Zugprobe unter wechselseitiger Belastung in zwei Dimensionen – Thermomecha-nische Kopplung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

7.5 Scherklotz unter wechselseitiger Belastung in zwei Dimensionen – Thermomecha-nische Kopplung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

8 Resümee 52

A Quelltexte 53

A.1 Scherklotz – 3D, isotherm, isotrope Verfestigung . . . . . . . . . . . . . . . . . . 53

A.2 Zugprobe – 3D, isotherm, isotrope Verfestigung . . . . . . . . . . . . . . . . . . . 55

A.3 Zugprobe – 2D, unter wechselseitiger Belastung – isotherm, isotrope Verfestigung 58

A.4 Zugprobe – 2D, unter wechselseitiger Belastung – isotherm, kinematische Verfes-tigung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

A.5 Zugprobe – 2D, unter wechselseitiger Belastung mit Temperaturentwicklung . . . 64

A.6 Scherklotz – 2D, unter wechselseitiger Belastung mit Temperaturentwicklung . . 68

Literaturverzeichnis 73

Abbildungsverzeichnis

2.1 Zur Bedeutung des Konfigurationsbegriffs und des Deformationsgradienten . . . 3

2.2 Zur Bedeutung des Verschiebungsvektors . . . . . . . . . . . . . . . . . . . . . . . 5

3.1 Zur Verfestigung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.2 Zur Bedeutung des plastischen Moduls h . . . . . . . . . . . . . . . . . . . . . . . 16

3.3 Zum Bauschinger-Effekt bei zyklischer Belastung . . . . . . . . . . . . . . . . . 19

4.1 Modellproblem Scherklotz – Verschiebungsrate als Randbedingung . . . . . . . . 24

4.2 Modellproblem Scherklotz – Scherspannungsrate als Randbedingung . . . . . . . 29

4.3 Modellproblem Zugstab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

7.1 Scherklotz in der numerischen Berechnung . . . . . . . . . . . . . . . . . . . . . . 45

7.2 Spannungs-Dehnungs-Diagramme – Analytische und numerische Lösungen. . . . 46

7.3 Zugprobe in der numerischen Berechnung . . . . . . . . . . . . . . . . . . . . . . 46

7.4 Dirichlet-Randbedingung am rechten Rand – ur1(t) . . . . . . . . . . . . . . . . 47

7.5 Rechengebiet mit verwendeten Randbedingungen . . . . . . . . . . . . . . . . . . 48

7.6 Zugprobe unter wechselseitiger Belastung, Fließen in beiden Belastungsrichtungen– Spannungs-Dehnungs-Diagramme für beide Verfestigungsmodelle . . . . . . . . 48

7.7 Spannungs-Dehnungs-Diagramm und Temperaturentwicklung der wechselseitigbelasteten Zugprobe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

7.8 Temperaturverteilung in der Probe zur Zeit t/tEnde = 0,2 . . . . . . . . . . . . . . 50

7.9 Spannungs-Dehnungs-Diagramm und Temperaturentwicklung des wechselseitigbelasteten Scherklotzes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

7.10 Ebener Scherklotz – numerische Lösungen für Verschiebungen und Temperatur . 51

I

1 Einführung

1 Einführung

Die Modellierung plastischer Verformungen, insbesondere von Metallen, findet in zahlreichenGebieten technische Anwendung. Sei es, um gezielt plastische Verformungen herbeizuführenwie bei der Blechumformung, beispielsweise an CNC-Biegepressen, oder um bei ungewollterplastischer Verformung Versagenshypothesen theoretisch zu stützen und das Ausmaß plastischerVerformungen in technischen Systemen abzuschätzen.

Obwohl die ersten experimentellen Erkenntnisse auf dem Gebiet der Plastizitätstheorie bereitsim 15. Jahrhundert gewonnen wurden, finden sich doch die Anfänge der kontinuumsmechani-schen Plastizitätstheorie, sprich der Theorie mit dem Ziel, ein mathematisches Materialmodellaufzustellen, erst zum Ende des 19. Jahrhunderts. Einige der ersten bekannten Aufzeichnungenüber Spannungs-Dehnungs-Zusammenhänge in inkrementeller Formulierung, die das Verhaltenvon Metallen über den Hookeschen linear elastischen Bereich hinaus beschreiben sollten, findensich in Frankreich. Als wegweisend gilt die Arbeit von de Saint-Venant (1870), der seiner-zeit zweidimensionale Probleme untersuchte, bei denen er jedoch noch keine Änderungen vonSpannungen nach Einsetzen des plastischen Fließens annahm. Lévy (1871) hat diese Theorieder sogenannten perfekten Plastizität auf den dreidimensionalen Fall erweitert. Dessen Arbeitwurde jedoch nicht außerhalb der Grenzen Frankreichs veröffentlicht und blieb dadurch langeZeit unbeachtet. So hat von Mises (1913) ohne Kenntnis von Lévys Arbeit dieselben Gesetz-mäßigkeiten postuliert. Auf dieser Theorie der perfekten Plastizität aufbauend haben zuerstPrandtl (1924) für den ebenen Fall und schließlich Reuss (1930) mit der dreidimensionalenVerallgemeinerung ein Modell präsentiert, das zusätzlich Verfestigungs- und Elastizitätseffektemiteinbezog und auf diese Weise sogar durch Lode (1926) experimentell verifiziert werden konn-te. Dieses Modell, welches als Prandtl-Reuss-Gleichungen in einem Großteil der Fachliteraturbehandelt wird, findet noch heute in modernen Simulationsprogrammen Anwendung.

Plastische Verformungen sind irreversible, dissipative Vorgänge. Somit sollen thermische Effek-te nicht ohne weiteres vernachlässigt werden. In dieser Arbeit soll zunächst die im damaligenKontext rein mechanisch betrachtete klassische Plastizitätstheorie behandelt werden, um die-se dann mit thermodynamischen Prinzipien zu erweitern und so die implizite Annahme derIsothermie aufzuheben. Auf Bilanzgleichungen aufbauend soll ein Gleichungssystem aufgestelltwerden, welches das Temperaturfeld als gesuchte Größe miteinbezieht. Schließlich sollen mithilfeder Methode der finiten Elemente Lösungen für die Feldgrößen Verschiebung und Temperaturin ausgewählten Fällen ermittelt werden.

1

2 Kontinuumstheoretische Grundlagen


Vor dem Einstieg in die Thematik soll an dieser Stelle zunächst eine kurze Einführung in ge-nutzte Schreibweisen sowie kontinuumstheoretische Grundlagen der Thermomechanik erfolgen.Tensoren beliebiger Stufe größer Null werden in koordinatenfreier Notation fettgedruckt darge-stellt, etwa σ. Die Plastizitätstheorie wurde in ihren Anfängen im frühen 20. Jahrhundert stetsin einem Orthonormalsystem, speziell in der kartesischen Standardbasis ek = e1, e2, e3 ,betrieben. Aus diesem Grund werden wir uns in dieser Arbeit ebenfalls auf die kartesische Stan-dardbasis festlegen. Ausführlich wird für einen beliebigen Tensor n-ter Stufe in dieser Basis

A = Ai1i2...inei1 ⊗ ei2 ⊗ . . .⊗ ein (2.1)

geschrieben, wobei mit a ⊗ b das dyadische Produkt zwischen zwei Vektoren a und b gemeintist und die Einsteinsche Summenkonvention verwendet wird.1 Wir wollen aus Gründen derÜbersichtlichkeit in dieser Arbeit auf die explizite Notation der Basisvektoren verzichten, dakeinerlei Mehrdeutigkeiten entstehen können. So werden wir stellenweise auch die kartesischenKomponenten Ai1i2...in eines Tensors A der Kürze halber als Tensor bezeichnen.

Unser Ziel ist die numerische Berechnung des Verschiebungsfeldes u und der Temperatur T fürsich elastisch-plastisch verformende Festkörper. Bei beiden Größen handelt es sich um Feld-größen, für welche sich Feldgleichungen aufstellen lassen. Diese Feldgleichungen werden aufGrundlage von materialunabhängigen Bilanzgleichungen konstruiert, welche mit entsprechen-den Materialgleichungen auf das interessierende Material spezifiziert werden.

Vor Einführung der Bilanzgleichungen soll jedoch zur Erklärung einiger der in ihnen vorkom-menden Größen zunächst ein kurzer Abschnitt über die Kinematik von Kontinua folgen.

2.1 Kinematik des Kontinuums

In der Kontinuumstheorie wird zwischen zwei grundlegenden Darstellungen unterschieden: derLagrangeschen Darstellung, in welcher Größen auf einen Ausgangszustand und der Euler-schen Darstellung, in der Größen auf den Momentanzustand bezogen werden. In der Festkör-permechanik ist es üblich, Felder im Lagrangeschen Sinne zu beschreiben. In der linearenTheorie, die bei der Annahme kleiner Verformungen verwendet wird, ist jedoch der Unterschiedzwischen beiden Darstellungen und infolgedessen auch zwischen den entsprechenden kinema-tischen Größen vernachlässigbar. Da sich die klassische Plastizitätstheorie stets innerhalb derTheorie kleiner Deformationen bewegte, soll an dieser Stelle eine kurze Motivation der verwen-deten Formulierungen über eine Einführung in die nichtlineare Kontinuumskinematik erfolgenund die zugehörigen Annahmen zur linearen Theorie vorgestellt werden.

1 Über wiederholte Indizes in einem Produkt- oder Differentialausdruck wird von eins bis drei summiert, etwa

aiai =3∑

i=1aiai = a1a1 + a2a2 + a3a3.

2


2.1.1 Deformations- und Verzerrungsmaße

Wir betrachten ein undeformiertes Kontinuum, dessen Teilchen jeweils eine eindeutige Aus-gangsposition Xi zu einem festen Zeitpunkt t0 in einem festen kartesischen Koordinatensystembesitzen (Abbildung 2.1 links). Dieser Zustand wird als Referenzkonfiguration bezeichnet. Wäh-rend eines beliebigen Prozesses wird das Kontinuum nun im Raum bewegt und deformiert. DieTeilchen haben ihre Positionen sowohl im Raum als auch relativ zueinander verändert. Die neuePosition des Teilchens mit ursprünglicher Lage Xi in Referenzkonfiguration wird mit xi(Xj , t)bezeichnet. Dieser Zustand zu einem beliebigen Zeitpunkt t wird als Momentankonfigurationbezeichnet. Die in diesem Kontext besonders interessierenden Deformationen eines Kontinuumslassen sich über Längenänderungen kleiner Linienelemente dxi beim Übergang von Referenz-zu Momentankonfiguration quantifizieren (Wilmański 2008). Hierzu führen wir den Deforma-

Abbildung 2.1: Zur Bedeutung des Konfigurationsbegriffs und des Deformationsgradienten

tionsgradienten F ein. Dieser schreibt sich in kartesischen Komponenten als

Fij = ∂xi∂Xj

(2.2)

und stellt eine Abbildung von Linienelementen dXi der Referenzkonfiguration auf Linienelementedxi in der Momentankonfiguration dar (siehe Abbildung 2.1), derart dass

dxi = FijdXj . (2.3)

Auch die Determinante des Deformationsgradienten spielt eine wesentliche Rolle, denn sie dientzur Transformation von Volumenelementen aus der Referenz- in die Momentankonfiguration:

dV = JdV0, J = det(F ). (2.4)

Notieren wir nun die quadratische Länge eines Linienelements in Momentankonfiguration alsFunktion der quadratischen Länge in Referenzkonfiguration mit

(dl)2 = dxidxi = FijdXjFikdX0k = dXjCjkdXk, (2.5)

3


so können wir mithilfe des Deformationsgradienten zunächst den rechten Cauchy-Green-De-formationstensor

Cjk = FijFik (2.6)

definieren. Dieser Deformationstensor ermöglicht zwar die Quantifizierung von Längenänderun-gen, ist aber als Verzerrungsmaß ungeeignet, da er bei Ausbleiben jeglicher Deformation nichtzum Nulltensor wird – ein Zusammenhang, welcher die nachfolgenden Überlegungen motiviert.Das Kontinuum gilt als unverzerrt, wenn, wie Ogden (1984) es formuliert, die Länge eines be-liebigen Linienelements in Momentankonfiguration gleich der Länge desselben Linienelements inReferenzkonfiguration ist, also in Formeln mit der Definition in Gleichung (2.6), wenn gilt:

dxjdxj − dXidXi = dXjCjkdXk − dXidXi = dXj (Cjk − δjk)dXk = 0. (2.7)

Für beliebige Linienelemente dXk 6= 0 ist also

Cjk − δjk = 0, (2.8)

zu setzen, wobei δjk das Kronecker-Delta ist, mit

δjk =

1, für j = k,

0, für j 6= k.(2.9)

Gleichung (2.7) schließt Starrkörperrotationen sowie -translationen ein und ist demnach erfüllt,solange das Kontinuum unverzerrt ist. Um bei Ausbleiben einer Deformation eine Verzerrungvon Null zu erhalten, führen wir den Green-Lagrangeschen Verzerrungstensor

Ejk = 12(Cjk − δjk) (2.10)

ein. Insbesondere quantifiziert der Green-Lagrangesche Verzerrungstensor relative Längen-änderungen bezogen auf die Referenzkonfiguration, denn es gilt mit Gleichung (2.7)

12

(dl)2 − (dl0)2

(dl0)2 = dXjEjkdXk

dXidXi. (2.11)

2.1.2 Lagrangesche und Eulersche Darstellung der Bewegung

Im Folgenden sollen ein Zusammenhang zwischen dem eingeführten Verzerrungsmaß und derVerschiebung

ui = xi −Xi (2.12)

4


gefunden und die Begriffe der Eulerschen bzw. Lagrangeschen Darstellung eingeführt wer-den. Die Verschiebung spielt eine wesentliche Rolle, da, wie noch ersichtlich werden wird, dasVerschiebungsfeld als primäre gesuchte Größe eingeführt und so der Bogen von der geometri-schen zur physikalischen Betrachtung gespannt werden kann. Wie in Abbildung 2.2 angedeutet,

Abbildung 2.2: Zur Bedeutung des Verschiebungsvektors

misst die Verschiebung die vektorielle Auslenkung eines Teilchens aus seiner Ausgangspositi-on. Mit den Definitionsgleichungen (2.2) und (2.12) kann nachgerechnet werden, dass für denGradienten der Verschiebung nach den Ortskomponenten der Referenzkonfiguration

∂ui∂Xj

= ∂xi∂Xj

− ∂Xi

∂Xj= Fij − δij (2.13)

gilt und so für den Deformationsgradienten der Zusammenhang

Fij = ∂ui∂Xj

+ δij (2.14)

aufgestellt werden kann. Mithin kann auch der rechte Cauchy-Green-Deformationstensor ausGleichung (2.6) in Abhängigkeit des Verschiebungsvektors ausgedrückt werden als

Cij =(∂ul∂Xi

+ δli

)(∂ul∂Xj

+ δlj

)= ∂ui∂Xj

+ ∂uj∂Xi

+ ∂ul∂Xi

∂ul∂Xj

+ δij . (2.15)

Setzen wir dies in die Definition des Green-Lagrangeschen Verzerrungstensors aus Glei-chung (2.10) ein, so erhalten wir diesen schließlich in Abhängigkeit der Verschiebungen als

Eij = 12

(∂ui∂Xj

+ ∂uj∂Xi

+ ∂ul∂Xi

∂ul∂Xj

). (2.16)

In dieser Darstellung wird direkt ersichtlich, weshalb dieser Verzerrungstensor den Namen La-granges trägt. Offenbar ist er in Lagrangescher Beschreibungsweise definiert. Felder, die imLagrangeschen Sinne beschrieben werden, werden auf die Referenzkonfiguration bezogen de-finiert, also beispielsweise in Abhängigkeit der ursprünglichen Teilchenpositionen Xi, wie sie indiesem Fall in den Gradienten auftreten. Insbesondere in der Festkörpermechanik wird daher

5


fast ausschließlich die Lagrangesche Darstellung verwendet, da es sich in der Anwendung alspraktikabel erweist, einen undeformierten Festkörper als Referenz zu deklarieren und gegenüberdiesem sämtliche Deformationen zu messen. Der Lagrangeschen Beschreibung steht die Eu-lersche Schreibweise gegenüber, in welcher Größen in Momentankonfiguration definiert werden.Damit kann analog zu Gleichung (2.6) über Längenänderungen der Cauchy- Deformationsten-sor cjk definiert werden:

(dL)2 = dSidSi = F−1ij dsjF−1

ik dsk, F−1ij = ∂Xi

∂xj, cjk = F−1

ij F−1ik . (2.17)

Dieser kann nach selber Argumentation wie zuvor zur Definition des Euler-Almansi-Ver-zerrungstensors

eijdef.= 1

2 (δij − cij) (2.18)

dienen und wird mit Verschiebungen ausgedrückt in die Form

eij = 12

(∂ui∂xj

+ ∂uj∂xi− ∂uk∂xi

∂uk∂xj

)(2.19)

gebracht, also zu einem Verzerrungsmaß in Bezug auf die Momentankonfiguration, sprich inEulerscher Darstellung. Es existieren zahlreiche weitere Definitionen von Verzerrungsmaßen,die hier jedoch nicht behandelt werden sollen.

2.2 Kinematik kleiner Deformationen

Die vorgestellten Verzerrungsmaße erlauben die Beschreibung beliebig großer Deformationen. Inder Theorie kleiner Verformungen werden in den verwendeten Verzerrungsmaßen nichtlineareTerme vernachlässigt, sodass der in Gleichung (2.16) vorgestellte Green-Lagrangesche Ver-zerrungstensor mit

Eij ≈12

(∂ui∂Xj

+ ∂uj∂Xi

)(2.20)

und der Euler-Almansi-Verzerrungstensor mit

eij ≈12

(∂ui∂xj

+ ∂uj∂xi

)(2.21)

abgeschätzt wird. Ferner lässt sich die Annahme kleiner Deformationen über den Deformations-gradienten aus Gleichung (2.2) als

∂ui∂Xj

δij ⇒ Fij = ∂ui∂Xj

+ δij ≈ δij (2.22)

6


ausdrücken, woraus der zusätzliche Zusammenhang

J = det(F ) ≈ det(δ) = 1 (2.23)

folgt. Mithin ist auch der Unterschied zwischen Momentan- und Referenzkonfiguration vernach-lässigbar und die beiden Verzerrungsmaße lassen sich als

Eij ≈ eij = εij = 12

(∂ui∂xj

+ ∂uj∂xi

)(2.24)

schreiben. Der in obiger Gleichung definierte sogenannte lineare Verzerrungstensor ε ist dasVerzerrungsmaß, das in der Theorie kleiner Verformungen verwendet wird und so auch in derklassischen Plastizitätstheorie verwendet wurde. Aus diesem Grund werden wir in dieser Arbeitausschließlich mit diesem Verzerrungsmaß arbeiten.

Besonderes Augenmerk ist abschließend auf den Geschwindigkeitsgradienten in Momentankon-figuration, also Eulerscher Darstellung, zu richten, der sich in kartesischen Koordinaten als

(gradv)ij = ∂vi∂xj

(2.25)

schreiben lässt. Dieser wird im anschließenden Abschnitt über Bilanzgleichungen auftreten, wes-halb wir bestrebt sind, einen Zusammenhang zu Verschiebungen oder Verzerrungen herzustel-len. Formulieren wir den linearen Verzerrungstensor nach Gleichung (2.20) in LagrangescherSchreibweise2 über die Definition der Verschiebung in Gleichung (2.12), so erhalten wir zunächst

εij (X, t) = 12

(∂xi∂Xj

+ ∂xj∂Xi

)− δij , (2.26)

sodass sich die totale Zeitableitung des linearen Verzerrungstensors zu

˙εij (X, t) = ∂ε

∂Xk

dXk

dt + ∂ε

∂t

dtdt = ∂ε

∂t= 1

2∂

∂t

(∂xi∂Xj

+ ∂xj∂Xi

)= 1

2

(∂vi∂Xj

+ ∂vj∂Xi

)(2.27)

≈ 12

(∂vi∂xj

+ ∂vj∂xi

)= ˙εij = εij (2.28)

bestimmen lässt und wir den symmetrischen Anteil des Geschwindigkeitsgradienten als zeitlicheAbleitung des linearen Verzerrungstensors verstehen können. Die Unterscheidung der Darstel-lungen spielt hier also nur in der geometrisch nichtlinearen Theorie eine Rolle.

2.3 Bilanzgleichungen

In der Literatur werden Bilanzgleichungen auf sich teils stark unterscheidende Arten eingeführt.Es gibt die strikt mathematisch formale Herangehensweise, wie sie beispielsweise bei Bertram2 An Stellen, an denen eine Unterscheidung zwischen den Darstellungen zweckmäßig ist, werden Felder in

Langrangescher Darstellung mit einer Tilde, in Eulerscher Darstellung mit einem Zirkumflex gekennzeichnet.

7


(2011) oder Wilmański (1998) praktiziert wird, aber auch die ergebnisorientierte ingenieurwis-senschaftliche Herangehensweise wie sie z. B. bei Müller (2011) oder Huang u. Khan (1995) zulesen ist. In dieser Arbeit werden wir uns an Letzterer orientieren, da wir schließlich das Ziel ver-folgen, konkrete Berechnungen für Anwendungen vorzunehmen. Folglich bauen wir nun auf derFünffeldtheorie der Kontinuumsthermomechanik auf, welche axiomatisch von fünf Feldgrößenausgeht, durch deren Bestimmung sich ein thermomechanischer Prozess vollständig beschreibenlässt. Diese sind durch die Massendichte ρ, die Geschwindkeit in drei Raumrichtungen v und dieTemperatur T gegeben. Bestimmbar werden jene Größen durch assoziierte bilanzierbare Grö-ßen: Die Geschwindigkeit v kann über die Impulsdichte ρv, die Temperatur über die spezifischeinnere Energie u bestimmt werden. Die Massendichte ρ muss nicht aus einer anderen Größe ab-geleitet, sondern kann direkt bilanziert werden. Neben diesen fünf Größen existiert eine weiterebilanzierbare, jedoch gemeinhin nicht als Basis-Feld interpretierte Zustandsgröße: die spezifischeEntropie s. Im Folgenden sollen die entsprechenden Bilanzgleichungen vorgestellt und an die indieser Arbeit verwendete Komponentendarstellung in einer kartesischen Basis angepasst werden,wobei auch die Annahme der geometrischen Linearität eingebracht wird.

2.3.1 Massenbilanz

Die GesamtmasseM in einem materiellen Volumen V (t) bleibt stets erhalten. Drücken wir dieseBedingung über das Feld der Massendichte ρ aus, dürfen wir schreiben:

dMdt = d

dt

∫V (t)

ρ dV = 0. (2.29)

Dies ist eine Bilanzgleichung für die Masse in globaler Form für ein materielles Volumen, dieaufgrund der Zeitableitung des zeitabhängigen Integrationsvolumens jedoch schwierig zu hand-haben ist. Für die Umformung in eine geeignetere Darstellung der globalen Bilanzgleichungund zum Überführen in eine lokale Form muss die Zeitableitung des zeitlich veränderlichen Vo-lumens über das Reynoldssche Transporttheorem3 ausgewertet werden. Wir erhalten damitunter Verwendung des Gaussschen Integralsatzes aus Gleichung (2.29) die Form

∫V (t)

(∂ρ

∂t+∇ · (ρv)

)dV = 0, (2.30)

wobei v die Geschwindigkeit ist und ∇ den Nabla-Operator darstellt. Da die linke Seite obigerGleichung für beliebige Volumen identisch Null werden muss, wird argumentiert, dass demnachder Integrand verschwinden muss. Wir erhalten so die lokale Form der Massenbilanz mit

∂ρ

∂t+∇ · (ρv) = 0, (2.31)

3 Für eine ausführliche Herleitung siehe zum Beispiel Müller (2011, Kapitel 3).

8


oder in kartesischen Komponenten notiert:

∂ρ

∂t+ ∂

∂xj(ρvj) = 0. (2.32)

Die Nutzung der Funktionaldeterminante J erlaubt eine wichtige Beziehung zwischen Massen-dichte in Referenzkonfiguration ρ0 und Massendichte in Momentankonfiguration ρ (nach Ogden1984):

ρ0 = Jρ, (2.33)

ein Ausdruck, der unter der Annahme kleiner Verformungen nach Gleichung (2.23) die Gleichung

ρ0 ≈ ρ (2.34)

liefert, sodass wir im Folgenden von einer konstanten Massendichte ρ0 und damit von eineridentisch erfüllten Massenbilanz ausgehen können.

2.3.2 Impulsbilanz

Nach Newtons Lex Secunda wird der Impuls P = Mv in einem System durch Krafteinwirkunggeändert:

dPdt = K = T + F . (2.35)

Der Kraftvektor K wird additiv in Oberflächenkräfte T und Volumenkräfte F zerlegt. In inte-graler Darstellung mit Impulsdichten p = ρv, Spannungsvektoren t und spezifischen Volumen-kräften f in einem materiellen Volumen schreiben wir demnach

ddt

∫V (t)

ρv =∫

∂V (t)

t dA+∫V (t)

ρf dV. (2.36)

Hierbei bezeichnet ∂V die Oberfläche des Volumens V . Nach Cauchy lässt sich der Span-nungsvektor t als lineare Abbildung des Oberflächeneinheitsnormalenvektors n deuten, wobeidie lineare Abbildung über den Cauchyschen Spannungstensor σ durch

t = n ·σ (2.37)

definiert ist. Ersetzen wir diesen Ausdruck auf der rechten Seite der Gleichung (2.36), wendenanschließend links analog zur Bilanzgleichung der Masse das Reynoldssche Transporttheoremund schließlich den Gaussschen Integralsatz an, so erhalten wir∫

V (t)

(∂

∂t(ρv) +∇ · (ρvv)− ρf

)dV =

∫∂V (t)

n ·σ dA. (2.38)

9


Eine abermalige Anwendung des Gaussschen Integralsatzes liefert nun die globale Bilanzglei-chung des Impulses der Form∫

V (t)

(∂

∂t(ρv) +∇ · (ρvv − σ)− ρf

)dV = 0. (2.39)

Nach derselben Argumentation, die bereits bei der Massenbilanz Anwendung fand, muss derIntegrand verschwinden und wir erhalten die lokale Bilanzgleichung des Impulses durch

∂

∂t(ρv) +∇ · (ρvv − σ)− ρf = 0, (2.40)

oder in Komponenten der kartesischen Basis geschrieben:

∂

∂t(ρvi) + ∂

∂xj(ρvivj − σji)− ρfi = 0. (2.41)

Da wir in dieser Arbeit an der Lösung statischer Prozesse interessiert sind und zusätzlich denEinfluss von Volumenkräften wie der Gravitationskraft vernachlässigen wollen, streichen wir dieentsprechenden Terme in obiger Gleichung und erhalten die Impulsbilanz der Statik

∂σji∂xj

= 0. (2.42)

2.3.3 Bilanz der inneren Energie

Die spezifische Gesamtenergie e in einem System, in dem keine chemische Reaktionen stattfin-den, setzt sich aus spezifischer innerer Energie u und spezifischer kinetischer Energie 1/2 (v ·v)zusammen. In Form von Energiedichten notiert schreibt sich die Gesamtenergiedichte ρe als

ρe = ρ

(u+ 1

2v ·v). (2.43)

Die Gesamtenergie in einem System bleibt erhalten und kann sich durch Wärmeaustausch mitder Umgebung über den Wärmeflussvektor q, Strahlungszufuhr über die spezifische Strahlungs-dichte r, Leistung durch Oberflächenkräfte T ·v und Leistung durch Volumenkräfte F ·v ändern.Die globale Bilanz der Gesamtenergie schreibt sich dann im Volumen:

ddt

∫V (t)

ρ

(u+ 1

2v ·v)

dV = −∫

∂V (t)

n · q dA+∫V (t)

ρr dV+

+∫

∂V (t)

n ·σ ·v dA+∫V (t)

ρf ·v dV.(2.44)

Wird nun die Bilanzgleichung des Impulses in Gleichung (2.40) skalar mit der Geschwindigkeit vmultipliziert und von der Gesamtenergiebilanz in Gleichung (2.44) abgezogen, ergibt sich daraus

10


die globale Bilanzgleichung der inneren Energie:4

ddt

∫V (t)

ρu dV = −∫

∂V (t)

q ·n dA+∫V (t)

σ · ∇v dV +∫V (t)

ρr dV. (2.45)

Unter Verwendung des Reynoldsschen Transporttheorems, des Gaussschen Integralsatzes undunter Vernachlässigung des Strahlungsterms ergibt sich die lokale Form der Bilanz der innerenEnergie zu

d(ρu)dt +∇ · q − σ · ∇v − ρf ·v = 0. (2.46)

Abermals argumentieren wir hierfür mit der Gültigkeit für beliebig kleine Volumen. Mithilfe derMassenbilanz aus Gleichung (2.32) und unter Vernachlässigung von Volumenkräften lässt sichdiese Bilanzgleichung auf die Form

ρdudt +∇ · q − σ · ∇v = 0 (2.47)

bringen, was in kartesischen Komponenten folgende Formulierung liefert:

ρdudt + ∂qi

∂xi− σji

∂vi∂xj

= 0. (2.48)

Aufgrund der Symmetrie des Cauchy-Spannungstensors geht in die doppelte Kontraktion zwi-schen Spannungstensor und Geschwindigkeitsgradient in obiger Gleichung nur der symmetrischeAnteil des Geschwindigkeitsgradienten ein. Formal können wir dies mit

σji∂vi∂xj

= σji12

(∂vi∂xj

+ ∂vj∂xi

)(2.49)

festhalten und können so den in Gleichung (2.28) aufgestellten Zusammenhang zum linearenVerzerrungstensor herstellen, sodass sich mit Gleichung (2.50) die Bilanz der inneren Energie inder in dieser Arbeit verwendeten Form zu

ρu+ ∂qi∂xi− σjiεij = 0 (2.50)

ergibt.

2.3.4 Entropiebilanz

Die Bilanzgleichung der Entropie anschaulich einzuführen stellt eine große Schwierigkeit dar, anwelcher wir uns in dieser Arbeit nicht versuchen wollen. Der rationale Weg führt über die so-genannte Clausius-Duhem-Ungleichung, die eine axiomatische Einführung durch das Postulat

4 In dieser Arbeit wird das Malzeichen · für das innere Produkt des jeweiligen Raumes verwendet. Also ist z.B.für zwei zweistufige Tensoren A und B in kartesischen Komponenten A ·B = AijBij .

11


des zweiten Hauptsatzes der Thermodynamik in der Formulierung

∫V (t)

d(ρs)dt dV ≥

∫V (t)

r

TdV −

∫∂V (t)

1Tq ·n dA. (2.51)

erlaubt (nach Truesdell u. Noll 1965). Hierbei bezeichnet s die spezifische Entropie. Nacherneuter Vernachlässigung des Strahlungsterms und unter Verwendung des Gaussschen Inte-gralsatzes können wir das Postulat auch zu

∫V (t)

(d(ρs)dt +∇ ·

(q

T

))dV = Σ,

Σ =∫V (t)

σ dV ≥ 0(2.52)

umformulieren und bezeichnen mit Σ die Entropieproduktion, bzw. mit σ die Produktionsdich-te der Entropie (nach Müller 2011). So können wir nach analoger Argumentation über dieGültigkeit für beliebige Volumina eine lokale Bilanzgleichung für die Entropie finden:

d(ρs)dt +∇ ·

(q

T

)= σ, (2.53)

oder in kartesischen Komponenten notiert:

d(ρs)dt + ∂

∂xj

(qjT

)= σ. (2.54)

Durch Verwendung der lokalen Massenbilanz in Gleichung (2.32) und der Produktregel könnenwir diesen Ausdruck weiter zu

ρTdsdt = − ∂qi

∂xi+ qiT

∂T

∂xi+ Tσ (2.55)

umformen.

12

3 Klassische Plastizität im mathematischen Modell


Wie eingangs bereits diskutiert, ist unser Ziel das Aufstellen von Feldgleichungen in den gesuch-ten Größen Verschiebung und Temperatur. Die vorgestellten Bilanzgleichungen sind unabhängigvom betrachteten Material gültig. Insbesondere der Spannungstensor muss daher über ein Ma-terialgesetz auf das interessierende Material spezifiziert, durch gesuchte Feldgrößen ausgedrücktund schließlich in die Bilanzgleichungen eingearbeitet werden, um letztlich ein Feldgleichungs-system zu erhalten, das nur die gesuchten Feldgrößen enthält. In der Festkörpermechanik wer-den Materialgesetze meist in Spannungs-Dehnungs-Zusammenhängen formuliert. Auch in derklassischen Plastizitätstheorie ist das Ziel eine mathematische Modellierung der zunächst nurexperimentell bekannten Spannungs-Dehnungs-Zusammenhänge. In diesem Abschnitt sollen dieklassischen Ideen und Modelle vorgestellt und diskutiert werden, um letztlich das Problem durcheine Kombination von Bilanzgleichungen und Materialgesetz in vollständiger Formulierung zustellen.

3.1 Die Prandtl-Reuss-Gleichungen

Die Prandtl-Reuss-Gleichungen werden auch als Grundgleichungen der zeitunabhängigenPlastizitätstheorie bezeichnet (Müller 2011) und bilden so die Basis der Modellierung plasti-scher Verformungsvorgänge. Die Überlegungen und Annahmen auf dem Weg zum vollständigenMaterialgesetz sollen im Folgenden vorgestellt werden.

3.1.1 Vorbetrachtungen

Zunächst wird die grundlegende Annahme getroffen, dass sich der Verzerrungstensor additiv ineinen elastischen und einen plastischen Anteil zerlegen lässt. Diese Annahme wurde schon durchdie Vorreiter der Plastizitätstheorie Anfang des 20. Jahrhunderts getroffen, jedoch erst durchLubliner (1972) thermodynamisch begründet. In Formeln schreiben wir

εij = εelij + εplij . (3.1)

Eine alternative Interpretation dieser Annahme liefert eine Definition des elastischen Verzer-rungstensors als

εelij = εij − εplij . (3.2)

Diese Darstellung erweist sich als nützlich, wenn man das Hookesche Gesetz wie folgt in Ra-tenform notiert:

σij = Cijklεelkl = Cijkl

(εkl − εplkl

). (3.3)

13


Aus dieser Form wird ersichtlich, dass das Ziel nachfolgender Diskussion als Beschreibung derplastischen Verzerrungsrate εplkl formuliert werden kann. Die Wahl der Ratenformulierung liegtdarin begründet, dass ein momentaner Spannungszustand bei plastischem Fließen nicht alleindurch den zugehörigen momentanen Verformungszustand bestimmt ist, wie es in der Elastizi-tätstheorie der Fall ist, sondern die gesamte Verformungsgeschichte eine Rolle spielt.

3.1.2 Assoziierte Fließregel

In der klassischen Plastizitätstheorie wird nun von der Existenz einer sog. assoziierten Fließregelausgegangen. Diese besagt, dass sich die plastische Dehnrate εplkl bis auf einen Faktor als Gradienteiner zeitunabhängigen skalaren Funktion f schreiben lässt, also formal:

εplkl = Λ ∂f

∂σkl. (3.4)

In der Literatur wird die Funktion f auch naheliegend als Fließpotential bezeichnet. Die ein-heitliche mathematische Theorie der plastischen Verformung kristalliner Werkstoffe ist auf vonMises (1928) zurückzuführen und erleichtert uns die formale Herleitung der Materialgleichun-gen. Er vereinte unter dem Konzept des Fließpotentials Theorien zur Plastizität, die schonfrüher existierten, wie die von Lévy (1871) oder Prandtl (1924), welche auf experimentellenBeobachtungen basierten.

Zur weiteren Betrachtung sind demnach das Fließpotential f und der Faktor Λ zu ermitteln.Zunächst soll das Potential f behandelt werden.

3.1.3 Fließkriterium

Aus experimentell gewonnenen Erkenntnissen ist bekannt, dass plastisches Fließen in einer bean-spruchten Probe erst ab einer gewissen Belastung einsetzt. Diese Tatsache motiviert das Ein-führen eines sogenannten Fließkriteriums. Ist dieses Kriterium erfüllt, so findet plastische Ver-formung statt. Es existieren einige Ansätze, ein solches Kriterium zu formulieren. Bei der Be-schränkung auf isotrope Medien bietet sich das von-Mises-Fließkriterium an, welches das amweitesten verbreitete Modell darstellt. Nach diesem Modell setzt Fließen ein, wenn gilt

f = 0, mit f(σ, σf) = 13σ

2vm −

13σ

2f = 1

2ΣrsΣrs −13σ

2f , (3.5)

wobei σf die momentane Fließspannung, σvm =√

32ΣrsΣrs die von-Mises-Vergleichsspannung

und Σ der Spannungsdeviator ist, der in kartesischen Komponenten durch

Σij = σij −13σrr (3.6)

definiert ist. Die im Fließkriterium verwendete Funktion f wird nun in Gleichung (3.4) alsFließpotential gesetzt. Zur Bestimmung der noch unbekannten Funktion Λ wird anschließend die

14


sogenannte Konsistenzbedingung ausgewertet. Diese stellt sicher, dass das Fließkriterium (3.5)nach dem Einsetzen plastischen Fließens zu allen Zeiten erfüllt ist und somit Änderungen in fverschwinden. Die Auswertung des totalen Differentials von f liefert dann

f(σ, σf) = ∂f

∂σijσij + ∂f

∂σfσf = 0. (3.7)

Da für f das von-Mises-Fließkriterium in Gleichung (3.5) gewählt wurde, lässt sich obige Glei-chung direkt unter Nutzung der Definition des Spannungsdeviators in Gleichung (3.6) auswerten:

0 = f(σ, σf) = ∂f

∂σijσij + ∂f

∂σfσf

= ∂f

∂Σkl

∂Σkl

∂σijσij + ∂f

∂σfσf

= ∂

∂Σkl

(12ΣrsΣrs −

13σ

2f

)∂

∂σij

(σkl −

13δklσrr

)σij −

13∂σ2

f∂σf

σf

=(1

22Σrsδkrδsl

)(δkiδlj −

13δinδjnδkl

)σij −

23σfσf

= Σkl

(δkiδlj −

13δinδjnδkl

)σij −

23σfσf

= (Σij − Σklδkl) σij −23σfσf

= (Σij − Σrr)σij −23σfσf

= Σij σij −23σfσf (3.8)

3.1.4 Verfestigungsmodell – isotrope Verfestigung

Um ausgehend vom letzten Ergebnis fortfahren zu können, muss eine Annahme über die Funk-tion getroffen werden, die die Fließspannung σf beschreibt. Modelliert wird damit die sog. Ver-festigung des Materials bei plastischer Verformung. Die Fließgrenze des Materials steigt mitplastischer Verformung, sodass eine erneut belastete Probe erst bei einer höheren Spannung, zufließen beginnt, wie in Abbildung 3.1 angedeutet. Auch hier existieren mehrere Ansätze zur Mo-dellierung, wobei die erstmals von Odqvist (1933) angenommene lineare isotrope Verfestigungdas einfachste Modell darstellt, das über die perfekte Plastizität hinausgeht. Hier wird von einemlinearen Anstieg der Fließspannung mit akkumulierter plastischer Verformung ausgegangen. InFormeln bedeutet dies:

σf(εpl)

= σf,0 + hεpl. (3.9)

Hierbei ist εpl die akkumulierte plastische Dehnung, welche der über die vergangene Zeit inte-grierten äquivalenten plastischen Dehnung nach von Mises entspricht und im eindimensionalenZugversuch gleich der totalen plastischen Dehnung εpl ist. Im allgemeinen dreidimensionalen

15


Abbildung 3.1: Zur Verfestigung – Zugprobe unter erneuter Belastung nach Entlastung im plas-tischen Bereich. Die Probe beginnt erst bei einer höheren Zugspannung σ > σf,0zu fließen.

Fall schreibt sich diese als

εpl =∫ √2

3 εplij ε

plij dt. (3.10)

Ferner ist h der sog. plastische Modul5. Dieser gibt die Steigung der Kurve σ(εpl) im eindi-mensionalen Zugversuch an, wie in Abbildung 3.2 dargestellt. Mit Gleichung (3.9) kann nun dieAbleitung σf in Gleichung (3.8) zu

σf = hεpl (3.11)

bestimmt werden. Dies ermöglicht eine weiterführende Auswertung der Konsistenzbedingung inGleichung (3.8):

Σij σij −23σfhε

pl = 0 ⇒ εpl = 3

2Σij σijhσf

. (3.12)

Abbildung 3.2: Zur Bedeutung des plastischen Moduls h

5 engl. plastic modulus z.B. in Simo u. Hughes (1998) oder Lubliner (1990)

16


Mit dem Hookeschen Gesetz aus Gleichung (3.3) ergibt sich dann hieraus

εpl = 3

2hσfΣijCijkl

(εkl − εplkl

). (3.13)

Setzt man hier die Fließregel aus Gleichung (3.4) ein und wertet anschließend die bereits zuvorexerzierte Ableitung des Fließpotentials aus, ergibt sich

εpl = 3

2hσfΣijCijkl

(εkl − Λ ∂f

∂σkl

)= 3

2hσfΣijCijkl

(εkl − ΛΣkl

). (3.14)

Um die Zeitableitung der äquivalenten plastischen Dehnrate zu eliminieren, wird die Bedingung

σklεplkl = σvmε

pl (3.15)

verwendet (vgl. etwa Huang u. Khan 1995, S. 136). Die äquivalente plastische Dissipationsleis-tung, also diejenige im eindimensionalen Zugversuch, soll gleich der tatsächlichen Dissipations-leistung sein. Da die Fließbedingung (3.5) erfüllt sein muss, wird hieraus

σklεplkl = σfε

pl. (3.16)

Mit der assoziierten Fließregel (3.4) ergibt sich so

εpl = 1

σfΛ∂f∂σ kl

σkl = 1σf

ΛΣklΣkl = 1σf

Λ23σ

2f = 2

3σfΛ. (3.17)

Dieser Zusammenhang kann nun in Gleichung (3.14) eingesetzt werden und wir erhalten

23σfΛ = 3

2hσfΣijCijkl

(εkl − ΛΣkl

)(3.18)

⇔ Λ = ΣijCijklεkl

ΣmnCmnopΣop + 49hσ

2f. (3.19)

Mit Gleichung (3.4) folgt dann für die plastische Dehnrate

εplkl = ΣrsCrstuεtu


2f

Σkl. (3.20)

Schließlich können wir mit Gleichung (3.3) für das Spannungsinkrement

σij =(Cijkl −

CijtuΣtuΣrsCrskl


2f

)εkl (3.21)

schreiben und erhalten so einen Spannungs-Dehnungszusammenhang in Ratenformulierung.

Gleichung (3.21) können wir nun in eine kompaktere Form überführen, indem wir den Stei-figkeitstensor C unter Zuhilfenahme der Lamé-Konstanten λ und µ auf isotrope Materialien

17


spezifizieren, sodass

Cijkl = λδijδkl + µδikδjl + µδilδjk. (3.22)

Setzen wir diese Vereinfachung in Gleichung (3.21) ein, ergibt sich

σij = (λδijδkl + µδikδjl + µδilδjk) εkl−

− (λδijδtu + µδitδju + µδiuδjt) ΣtuΣrs (λδrsδkl + µδrkδsl + µδrlδsk)Σmn (λδijδkl + µδikδjl + µδilδjk) Σop + 4

9hσ2f

.(3.23)

Da der Spannungsdeviator spurfrei und symmetrisch ist, also Σkk = 0 und Σtu = Σut, ergibtsich hieraus

σij = λεrrδij + 2µεij −2µΣijΣklεkl

ΣmnΣmn + 29hµσ

2f. (3.24)

Im nächsten Schritt wird die Fließbedingung aus Gleichung (3.5) in der Form

ΣijΣij = 23σ

2f (3.25)

eingesetzt. Dies liefert eine kompakte Formulierung mithilfe der Lamé-Konstanten der Form

σij = λεrrδij + 2µεij −3µ

1 + h3µ

ΣijΣkl

σ2f

εkl. (3.26)

Auch kann mit Hilfe des Elastizitätsmoduls E und der Poisson-Zahl ν unter Nutzung derZusammenhänge

λ = Eν

(1 + ν)(1− 2ν) ,

µ = E

2(1 + ν) ,(3.27)

folgende Formulierung aus Gleichung (3.26) gewonnen werden:

σij = E

1 + ν

(εij + ν

1− 2ν εrrδij −32

3E3E + 2h(1 + ν)

ΣijΣkl

σ2f

εkl

). (3.28)

Die drei äquivalenten Formulierungen des Prandtl-Reuss-Materialgesetzes in den Gleichun-gen (3.21), (3.26) und (3.28) erlauben es, Spannungsinkremente eines plastisch fließenden Mate-rials zu ermitteln. Wir fassen zusammen und ergänzen die eben genannten Materialgleichungenfür die plastische Verformung um das Hookesche Gesetz zur Modellierung des elastischen Be-reichs und notieren zunächst in Schreibweise mit Lamé-Parametern aus Gleichung (3.26) unter

18


Nutzung des von-Mises-Fließkriteriums in Gleichung (3.5)

σij =

λεrrδij + 2µεij für f < 0,

λεrrδij + 2µεij −3µ

1 + h3µ

ΣijΣkl

σ2f

εkl für f = 0.(3.29)

Dies sind die Prandtl-Reuss-Gleichungen zur Modellierung elastisch-plastischer Verformungs-vorgänge. Enthalten sind das von-Mises-Fließkriterium und ein lineares isotropes Verfestigungs-modell.

3.1.5 Grenzen des Modells – der Bauschinger-Effekt

Die Prandtl-Reuss-Gleichungen ermöglichen zwar die korrekte Vorhersage der wachsendenFließgrenze bei monotoner Belastung. Wird jedoch eine beispielsweise bereits durch Zug plas-tisch verformte Probe unter Druckbeanspruchung erneut plastisch verformt, so ist die Initial-fließspannung entgegen der Voraussage des linear-isotropen Verfestigungsmodells betragsmäßigkleiner als im ersten Verformungsvorgang. Dieses Phänomen wurde zuerst von Bauschinger(1886) experimentell untersucht und wird daher auch als Bauschinger-Effekt bezeichnet. InAbbildung 3.3(a) sind Versuchsdaten von Žerovnik u. a. (2010) dargestellt, die den Bauschin-ger-Effekt gut widerspiegeln. Diese sind dem schematisch dargestellten isotropen Verfestigungs-modell unter zyklischer Belastung in Abbildung 3.3(b) gegenübergestellt. Trotz der fehlerhaften

(a) Bauschinger-Effekt im Experiment (Quelle:Žerovnik u. a. 2010)

(b) Isotropes Verfestigungsmodell – idealisierter1D-Zugversuch

Abbildung 3.3: Zum Bauschinger-Effekt bei zyklischer Belastung

Modellierung bei wechselseitiger Belasting wird Letzteres häufig insbesondere wegen seiner Ein-fachheit genutzt, dann jedoch nur zur Modellierung monoton belasteter Systeme (Doghri 2000).Um den Bauschinger-Effekt korrekt zu modellieren bedarf es eines erweiterten Verfestigungs-modells.

19


3.2 Der Bauschinger-Effekt im Modell

Betrachten wir die Fließregel aus Gleichung (3.5) und setzen das lineare isotrope Verfestigungs-modell aus Gleichung (3.9) ein, erhalten wir das in den Prandtl-Reuss-Gleichungen verwen-dete Fließpotential f in der Form

f(σij , ε

pl)

= 32ΣijΣij −

(σf,0 + hεpl

)(3.30)

und können die allgemeine Struktur dieser Funktion als

f(σij , ε

plij

)= F (σij)− k

(εpl)

(3.31)

notieren. Gleichung (3.31) stellt die allgemeine Struktur eines isotropen Verfestigungsmodellsdar (Lubliner 1990). Deutlich ist, dass die Funktion k unabhängig von der Belastung desMaterials während plastischen Fließens monoton anwachsen muss um das Fließkriterium f = 0zu erfüllen. So ergibt sich das zuvor genannte und in Abbildung 3.3b dargestellte Phänomender stets wachsenden Fließgrenze. Um dieser Problematik zu begegnen, hat Melan (1938) einFließpotential in der Form

f(σij , ε

plij

)= 1

2

(Σij − αij

)(Σij − αij

)− 1

3σ2f,0 (3.32)

vorgeschlagen. Insbesondere geht das einfachste kinematische Verfestigungsmodell hierauf auf-bauend auf Prager (1955) zurück, der für den spurfreien sogenannten Rückspannungstensor6

den Zusammenhang

αij = cεplij (3.33)

vorschlug. Der Begriff der kinematischen Verfestigung wurde ebenfalls erst durch Prager (1955)eingeführt, der fiktive technische Apparaturen und ihre kinematischen Zusammenhänge zur Ver-anschaulichung des mathematischen Formalismus nutzte.7 Auf Pragers Modell aufbauend, hatZiegler (1959) nach eingehender Analyse (Shield u. Ziegler 1958) eine Modifikation derForm

αij = (σij − αij) Γ, Γ > 0, (3.34)

vorgeschlagen. Aufbauend hierauf soll im Folgenden das entsprechende Materialgesetz für daskinematische Verfestigungsmodell entwickelt werden.

Die Vorgehensweise zum Auffinden des Materialgesetzes mit implementierter kinematischerVerfestigung ähnelt der Vorgehensweise, die im vorigen Abschnitt zu den Prandtl-Reuss-Gleichungen vorgestellt wurde. Ein wichtiger Unterschied besteht in der Wahl des Fließpotenti-

6 engl.: back stress7 „To avoid lengthy though elementary algebra [...]“, in Prager (1955)

20


als, das sich aufgrund des veränderten Verfestigungsmodells ebenfalls ändert. Wie zuvor startenwir bei der assoziierten Fließregel

εplkl = Λ ∂f

∂σkl, (3.35)

welche in Gleichung (3.4) vorgestellt wurde, verwenden hier jedoch den in der Fließregel ausGleichung (3.32) verwendeten Ausdruck als Fließpotential:

f = 12

(Σij − αij

)(Σij − αij

)− 1

3σ2f,0, (3.36)

mit der zusätzlichen Entwicklungsgleichung nach Ziegler aus Gleichung (3.34). Weiterhin stel-len wir wie im vorigen Abschnitt die Konsistenzbedingung auf, dieses Mal jedoch, um den ska-laren Faktor Γ zu bestimmen:

f = ∂f

∂σijσij + ∂f

∂αklαkl = 0. (3.37)

Die sowohl in der Fließregel in Gleichung (3.35) als auch in der Konsistenzbedingung in Glei-chung (3.37) auftretende Ableitung nach dem Spannungstensor berechnet sich zu

∂f

∂σij= ∂f

∂Σkl

∂Σkl

∂σij

= (Σkl − αkl)(δkiδlj −

13δklδij

)= Σij − αij .

(3.38)

Hierbei wurde beachtet, dass sowohl Σij als auch αij spurfrei ist, also Σii = αkk = 0. Setzen wirdieses Ergebnis in Gleichung (3.35) ein, so erhalten wir

εplij = Λ (Σij − αij) . (3.39)

Weiter ergibt sich mit

∂f

∂αkl= (Σij − αij) (−δikδjl) = − (Σkl − αkl) (3.40)

aus der Konsistenzbedingung in Gleichung (3.37) der Ausdruck

f = (Σij − αij) σij − (Σkl − αkl) αkl = 0. (3.41)

Setzen wir hier die Entwicklungsgleichung für αkl aus Gleichung (3.34) ein, finden wir

(Σij − αij) σij = (Σkl − αkl) (σkl − αkl) Γ (3.42)

21


und bestimmen hieraus unmittelbar die Ratengleichung

Γ = (Σij − αij) σij(Σkl − αkl) (σkl − αkl)

. (3.43)

Beachten wir nun, dass in die Kontraktion eines Tensors mit einem spurfreien Tensor nur derdeviatorische Anteil Ersteren eingeht, dürfen wir den Nenner in obiger Gleichung umschreibenund erhalten

Γ = (Σij − αij) σij(Σkl − αkl) (Σkl − αkl)

. (3.44)

Hier identifizieren wir im Nenner die modifizierte Vergleichsspannung aus Gleichung (3.32),sodass wir obigen Ausdruck für den Fall f = 0, also für plastisch fließendes Material weiterumformen können zu

Γ = (Σij − αij) σij23σ

2f,0

. (3.45)

Da diese Ratengleichung die Erfüllung der Konsistenzbedingung sicherstellt, kann die skalareFunktion Λ frei gewählt werden, sodass hier das Resultat für Pragers ursprüngliches Rück-spannungsmodell in Gleichung (3.33) verwendet werden darf, welches sich zu

Λ = (Σij − αij)49hkσ

2f,0

σij (3.46)

berechnen lässt (Ziegler 1959). Die Konstante c haben wir durch eine andere Konstante 23hk

ersetzt, um eine spätere Analogie zum isotropen Verfestigungsmodell herzustellen.

Um nun an Stelle der Abhängigkeit von der Spannungsrate eine Abhängigkeit von der Dehnrateim obigen Ausdruck zu erzeugen, verwenden wir das Hookesche Gesetz aus Gleichung (3.3) undermitteln so zunächst

Λ =

(Σij − αij

)Cijkl

(εkl − εplkl

)49hkσ

2f,0

. (3.47)

Setzen wir für die plastische Dehnrate wieder die assoziierte Fließregel aus Gleichung (3.39) ein,erhalten wir

Λ =

(Σij − αij

)Cijkl

(εkl − Λ (Σkl − αkl)

)49hkσ

2f,0

. (3.48)

Nach kurzer Rechnung lässt sich die Funktion Λ isolieren und wir erhalten

Λ = (Σij − αij)Cijklεkl49hkσ

2f,0 + (Σmn − αmn)Cmnop (Σop − αop)

, (3.49)

was in die assoziierte Fließregel aus Gleichung (3.39) eingesetzt die plastische Dehnrate in Ab-

22


hängigkeit von der totalen Dehnrate liefert mit

εplij = (Σrs − αrs)Crstu (Σij − αij)49hkσ


εtu. (3.50)

Mit dem Hookeschen Gesetz aus Gleichung (3.3) können wir so letztlich den inkrementel-len Spannungs-Dehnungszusammenhang inklusive eines linearen kinematischen Verfestigungs-modells wie folgt formulieren:

σij =(Cijkl −

Cijtu (Σtu − αtu) (Σrs − αrs)Crskl49hkσ


)εkl. (3.51)

Ergänzen wir nun das in obiger Gleichung beschriebenene Verhalten eines plastisch fließendenMaterials analog zu Gleichung (3.29) um den linear elastischen Hookeschen Bereich, so könnenwir schreiben:

σij =

Cijklεkl für f < 0,(Cijkl −



)εkl für f = 0,

αij = (σij − αij) Γ, Γ = (Σij − αij) σij23σ

2f,0

.

(3.52)

23

4 Analytische Berechnung von Spannungs-Dehnungs-Beziehungen

4 Analytische Berechnung von Spannungs-Dehnungs-Beziehun-gen

Im Folgenden werden die Prandtl-Reuss-Gleichungen unter Annahme linearer isotroper Ver-festigung, sprich in der Form der Gleichungen (3.21), (3.26) und (3.28), auf die zwei Mo-dellprobleme Scherklotz und 1D-Zugstab angewendet, um auf analytischem Wege Spannungs-Dehnungs-Zusammenhänge zu ermitteln und diese mit späteren numerischen Lösungen zu ver-gleichen. Wenngleich für die Berechnung eines Spannungs-Dehnungs-Zusammenhangs keine Lö-sung der Bilanzgleichung des Impulses gefunden werden muss, wird hier für den Scherklotzzunächst über einen semi-inversen Ansatz ein Verschiebungsfeld als Lösung der Impulsbilanzermittelt und daraus schließlich der Spannungs-Dehnungszusammenhang abgeleitet. Anschlie-ßend wird sowohl für den Scherklotz als auch für den 1D-Zugstab die direkte Berechnung unterNutzung des Prandtl-Reuss-Materialgesetzes vorgestellt.

4.1 Modellproblem Scherklotz – Lösung der Impulsbilanz unter semi-inver-sem Ansatz

Mit dem Ziel einen Spannungs-Dehnungszusammenhang der Form σ12(ε12) zu finden, unter-suchen wir den Scherklotz in Abbildung 4.1. Dieser ist einer gegebenen zeitlich konstantenhorizontalen Verschiebungsrate uoben1 auf der oberen Seite ausgesetzt. Gestartet wird bei der

Abbildung 4.1: Modellproblem Scherklotz – Verschiebungsrate als Randbedingung

plastischen Dehnrate nach Prandtl und Reuss aus Gleichung (3.20). Mit der Einschränkungauf ein isotropes Material ergibt sich dann mit Gleichung (3.22) zunächst

εplkl = Σrs (λδrsδtu + µδrtδsu + µδruδst) εtuΣkl

Σop (λδopδmn + µδomδpn + µδonδpm) Σmn + 49hσ

2f. (4.1)

Eine weitere Auswertung unter Nutzung der Definition des Kronecker-Delta liefert

εplkl = (λΣkkεmm + µΣtuεtu + µΣutεtu) Σkl

λΣkkΣmm + µΣmnΣmn + µΣnmΣmn + 49hσ

2f. (4.2)

24


Da der Spannungsdeviator spurfrei und symmetrisch ist, also Σkk = 0 und Σtu = Σut, ergibtsich hieraus

εplkl = 2µΣtuεtuΣkl

2µΣmnΣmn + 49hσ

2f. (4.3)

Im nächsten Schritt wird die Fließbedingung aus Gleichung (3.5) in der Form

f = 0 ⇔ σ2f = 3

2ΣijΣij (4.4)

eingesetzt. Dies liefert

εplkl = 2µΣtuεtuΣkl

2(µ+ 1

3h)

ΣmnΣmn

. (4.5)

Wird nun die Definition des Spannungsdeviators eingesetzt, erhalten wir

εplkl =2µ(σtu − 1

3σrrδtu)εtu(σkl − 1

3σrrδkl)

2(µ+ 1

3h) (σmnσmn − 1

3 (σrr)2)

=2µ(σtuεtu − 1

3σrrεii) (σkl − 1

3σrrδkl)

2(µ+ 1

3h) (σmnσmn − 1

3 (σrr)2) . (4.6)

Motiviert durch den skizzierten Verformungszustand (rot) wird nun für das Verschiebungsfeldui der semi-inverse Ansatz

ui ≡

u1(x2, t)

00

(4.7)

gewählt. Für den Verzerrungstensor ergibt sich hieraus

εij ≡

0 1

2∂u1∂x2

012∂u1∂x2

0 00 0 0

(4.8)

und für die Dehnrate

εij ≡

0 1

2∂u1∂x2

012∂u1∂x2

0 00 0 0

. (4.9)

Die Spur der Dehnrate wird somit εkk = 0. Auch ist die Spur des Spannungstensors σrr = 0,

25


denn der vorliegende Spannungszustand wird als kollinear angenommen:

σkl ≡

0 σ12 0σ12 0 00 0 0

. (4.10)

Diese Zusammenhänge bedeuten, dass sich Gleichung (4.6) zu

εplkl = µσ12ε12σkl(µ+ 1

3h)

(σ12)2(4.11)

vereinfachen lässt. Schreibt man nun das Hookesche Gesetz (3.3) ebenfalls unter Nutzung vonGleichung (3.22) um, so erhalten wir zunächst

σkl = λεjjδkl − λεpljjδkl + 2µεkl − 2µεplkl. (4.12)

Aus den Gleichungen (4.10) und (4.11) wird ersichtlich, dass neben der Spur der Dehnrate auchdie Spur der plastischen Dehnrate Null ist, sodass sich der Ausdruck

σkl = 2µεkl − 2µεplkl. (4.13)

ergibt. Hier setzen wir nun die plastische Dehnrate aus Gleichung (4.11) ein. Es folgt also, dass

σkl = 2µεkl −2µ2σ12ε12σkl(µ+ 1

3h)

(σ12)2. (4.14)

Von dieser Form ausgehend spezifizieren wir weiter und werten für k = 1 und l = 2 aus:

σ12 = 2 µh

3µ+ hε12. (4.15)

Nun wird hier die Dehnrate aus Gleichung (4.9) eingesetzt, sodass wir

σ12 = µh

3µ+ h

∂u1∂x2

(4.16)

erhalten. Diese Form des Spannungstensors soll die Impulsbilanz aus Gleichung (2.42) erfüllen.Demnach folgt, dass

∂σji∂xj

= ∂σ1i∂x1

+ ∂σ2i∂x2

= 0. (4.17)

Für i = 1, 2 ausgewertet ergeben sich die beiden Gleichungen

∂σ12∂x2

= 0 und ∂σ12∂x1

= 0. (4.18)

26


Wird hier nun die Scherrate aus Gleichung (4.16) eingesetzt, finden wir

µh

3µ+ h

∂2u1∂x2∂x2

= 0, (4.19)

µh

3µ+ h

∂2u1∂x1∂x2

= 0. (4.20)

Gleichung (4.20) ist wegen des Ansatzes in Gleichung (4.7) identisch erfüllt. Es verbleibt die zulösende Gleichung (4.19). Diese impliziert, dass die Zeitableitung der Verschiebung u1 linear inx2 sein muss:

u1 = C1x2 + C2. (4.21)

Eine der beiden nun benötigten Randbedingungen ist offensichtlich

u1(x2 = 0, t

)= 0. (4.22)

Geben wir nun zusätzlich eine konstante Geschwindigkeit vor, mit welcher am oberen Rand desKlotzes geschert wird, so können wir daraus die zweite Randbedingung ableiten:

u1(x2 = b, t

)= uoben1 . (4.23)

So können wir die Konstanten in Gleichung (4.21) bestimmen und finden

u1 = uoben1b

x2. (4.24)

Setzen wir diese Verschiebung nun zurück in Gleichung (4.16) ein, erhalten wir

σ12 = µh

3µ+ h

uoben1b

. (4.25)

Die Integration in der Zeit liefert dann

σ12 = µh

3µ+ h

uoben1b

t+ C3. (4.26)

Geben wir nun vor, dass plastisches Fließen zum ersten Mal zum Zeitpunkt t = tfließ einsetzt,so erhalten wir aus der Bedingung

σ12(t = tfließ

)= σfl,012 (4.27)

die Spannungsfunktion in Abhängigkeit der Zeit:

σ12(t) = µh

3µ+ h

uoben1b

(t− tfließ) + σfl,012 , (4.28)

wobei σfl,012 diejenige Scherspannung ist, bei der Fließen zum ersten Mal einsetzt. Diese kann mit

27


Gleichung (4.10) und der Fließbedingung (4.4) zu

σfl,012 = 1√3σ0f (4.29)

umformuliert werden. Gleichermaßen können wir Gleichung (4.24) in Gleichung (4.9) einsetzen,um die Verzerrungsfunktion zu bestimmen. Zunächst erhalten wir

ε12 = uoben12b . (4.30)

Die Integration liefert schließlich

ε12 = uoben12b t+ C4, (4.31)

sodass wir mit der Bedingung

ε12(t = tfließ

)= εfl,012 (4.32)

die Verzerrungsfunktion zu

ε12(t) = uoben12b (t− tfließ) + εfl,012 (4.33)

bestimmen können. Hierbei ist εfl,012 diejenige Verzerrung, bei welcher zum ersten Mal Fließeneinsetzt. Obige Gleichung können wir nun in Gleichung (4.28) einsetzen, sodass wir mit Glei-chung (4.29) die Funktion

σ12 (ε12) = 2µh3µ+ h

(ε12 − εfl,012

)+ 1√

3σ0f (4.34)

erhalten. Der plastische Bereich des Spannungs-Dehnungs-Diagramms ist nun mit der letztenGleichung beschrieben. Für ein vollständiges Spannungs-Dehnungs-Diagramm verwenden wirfür den Bereich linearer Elastizität den bekannten Zusammenhang

σ12 (ε12) = 2µε12. (4.35)

Der volle Spannungs-Dehnungszusammenhang kann nun als

σ12 (ε12) =

2µε12, für σ12 ≤

1√3σ0f ,

2µh3µ+ h

(ε12 − εfl,012

)+ 1√

3σ0f , für σ12 >

1√3σ0f

(4.36)

geschrieben werden.

28


4.2 Modellproblem Scherklotz – Direkte Berechnung

Nach der Lösung der Impulsbilanz im vorigen Abschnitt werden wir nun denselben Scherklotzerneut untersuchen, jedoch ohne die Impulsbilanz zu berücksichtigen. Der Einfachheit halberwurden die Randbedingungen angepasst, sodass wie in Abbildung 4.2 angedeutet eine Scher-spannungsrate τ am oberen Rand vorgegeben ist. Für die tensoriellen Spannungs- und Dehnratennotieren wir

σij ≡

0 τ 0τ 0 00 0 0

, εij ≡

ε11 ε12 ε13

ε21 ε22 ε23

ε31 ε32 ε33

. (4.37)

Der Spannungsdeviator ist wegen der Spurfreiheit des Spannungstensors identisch mit selbigem

Abbildung 4.2: Modellproblem Scherklotz – Scherspannungsrate als Randbedingung

und lässt sich nach Zeitintegration bestimmen:

Σij = σij ≡

0 τ 0τ 0 00 0 0

, τ =t∫

t=0

τ(t)

dt. (4.38)

Unter Nutzung der Symmetrieeigenschaft des Dehnratentensors, also ε12 = ε21, finden wir

Σklεkl = 2τ ε12 (4.39)

und erhalten so mit dem Prandtl-Reuss-Materialgesetz aus Gleichung (3.26) nach Auswertungfür i = 1 und j = 2 die Gleichung

τ =(

2µ− 3µ1 + h

3µ

2τ2

σ2f

)ε12. (4.40)

Da das Material nach Voraussetzung bereits plastisch deformiert ist, können wir mithilfe derFließbedingung aus Gleichung (3.5) finden, dass die Identität

σ2f = 3

2ΣklΣkl = 3τ2 (4.41)

29


gelten muss und sich Gleichung (4.40) mithin zu

τ = 2µh3µ+ h

ε12 (4.42)

vereinfachen lässt. Wird diese Gleichung in der Zeit integriert und um den linear-elastischenBereich unter Nutzung des Hookeschen Gesetzes ergänzt, ergibt sich mit dem Fließkriteriumin Gleichung (4.41) der vollständige Spannungs-Dehnungszusammenhang zu

τ (ε12) =

2µε12, für τ < 1√

3σ0f ,

2µh3µ+ h

(ε12 − εfl,012

)+ 1√

3σ0f , für τ = 1√

3σf.

(4.43)

Hierbei ist εfl,012 der zur Initialfließscherspannung τ = 1√3σ

0f gehörige Dehnungswert. In Ab-

schnitt 7 wird eine Darstellung der ermittelten Lösungen geliefert und einer numerischen Lösunggegenübergestellt.

4.3 Modellproblem Zugstab

Betrachtet wird nun der Zugstab in Abbildung 4.3. Dieser unterliege einem eindimensiona-len Zugspannungszustand durch die vorgegebene zeitabhängige Spannungsrate σ(t) und ha-be bereits die Initialfließspannung überschritten. Gesucht ist wieder ein Spannungs-Dehnungs-

Abbildung 4.3: Modellproblem Zugstab

Zusammenhang, diesmal der Form σ(ε11). Für die tensoriellen Spannungs- und Dehnraten schrei-ben wir

σij ≡

σ 0 00 0 00 0 0

, εij ≡

ε11 ε12 ε13

ε21 ε22 ε23

ε31 ε32 ε33

. (4.44)

30


Mit dem Prandtl-Reuss-Materialgesetz aus Gleichung (3.26) lässt sich so das folgende lineareGleichungssystem durch Auswertung für ij = 11, 22, 33 aufstellen:

σ = λ (ε11 + ε22 + ε33) + 2µε11 −23

3µ1 + h

3µ

(σ

σf

)2 (23 ε11 −

13 ε22 −

13 ε33

), (4.45a)

0 = λ (ε11 + ε22 + ε33) + 2µε22 + 13

3µ1 + h

3µ

(σ

σf

)2 (23 ε11 −

13 ε22 −

13 ε33

), (4.45b)

0 = λ (ε11 + ε22 + ε33) + 2µε33 + 13

3µ1 + h

3µ

(σ

σf

)2 (23 ε11 −

13 ε22 −

13 ε33

). (4.45c)

Hierfür wurde der Spannungsdeviator zu

Σij = σij −13σkkδij ≡

23σ 0 00 −1

3σ 00 0 −1

3σ

, σ =t∫

t=0

σ(t)

dt (4.46)

berechnet und eingesetzt. Beachtet man nun, dass nach dem Fließkriterium in Gleichung (3.5)während des Fließens

σ = σf (4.47)

gelten muss und nutzen die Abkürzung

ξ = 3µ1 + h

3µ, (4.48)

vereinfacht sich das lineare Gleichungssystem (4.45) zu

σ = λ (ε11 + ε22 + ε33) + 2µε11 −23ξ(2

3 ε11 −13 ε22 −

13 ε33

), (4.49a)

0 = λ (ε11 + ε22 + ε33) + 2µε22 + 13ξ(2

3 ε11 −13 ε22 −

13 ε33

), (4.49b)

0 = λ (ε11 + ε22 + ε33) + 2µε33 + 13ξ(2

3 ε11 −13 ε22 −

13 ε33

). (4.49c)

Aus Gleichungen (4.49b) und (4.49c) folgt unmittelbar der Zusammenhang

ε22 = ε33. (4.50)

So lässt sich durch Einsetzen in Gleichungen (4.49a) und (4.49b) finden, dass

σ =(λ+ 2µ− 4

9ξ)ε11 +

(2λ+ 4

9ξ)ε22, (4.51a)

0 =(λ+ 2

9ξ)ε11 +

(2λ+ 2µ− 2

9ξ)ε22. (4.51b)

31


Durch Kombination beider Gleichungen erhalten wir

ε22 = −λ+ 2

9ξ

2λ+ 2µ− 29ξε11. (4.52)

Dies in Gleichung (4.51a) eingesetzt liefert letztlich

σ = µ(3λ+ 2µ)1− 1

1 + h/3µ

λ+ µ

(1− 1

31

1 + h/3µ

) ε11. (4.53)

Wie schon am Beispiel des Scherklotzes im vorigen Abschnitt exerziert kann nun mithilfe desHookeschen Gesetzes der volle Spannungs-Dehnungs-Zusammenhang geschrieben werden als

σ =

Eε11, für σ < σ0f ,

µ(3λ+ 2µ)1−

11 + h/3µ

λ+µ(

1−13

11 + h/3µ

) (ε11 − εfl,011

)+ σ0

f . für σ = σf.(4.54)

Auch für den Zugstab wird in Abschnitt 7 eine Darstellung der ermittelten Lösung geliefert undeiner numerischen Lösung gegenübergestellt.

32

5 Thermische Erweiterung


Bislang haben wir in Anlehnung an die klassische Plastizitätstheorie nur die rein mechanischeEbene der Plastizität diskutiert. Praktisch gesprochen haben wir also bislang ohne die Zu-standsgröße Temperatur gearbeitet. Die Erweiterung der bestehenden Gleichungen um einenTemperatureinfluss soll im folgenden Kapitel geschehen, mit dem Ziel, auch die Bilanz der inne-ren Energie auf ein plastisch fließendes Material zu spezifizieren und die bisher einzige gesuchteGröße Verschiebung um eine weitere gesuchte Feldgröße – die Temperatur – zu ergänzen.

5.1 Kopplung von innerer Energie und Entropie

Unser Ziel ist das Aufstellen einer partiellen Differentialgleichung zur Bestimmung der skala-ren Feldgröße Temperatur. Hierfür müssen sowohl die Bilanz der inneren Energie, als auch dieEntropiebilanz verwendet werden. Im Folgenden sollen beide Gleichungen gekoppelt und alleauftretenden Größen so ausgedrückt werden, dass letztlich nur noch Verschiebungen und dieTemperatur in einer partiellen Differentialgleichung auftreten.

5.1.1 Grundgleichungen und Annahmen

Zunächst wird eine Annahme bezüglich des Spannungsleistungterms σij εij in Energie- undEntropiebilanz getroffen. Nach Collins u. Houlsby (1997) existiert eine additive Zerlegungdes Spannungsleistungterms in einen rein reversiblen und einen rein irreversiblen Teil der Form

σij εij = σij εrevij︸︷︷︸

reversibel

+ σij εirrij .︸︷︷︸

irreversibel

(5.1)

Ferner nehmen wir an, dass der irreversible Dehnratenanteil εirr allein durch plastische Verfor-mung zustande kommt:

εirrkl = εplkl. (5.2)

Setzt man mit diesen Annahmen die lokale Entropiebilanz aus Gleichung (2.55) in die lokaleBilanz der inneren Energie in Gleichung (2.50) ein und teilt durch die Temperatur T , ergibt sich

1Tρu− ρs = − qi

T 2∂T

∂xi+ 1Tσij ε

revij + 1

Tσij ε

plij − σ. (5.3)

Für die Entropieproduktion kann nach Ziegler u. Wehrli (1987) die folgende Annahme ge-troffen werden:

σ = 1T

(−qiT

∂T

∂xi+ σij ε

plij

). (5.4)

33


Da nach Gleichung (2.52) die Entropieproduktionsdichte σ positiv sein muss, also auch im Fallenicht vorhandener plastischer Verformung εplij , können wir anhand obiger Gleichung das Fouri-ersche Wärmeleitungsgesetz in Gleichung (5.5) motivieren. So ist für

qi = −κ ∂T∂xi

, κ > 0 (5.5)

die Entropieproduktion nach Gleichung (5.4) für fernbleibende plastische Verformungen stetspositiv.

Zudem erhalten wir mit Gleichung (5.4) aus Gleichung (5.3) den Ausdruck

ρu = ρT s+ σij εrevij (5.6)

und hieraus in differentieller Notation die Gibbssche Gleichung

ρdu = ρTds+ σijdεrevij . (5.7)

Diese Gleichung beschreibt die Änderung der spezifischen inneren Energie u als Funktion derÄnderungen von spezifischer Entropie s und reversibler Dehnung εrev .

5.1.2 Zustandsraumwechsel und Materialgleichungen

Der Zustandsraum ist nach Gleichung (5.7) durch Z = s, εrev gegeben und beinhaltet damitdie schwer messbare Größe Entropie. Zudem ist noch keine Materialgleichung für den Span-nungstensor σ gefunden oder festgelegt worden. Um den Spannungsterm zu ersetzen, also eineMaterialgleichung zu finden, führen wir nun die Gibbssche freie Energiedichte

φ = ρu− ρTs− σijεrevij (5.8)

ein, deren Änderung

dφ = ρdu− ρdTs− ρTds− dσijεrevij − σijdεrevij (5.9)

mit Gleichung (5.7) durch Änderungen der Spannung σij und der Temperatur T ausgedrücktwerden kann:

dφ = −SdT − εrevij dσij . (5.10)

34


Durch diese Transformation haben wir den Zustandsraum zu Z = T,σ gewechselt und könnenso mit S = ρs die folgenden totalen Differentiale notieren:

dS = ∂S

∂T︸︷︷︸def.= Cp/T

dT + ∂S

∂σij︸︷︷︸def.= pij

dσij , (5.11a)

dεrevij =∂εrevij∂T︸︷︷︸

def.= αij

dT +∂εrevij∂σkl︸︷︷︸

def.= Sijkl

dσkl. (5.11b)

Die Wahl der in Gleichung (5.11a) getroffenen Definition

∂S

∂T= Cp

T(5.12)

findet ihre Begründung in der Thermostatik von Gasen, in der die spezifische Wärmekapazitätbei konstantem Druck als

Cp = ∂Q

∂T(5.13)

definiert wird. Über den ersten Hauptsatz der Thermodynamik für ein spannungsfreies Systemin der Form

1TdQ = dS (5.14)

und die totalen Differentiale

dS = ∂S

∂TdT (5.15)

dQ = ∂Q

∂TdT (5.16)

kann so ein Zusammenhang der Form

1TCpdT = ∂S

∂TdT (5.17)

gefunden werden, sodass

∂S

∂T= Cp

T. (5.18)

gilt. Um einen Ausdruck für den Spannungstensor zu finden und die Definitionen der partiellenAbleitungen in Gleichung (5.11b) zu motivieren, wird diese nun aus einem spannungs- unddehnungsfreien Referenzzustand mit der Temperatur Tref integriert und es entsteht der Ausdruck

εrevij = αij (T − Tref) + Sijklσkl. (5.19)

35


Wird nun Gleichung (5.19) nach dem Spannungstensor umgestellt, ergibt sich

σkl = Cijkl(εrevij − αij (T − Tref)

), (5.20)

wobei sich der Zusammenhang zwischen dem SteifigkeitstensorC und dem NachgiebigkeitstensorS durch die Umformung von Gleichung (5.19) nach Gleichung (5.20) ergibt. Hier lässt sich unterder geläufigen Annahme der Subsymmetrien des Steifigkeitstensors C, also

Cijkl = Cklij , (5.21)

das um thermische Spannungen erweiterte Hookesche Gesetz nach Duhamel-Neumann iden-tifizieren:

σij = Cijkl (εrevkl − αTδkl (T − Tref)) . (5.22)

In dieser Form können wir zudem mit Kenntnis des Hookeschen Gesetzes der Form

σij = Cijklεelkl (5.23)

die additive Zerlegung des reversiblen Verzerrungsanteils mit

εrevkl = εelkl + εthkl , εthkl = αTδkl (T − Tref) (5.24)

bestimmen. Im thermischen Verzerrungsanteil wurde hierbei die zusätzliche Annahme getrof-fen, dass es nur eines skalaren Wärmeausdehnungskoeffizienten αT bedarf. Das Material wirddemnach als isotrop festgelegt. Auch wird so angenommen, dass thermische Einflüsse keinerleiScherspannungen verursachen. Mit Gleichung (5.22) haben wir so eine Materialgleichung gefun-den, die Spannungen mit der Temperatur und den Dehnungen verknüpft.

Um die in Gleichung (5.11a) definierte Größe pij interpretieren zu können, stellen wir sogenannteMaxwell-Beziehungen auf. Zu diesem Zweck bestimmen wir zunächst das totale Differentialder Gibbsschen freien Energiedichte φ zu

dφ = ∂φ

∂TdT + ∂φ

∂σkldσkl. (5.25)

Mit Gleichung (5.10) können wir nun die partiellen Ableitungen in Gleichung (5.25) zu

∂φ

∂T= −S, (5.26a)

∂φ

∂σkl= −εrevkl (5.26b)

bestimmen. Leiten wir nun Gleichung (5.26a) nach σkl und Gleichung (5.26b) nach T ab, erhalten

36


wir

∂2φ

∂σkl∂T= − ∂S

∂σkl,

∂2φ

∂T∂σkl= −∂ε

revkl

∂T,

(5.27)

sodass wir hiermit nach dem Satz von Schwarz und mit Gleichung (5.11a)

− ∂S

∂σkl= −∂ε

revkl

∂T⇔ pij = αij = αTδij (5.28)

schreiben können. So können wir die unbekannte tensorielle Größe pij durch den messbarenskalaren Ausdehnungskoeffizienten αT ersetzen.

5.2 Aufstellen der Feldgleichung

Um die erhaltenen Ergebnisse zusammenzufassen, schreiben wir Gleichung (5.11a) unter Nut-zung der Gleichungen (5.12) und (5.28) um und erhalten

ρds = CpT

dT + αTdσkk (5.29)

bzw. in Ratennotation

ρs = CpTT + αTσkk. (5.30)

Setzen wir diesen Ausdruck in die Gibbssche Gleichung (5.6) ein und ersetzen dort zudem denAusdruck auf der linken Seite durch die Bilanz der inneren Energie aus Gleichung (2.50), erhaltenwir als Zwischenergebnis

− ∂qi∂xi

+ σklεplkl − CpT − αTT σrr = 0. (5.31)

Es verbleibt die Aufgabe, die Materialgleichungen für den Spannungstensor und den Wärme-fluss, also das Hookesche Gesetz aus Gleichung (5.22) sowie die Fouriersche Wärmeleitungaus Gleichung (5.5) einzuarbeiten. Mit der additiven Zerlegung der Dehnungen nach den Glei-chungen (5.1) und (5.24)

εkl = εrevkl + εplkl = εelkl + εthkl + εplkl (5.32)

und der Definitionsgleichung (2.24) für die Komponenten des Verzerrungstensors können wirzunächst das Hookesche Gesetz zu

σij = Cijkl

(12

(∂uk∂xl

+ ∂ul∂xk

)− αT (T − Tref) δkl − εplkl

)(5.33)

37


umformen. Da wir uns bereits auf isotrope Materialen festgelegt haben, dürfen wir den Steifig-keitstensor C abermals über Gleichung (3.22) durch die Lamé-Konstanten ausdrücken. So wirdaus Gleichung (5.33) unter Berücksichtigung der bereits in Abschnitt 3 ermittelten Spurfreiheitder plastischen Dehnungen εpl und mit der Definition ∆T = T − Tref der Ausdruck

σij = δijλ

(∂uk∂xk− 3αT∆T

)+ 2µ

(12

(∂ui∂xj

+ ∂uj∂xi

)− αT∆Tδij − εplij

). (5.34)

Setzen wir diese Materialgleichung zusammen mit dem Fourierschen Wärmeleitungsgesetz ausGleichung (5.5) in Gleichung (5.31) ein, erhalten wir

κ∂2T

∂xi∂xi− αTT (3λ+ 2µ)

(∂ui∂xi− 3αT T

)+

+ 2µ(

12

(∂ui∂xj

+ ∂uj∂xi

)− εplij

)εplij − CpT = 0.

(5.35)

In dieser Gleichung muss nun noch einer der in Abschnitt 3 gefundenen Zusammenhänge fürdie plastische Dehnrate eingesetzt werden. Für ein sich linear isotrop verfestigendes Materialwird also Gleichung (3.20) verwendet, für ein sich linear kinematisch verfestigendes MaterialGleichung (3.50). Zudem wird die plastische Verzerrung über die Dehnraten durch

εplij =t∫

t=0

εplij dt (5.36)

ausgedrückt.

38

6 Numerische Behandlung


Die numerische Lösung der in Abschnitt 2.3 vorgestellten Bilanzgleichungen erfolgt mit derMethode der finiten Elemente unter Nutzung von frei verfügbaren interagierenden Programm-bausteinen, die im Rahmen des FEniCS Project8 zur Verfügung gestellt und weiterentwickeltwerden. Partielle Differentialgleichungen werden hierfür in Variationsprobleme überführt, welchedann über entsprechende in FEniCS enthaltene Programmkomponenten diskretisiert und mit-tels Finiter-Elemente-Methode gelöst werden. Die diskutierten Bilanzgleichungen werden alsozunächst unter Berücksichtigung von Rand- und Anfangsbedingungen in die schwache Formu-lierung überführt, um ein Variationsproblem zu stellen.

6.1 Schwache Formulierung der Impulsbilanz

Um eine schwache Formulierung zu erhalten, wird die vektorwertige Impulsbilanz der Statikaus Gleichung (2.42) mit einer vektoriellen Testfunktion δui entsprechender Dimension (hierVerschiebung in m) skalar multipliziert:

∂σji∂xj

δui = 0. (6.1)

Der entstehende Ausdruck wird anschließend über das Gebiet Ω, auf dem die Bilanzgleichungerfüllt sein soll, integriert und wir erhalten∫

Ω

∂σji∂xj

δui dV = 0. (6.2)

Auf den Integranden wenden wir die Produktregel an und finden∫Ω

∂

∂xj(σjiδui) dV −

∫Ω

σji∂δui∂xj

dV = 0. (6.3)

Im ersten Term obiger Gleichung identifizieren wir einen Divergenzausdruck, sodass wir mit demGaussschen Integralsatz weiter zu∫

∂Ω

σjiδuinj dA−∫Ω

σji∂δui∂xj

dV = 0 (6.4)

umformen können. Zuletzt können wir im ersten Term den Cauchyschen Zusammenhang ausGleichung (2.37) anwenden und finden

∫∂Ω

tiδui dA−∫Ω

σji∂δui∂xj

dV = 0, (6.5)

8 http://fenicsproject.org/ (letzter Zugriff am 20. Mai 2013)

39

http://fenicsproject.org/


wobei wir mit ti das auf der Oberfläche wirkende Spannungsvektorfeld bezeichnen. Dies ist dieschwache Formulierung der Impulsbilanz der Statik.

6.2 Schwache Formulierung der Energiebilanz

Zur variationellen Formulierung der in Gleichung (5.35) aufgestellten Energiegleichung wird diesezunächst mit einer skalaren normierten Testfunktion δT/T multipliziert, sodass physikalisch alsLeistung interpretierbare Ausdrücke entstehen. Die entstehende Gleichung wird über das GebietΩ integriert, sodass wir den Ausdruck

∫Ω

δT

Tκ∂2T

∂xi∂xi− δTαT (3λ+ 2µ)


)+

+ δT

T2µ(

12

(∂ui∂xj

+ ∂uj∂xi

)− εplij

)εplij −

δT

TCpT

dV = 0

(6.6)

erhalten. Wieder können wir nach analoger Argumentation die zweite Ableitung im ersten Sum-manden obiger Gleichung über die Produktregel und den Gaussschen Integralsatz um eineOrdnung reduzieren:

∫Ω

κ∂2T

∂xi∂xi

δT

TdV =

∫Ω

κ∂

∂xi

(∂T

∂xi

δT

T

)dV −

∫Ω

κ∂T

∂xi

∂

∂xi

(δT

T

)dV

=∫∂Ω

κ∂T

∂xini︸︷︷︸

=hth

δT

TdA−

∫Ω

κ∂T

∂xi

∂

∂xi

(δT

T

)dV.

(6.7)

Hiermit finden wir die schwache Formulierung der Energiegleichung in der Form

∫Ω

− κ ∂T

∂xi

∂

∂xi

(δT

T

)− δTαT (3λ+ 2µ)


)+

+ δT

T2µ(

12

(∂ui∂xj

+ ∂uj∂xi

)− εplij

)εplij −

δT

TCpT

dV +

∫∂Ω

hthδT

TdA = 0.

(6.8)

6.3 Zeitdiskretisierung und Numerik der Materialgleichungen

Sowohl in der Energiegleichung 6.8 als auch in den Materialgleichungen sind Zeitableitungenenthalten, welche diskretisiert werden müssen. Zeitableitungen einer beliebigen Größe φ diskre-tisieren wir mit dem Euler-Verfahren (siehe etwa Press u. a. (1986)) und bezeichnen dabei dieZeitschrittweite mit ∆t. In Formeln schreiben wir

φ = 1∆t

(φ− φ0

), (6.9)

wobei eine hochgestellte 0 den Wert der Größe im letzten Zeitschritt bedeutet.

40


Die Energiegleichung (6.8) lässt sich dann in zeitdiskreter Form nach Multiplikation mit ∆t zurUmwandlung des Leistungsausdrucks in einen Energieausdruck als

∫Ω

−∆tκ ∂T

∂xi

∂

∂xi

(δT

T

)− δTαT (3λ+ 2µ)


)+

+ δT

T2µ(

12

(∂ui∂xj

+ ∂uj∂xi

)−(εpl,0ij + ∆tεplij

))εplij∆t−

δT

TCp(T − T 0

)dV+

+∫∂Ω

∆ththδT

TdA = 0

(6.10)

schreiben. Die plastische Dehnrate wird nicht diskretisiert, weil für diese Materialgleichungenin Ratenform vorliegen: Gleichung (3.20) für ein sich linear isotrop verfestigendes Material undGleichung (3.50) für ein sich linear kinematisch verfestigendes Material. Die dort auftretendenAbleitungen der totalen Dehnung behandeln wir analog und schreiben für diese

εij = 1∆t

(12

(∂ui∂xj

+ ∂uj∂xi

)− ε0

ij

). (6.11)

Der Spannungstensor in der Impulsgleichung (6.5) wird durch eine der Ratengleichungen (3.29)oder (3.52) ausgedrückt, wobei diese zusätzlich um die thermischen Spannungen aus Glei-chung (5.24) ergänzt werden müssen, da die Betrachtung der klassischen Plastizitätstheoriein Abschnitt 3 keine Berücksichtigung der Temperatur einschließt. Um die Ratengleichungenkompakter zu schreiben, wird eine skalare Feldgröße α wie folgt eingeführt:

α =

1 für f < 0,

0 für f ≥ 0,(6.12)

wobei f das verwendete Fließpotential ist. Unter Berücksichtigung des Thermospannungstermsergibt sich für die linear isotrope Verfestigung aus Gleichung (3.29)

σij =(Cijkl − α

CijtuΣtuΣrsCrskl49hσ

2f + ΣmnCmnopΣop

)εkl − CijklαT δklT . (6.13)

Für die linear kinematische Verfestigung findet sich mit Gleichung (3.52) analog

σij =(Cijkl − α



)εkl − CijklαT δklT . (6.14)

Die Nutzung des Euler-Verfahrens aus Gleichung (6.9) liefert die zeitdiskrete Form des Span-nungstensors zum Einsetzen in die Impulsgleichung mit

σij = σ0ij + ∆tσij . (6.15)

Hierbei ist σij durch die Ratengleichung für ein sich linear isotrop verfestigendes Material inGleichung (6.13) oder durch die Ratengleichung für ein sich linear kinematisch verfestigendes

41


Material in Gleichung (6.14) gegeben. Die in diesen Gleichungen auftretenden Zeitableitungenlassen sich dabei analog mit

εij = 1∆t

(12

(∂ui∂xj

+ ∂uj∂xi

)− ε0

ij

),

T = 1∆t(T − T

0)(6.16)

diskretisieren.

6.4 Gekoppelte variationelle Formulierung und Randbedingungen

Die Addition der Gleichungen (6.5) und (6.10) liefert unter Nutzung des in Gleichung (6.15)gegebenen Zusammenhangs die zeitdiskrete schwache Formulierung des gekoppelten Systemsaus Impuls- und Energiebilanz:

∫Ω

−∆tκ ∂T

∂xi

∂

∂xi

(δT

T

)− δTαT (3λ+ 2µ)


)+

+ δT

T2µ(

12

(∂ui∂xj

+ ∂uj∂xi


))εplij∆t−

δT

TCp(T − T 0

)dV−

−∫Ω

(σ0ij + ∆tσij

) ∂δui∂xj

dV+

+∫∂Ω

∆ththδT

TdA+

∫∂Ω

tiδui dA = 0.

(6.17)

Das Variationsproblem in Gleichung (6.17) ist erst durch Randbedingungen vollständig ge-stellt, wobei in dieser Arbeit nur auf Dirichlet- und Neumann-Randbedingungen eingegangenwerden soll. Neumann-Randbedingungen ergeben sich durch Festlegung der beiden Oberflä-chenterme hth und t in Gleichung (6.17) auf entsprechenden Teilgebieten ∂ΩN. Dirichlet-Randbedingungen müssen nicht in das Variationsproblem eingearbeitet werden, sondern wer-den an anderer Stelle im Programm berücksichtigt, in der Weise, als dass der Verschiebungui bzw. der Temperatur T auf entsprechenden Teilgebieten ∂ΩD mit zu setzenden Dirichlet-Randbedingungen direkt der geforderte feste oder zeitlich veränderliche Wert zugewiesen wird.Das bedeutet, wir können den Rand zweckmäßig in zwei Teilgebiete

∂ΩN + ∂ΩD = ∂Ω (6.18)

unterteilen. Auf dem Gebiet ∂ΩD wird explizit die Lösung vorgegeben, sodass die Testfunktionenauf diesen Teilgebieten den Wert Null annehmen. In anderen Worten: Ist die Lösung bekannt,

42


dann muss nicht getestet werden. Damit folgt für die schwache Formulierung

∫Ω

−∆tκ ∂T

∂xi

∂

∂xi

(δT

T

)− δTαT (3λ+ 2µ)


)+

+ δT

T2µ(

12

(∂ui∂xj

+ ∂uj∂xi


))εplij∆t−

δT

TCp(T − T 0

)dV−

−∫Ω

(σ0ij + ∆tσij

) ∂δui∂xj

dV+

+∫

∂ΩN,h

∆ththδT

TdA+

∫∂ΩN,σ

tiδui dA = 0

(6.19)

mit vorzugebendem gerichtetem Wärmefluss über die Oberfläche hth und vorzugebendem Span-nungsvektorfeld auf der Oberfläche t. Die in dieser Arbeit verwendeten Randbedingungen werdenim Zuge der später folgenden Diskussion konkreter Berechnungen vorgestellt. Zur Schließung derProblemstellung müssen letztlich die Materialgleichungen für die in Gleichung (6.17) enthalte-nen Größen εplij , ε

plij und σij für die numerische Auswertung aufbereitet und in geeigneter Weise

eingearbeitet werden.

6.5 Prinzipielle Struktur der Berechnung

Eine weitere zu behandelnde Problematik ergibt sich durch die implizite Struktur der Material-gleichungen (6.13) und (6.14). Hier hängt die Spannungsrate offenbar vom Spannungsdeviatorab und kann demnach nicht explizit und direkt berechnet werden. Im Umgang mit diesem Sach-verhalt können in der numerischen Berechnung unterschiedliche Verfahren verwendet werden. InSimo u. Hughes (1998) oder auch in Doghri (2000) werden beispielsweise verschiedene return-mapping-algorithms bzw. predictor-corrector-Verfahren vorgestellt, welche eine geratene Lösung(predictor) solange anpassen, bis die Gleichungen erfüllt sind (corrector). Diese Verfahren brin-gen jedoch einen großen zusätzlichen Rechenaufwand mit sich, sodass in dieser Arbeit von derVerwendung abgesehen wird. Stattdessen werden für noch unbekannte zur Berechnung benö-tigte Größen die bekannten Werte aus dem letzten Zeitpunkt verwendet, gekennzeichnet durcheine hochgestellte 0. Für genügend kleine Zeitschritte ∆t wird dabei von einem vernachlässigbarkleinen Fehler ausgegangen.

Konkret erhalten wir unter Verwendung der in diesem Abschnitt diskutierten Ergebnisse den

43


folgenden Satz Gleichungen:

Form =∫Ω

−∆tκ ∂T

∂xi

∂

∂xi

(δT

T

)− δTαT (3λ+ 2µ)


)+

+ δT

T2µ(

12

(∂ui∂xj

+ ∂uj∂xi


))εplij∆t−

δT

TCp(T − T 0

)dV−

−∫Ω

(σ0ij + ∆tσij

) ∂δui∂xj

dV +∫

∂ΩN,h

∆ththδT

TdA+

∫∂ΩN,σ

tiδui dA

σij =

(Cijkl − α0 CijtuΣ0

tuΣ0rsCrskl

49hσ

2f + Σ0

mnCmnopΣ0op

)εkl−

−CijklαTδkl1

∆T(T − T 0

),

σf = σf,0 + hεpl,

iso.,

Cijkl − α0 Cijtu(Σ0tu − α0

tu

) (Σ0rs − α0

rs

)Crskl

49hkσ

2f,0 + (Σ0

mn − α0mn)Cmnop

(Σ0op − α0

op

) εkl−

−CijklαTδkl1

∆T(T − T 0

),

αij = α0ij + ∆t

(σij − α0

ij

)Γ, Γ =

(Σ0ij − α0

ij

)σij

23σ

2f,0

,

kin.,

εplkl =

Σ0rsCrstuΣ0

kl

Σ0mnCmnopΣ0

op + 49hσ

2fεtu,

σf = σf,0 + hεpl,

iso.,

εplij =(Σ0rs − αrs

)Crstu

(Σ0ij − αij

)49hkσ

2f,0 + (Σ0

mn − αmn)Cmnop(Σ0op − αop

) εtu,αij = α0

ij + ∆t(σij − α0

ij

)Γ, Γ =

(Σ0ij − α0

ij

)σij

23σ

2f,0

,

kin.,

εkl = 1∆t

(12

(∂uk∂xl

+ ∂ul∂xk

)− ε0

kl

).

(6.20)

44

7 Konkrete numerische Berechnungen


Die im vorigen Kapitel vorgestellten Gleichungen können nun durch Festlegung von Randbe-dingungen und Geometrien auf konkrete Probleme spezifiziert und schließlich mittels finiterElementmethode gelöst werden.

7.1 Vergleich zu den analytischen Lösungen – isotherme Berechnung mitisotropem Verfestigungsmodell

An dieser Stelle sollen die in Abschnitt 4 vorgestellten Modellprobleme mit der Methode derfiniten Elemente nachgerechnet werden, um die korrekte Wiedergabe des Modells durch die nu-merische Berechnung zu überprüfen. Berechnet wurde die Lösung der schwachen Formulierung

(a) Scherklotz unverformt (b) Scherklotz verformt, 50fach vergrößerte Darstel-lung der Verschiebungen

Abbildung 7.1: Scherklotz in der numerischen Berechnung

der Impulsbilanz in drei Raumdimensionen unter verschiedenen Randbedingungen, die die Sys-teme in den Abbildungen 4.1 und 4.3 wiedergeben. Materialkonstanten sind in der Größenord-nung von für gewöhnlichen Stahl zu findenden Werten gewählt worden. Der Scherklotz in Abbil-dung 7.1(a) ist ein Würfel, bei dem die Verschiebungen auf Ober- und Unterseite mittels Dirich-let-Randbedingungen gesetzt werden. Auf der unteren Seite werden diese zu Null gezwungen.Auf der oberen Seite wird die Tangentialkomponente als linear monoton wachsend festgelegt,um eine Scherung zu erreichen. Die restlichen Ränder sind über Neumann-Randbedingungenals spannungsfrei festgelegt. In Abbildung 7.1(b) ist der Endverformungszustand dargestellt.Die Farbgebung stellt den Betrag der Verschiebungsvektoren dar. Anhand dieser lässt sich leichtdie Unabhängigkeit der Verschiebung von zwei der drei Koordinaten erkennen, wie es im semi-inversen Ansatz der analytischen Lösung in Gleichung (4.7) angenommen wurde. Die analytischeLösung für den Spannungs-Dehnungszusammenhang nach den Gleichungen (4.36) und (4.43) istin Abbildung 7.2(a) dargestellt und lässt eine deutliche Übereinstimmung mit der numerischen

45


Lösung9 erkennen. Auf ähnliche Weise wird die Zugprobe in Abbildung 7.3 behandelt. Die Ver-

0 0.0002 0.0004 0.0006 0.0008 0.001 0.0012 0.0014 0.00160

2

4

6

8

10

12

14

16

18x 10

7

σ12/Pa

ε12

AnalytischNumerisch

(a) Scherklotz

0 0.5 1 1.5 2

x 10−3

0

0.5

1

1.5

2

2.5

3x 10

8

σ11/Pa

ε11

AnalytischNumerisch

(b) Zugprobe

Abbildung 7.2: Spannungs-Dehnungs-Diagramme – Analytische und numerische Lösungen.

schiebungen auf der linken Seite werden mit Dirichlet-Randbedingungen zu Null gezwungenund auf der rechten Seite zeitabhängig als monoton ansteigend festgelegt. Die übrigen Rän-der sind wieder als spannungsfrei festgelegt. Wieder ist die Unabhängigkeit der Verschiebungen

(a) Zugprobe unverformt (b) Zugprobe verformt, 200fach vergrößerte Darstel-lung der Verschiebungen

Abbildung 7.3: Zugprobe in der numerischen Berechnung

von zwei der drei Ortskoordinaten zu erkennen. Das Spannungs-Dehnungs-Diagramm in Ab-bildung 7.2(b) stellt die analytische Lösung aus Gleichung (4.54) der numerischen Lösung10

gegenüber und lässt eine optimale Übereinstimmung erkennen.

7.2 Problematik bei kinematischer Verfestigung und thermischer Kopplung

Nach den Betrachtungen im vorigen Teilabschnitt können wir offenbar zumindest von einer hin-reichenden qualitativen sowie quantitativen Korrektheit der Lösung der Impulsbilanz mittelsfiniter Elementmethode für ein sich linear isotrop verfestigendes plastifizierendes Material aus-gehen. Bei dem Versuch, auf analoge Weise das kinematische Verfestigungsmodell – zunächst9 Vgl. Anhang A.110Vgl. Anhang A.2

46


unter Vernachlässigung thermischer Kopplung – numerisch zu behandeln, ist die zur Verfügungstehende Rechenkapizität11 schon bei der Erstellung der linearen Gleichungssysteme für die Ele-mente an ihre Grenzen geraten, sodass sich der Rechner schon vor einer konkreten Berechnungnicht mehr bedienen ließ und neugestartet werden musste. Erst die Reduktion des Gitters vondrei auf zwei Raumdimensionen führte zur erfolgreichen Berechnung. Jedoch muss dabei berück-sichtigt werden, dass bei der Reduktion der Dimension in jedem Fall Fehler einfließen, da diebehandelten Materialgesetze nur in drei Raumdimensionen die volle in der Herleitung gegebenephysikalische Gültigkeit besitzen.

Nachfolgend soll daher zunächst das kinematische Verfestigungsmodell dem bereits untersuchtenisotropen Verfestigungsmodell in der zweidimensionalen numerischen Behandlung gegenüberge-stellt werden und Fehlereinflüsse gesucht werden.

7.3 Zugprobe unter wechselseitiger Belastung in zwei Dimensionen – Verfes-tigungsmodelle im Vergleich

Als Beispielproblem dient die in Abbildung 7.5 dargestellte linksseitig festgehaltene und rechts-seitig weggesteuert gedehnte bzw. gestauchte Zugprobe. Der Zug bzw. Druck wird mit Dirich-let-Randbedingung am rechten Rand über eine zeitabhängige Verschiebungsfunktion ur(t) =ur1(t)e1 nach Abbildung 7.4 vorgegeben. Diese Wahl der Funktion ermöglicht die Beurteilung

0 0.2 0.4 0.6 0.8 1−1

0

1

ttEnde

ur 1(t)

uM

ax

Abbildung 7.4: Dirichlet-Randbedingung am rechten Rand – ur1(t)

der qualitativen Korrektheit in der Modellierung des in Abschnitt 3.1.5 diskutierten Bauschin-ger-Effekts und somit auch die Beurteilung der qualitativen Korrektheit bei der Berechnungmit dem kinematischen Verfestigungsmodell. In diesem Fall wird abermals nur die Lösung derImpulsbilanz bestimmt, also ein isothermer Prozess angenommen, da das Hauptaugenmerk hierauf den Verfestigungsmodellen liegt. Neben dieser Randbedingung sind zur Zugprobe im vori-gen Abschnitt analoge Randbedingungen gesetzt worden (siehe Abbildung 7.5). Die Spannungs-Dehnungs-Diagramme in Abbildung 7.6 stellen die Zugspannung σ11 in der Mitte der Probe11Intel R© CoreTM i3-3110M CPU @ 2.40GHz, 3742 MB DDR3 RAM

47


Abbildung 7.5: Rechengebiet mit verwendeten Randbedingungen

über der zugehörigen Verzerrungskomponente ε11 in diesem Punkt dar. In Abbildung 7.6(a) ist

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2

x 10−3

−3

−2

−1

0

1

2

3

4x 10

8

ε11

σ11/Pa

(a) Kinematisches Verfestigungsmodell

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2

x 10−3

−4

−3

−2

−1

0

1

2

3

4x 10

8

ε11

σ11/Pa

(b) Isotropes Verfestigungsmodell

Abbildung 7.6: Zugprobe unter wechselseitiger Belastung, Fließen in beiden Belastungsrichtun-gen – Spannungs-Dehnungs-Diagramme für beide Verfestigungsmodelle

die Simulation des in Abschnitt 3.1.5 beschriebenen Bauschinger-Effekts durch das Modell derkinematischen Verfestigung, also die Verminderung der Initialfließspannung bei Belastung in Ge-genrichtung, umittelbar zu erkennen. In Abbildung 7.6(b) ist das unter denselben Bedingungenmit dem isotropen Verfestigungsmodell berechnete Spannungs-Dehnungsdiagramm zu sehen, beiwelchem die Fließgrenze stets wächst, der erste Fließvorgang jedoch dem mit dem kinematischenVerfestigungsmodell berechneten identisch ist. Die qualitative Güte der Ergebnisse ist offenbarals hinreichend anzunehmen, insbesondere bei Vergleich mit den in Abbildung 3.3 dargestell-ten Phänomenen. Quantitativ betrachtet muss festgestellt werden, dass bereits die Steigung derFunktionen im linear-elastischen Bereich eine leichte Abweichung vom eingestellten Elastizitäts-modul zeigt. Da uns jedoch keine Alternative geboten ist, werden die folgenden, die thermischenBetrachtungen einschließenden Berechnungen vornehmlich auf qualitativer Ebene diskutiert undquantitativ nur im Rahmen von Größenordnungen beurteilt werden.

7.4 Zugprobe unter wechselseitiger Belastung in zwei Dimensionen – Ther-momechanische Kopplung

Betrachtet wird abermals die unter der in 7.4 vorgegebenen Belastung in Abbildung 7.5 dar-gestellte Zugprobe, für die in diesem Beispiel die vollständige Gleichung (6.20) gelöst werden

48


soll. Die thermische Kopplung wird also unter Einbezug der Bilanz der inneren Energie be-rücksichtigt, sodass das Temperaturfeld als zusätzliche gesuchte Größe erscheint12. Verwendetwird das kinematische Verfestigungsmodell, da das isotrope Verfestigungsmodell bei wechselsei-tigem Fließen physikalisch falsche Ergebnisse liefert. Die Randbedingungen für die Tempera-tur werden als Neumann-Randbedingungen formuliert, indem über einen Wärmefluss Energieüber den Rand mit der Umgebung ausgetauscht werden kann. Die Anfangstemperatur wird beiRaumtemperatur von T0 = 293,15 K festgelegt. In Abbildung 7.7 ist neben dem Spannungs-Dehnungs-Diagramm die Temperaturentwicklung über der Zeit aufgetragen. Bei Vergleich mit

−1.5 −1 −0.5 0 0.5 1 1.5 2

x 10−3

−3

−2

−1

0

1

2

3x 10

8

σ11/Pa

ε110 0.2 0.4 0.6 0.8 1

292.9

293

293.1

293.2

293.3

293.4

293.5

293.6

T/K

ttEnde

Abbildung 7.7: Spannungs-Dehnungs-Diagramm und Temperaturentwicklung der wechselseitigbelasteten Zugprobe

Abbildung 7.6(a) ist zunächst festzustellen, dass der Einbezug der Temperaturentwicklung nureinen insignifikanten Einfluss auf den Verlauf der Spannungen über den Dehnungen hat. Zudemkönnen Abweichungen auch aufgrund der notwendigerweise für die Berechnung der gekoppeltenGleichungen stark verringerten Elementanzahl im Gitter zustande kommen. In Anbetracht deskleinen Maßstabes, in dem sich die Temperaturentwicklung abspielt – Änderungen finden sicherst in der ersten Nachkommastelle – scheint diese Tatsache auch sehr plausibel. Die Tempera-turentwicklung für sich zeigt, dass sich die Probe bei Zug abkühlt. Dieses Phänomen kann durchdie im elastischen Bereich entropieabführende Wirkung der Volumenkontraktion begründet wer-den (Abali 2011). Die zeitliche Rate der Abkühlung nimmt nach dem Einsetzen des plastischenFließens bei etwa t/tEnde = 0,17 auf fast Null ab. Die Entropie nimmt somit nicht mehr im selbenMaße wie im elastischen Bereich ab, ein Indiz für die vorhandene Dissipation, wenn auch in sehrkleinen Größenordnungen. In Abbildung 7.8 ist das Temperaturfeld für t/tEnde = 0,2 dargestellt.Die Homogenität ist gut erkennbar. Ebenso lässt sich sehr gut feststellen, dass für das Beispieldes plastisch fließenden Zugstabes die Annahme der Isothermie gerechtfertigt ist.

12Vgl. Anhang A.5

49


Abbildung 7.8: Temperaturverteilung in der Probe zur Zeit t/tEnde = 0,2

7.5 Scherklotz unter wechselseitiger Belastung in zwei Dimensionen – Ther-momechanische Kopplung

Neben der Zugprobe soll auch der Scherklotz unter dem Aspekt der thermomechanischen Kopp-lung untersucht werden. Das gekoppelte Gleichungssystem entspricht dem in Abschnitt 7.4 fürdie Zugprobe genutzten System aus Gleichung (6.20). Die Randbedingungen haben wir analogzu den in Abschnitt 7.1 gesetzten Randbedingungen gewählt.13 In Abbildung 7.9 ist links dieSpannungskomponente σ12 über der Dehnungskomponente ε12 und rechts der Temperaturver-lauf über der Zeit dargestellt. An diesem Beispiel können deutlich die Fehlereinflüsse durch die

−1.5 −1 −0.5 0 0.5 1 1.5

x 10−3

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2x 10

8

σ12/Pa

ε120 0.2 0.4 0.6 0.8 1

293.12

293.14

293.16

293.18

293.2

293.22

293.24

293.26

T/K

ttEnde

Abbildung 7.9: Spannungs-Dehnungs-Diagramm und Temperaturentwicklung des wechselseitigbelasteten Scherklotzes

kapazitätsbedingt reduzierte numerische Genauigkeit erkannt werden. Beispielsweise entwickelnsich im Temperaturverlauf kleine Spitzen zu den Zeitpunkten, zu denen das plastische Fließeneinsetzt. Dies ist durch die Nutzung der Werte aus dem letzten Zeitschritt bedingt. Qualitativkönnen wir aber gleichermaßen das erwartete Ergebnis erkennen. Die Temperatur steigt durchdie plastische Verformung im Scherklotz stetig und kann nicht durch elastische Entropieab-fuhr reguliert werden, da die Volumenänderung im elastischen Bereich des Scherexperiments13Vgl. auch Anhang A.6

50


sehr klein ist. In Abbildung 7.10 sind die zu den betragsmäßig größten Scherungen zugehörigenVerschiebungs- und Temperaturfelder dargestellt. Die annähernde Konstanz des Volumens wirddeutlich sowie die leicht steigende Temperatur in der gesamten Probe.

Abbildung 7.10: Ebener Scherklotz – numerische Lösungen für Verschiebungen und Temperatur:Verformungen in 50fach vergrößerter Darstellung (links) und Temperaturfeld(rechts) zu den Zeiten t/tEnde = 0,2 (oben) und t/tEnde = 0,6 (unten). Die Ein-färbung der linken Abbildungen stellt die Beträge der Verschiebungsvektorendar.

51

8 Resümee

8 Resümee

Ziel dieser Arbeit war es, die rein mechanische klassische Plastizitätstheorie um thermischeEinflüsse zu ergänzen und mithilfe der Methode der finiten Elemente Temperatur- und Ver-schiebungsfelder zu berechnen. Wir haben zunächst die rein mechanische Theorie untersuchtund analytische Lösungen für ausgewählte nicht dreidimensionale Beispiele gefunden. Innerhalbder rein mechanischen Theorie gelang es, den allgemeinen dreidimensionalen Fall unter monoto-ner Belastung numerisch zu berechnen und die Resultate mit den analytischen zu vergleichen,was sehr gute Ergebnisse lieferte. Auf Grundlagen der kontinuumstheoretischen Thermodynamikaufbauend, konnten die klassischen Gleichungen der Plastizitätstheorie um thermische Effekteergänzt werden. Durch die stark nichtlineare und nicht-explizite Stuktur der Gleichungen imkomplizierteren kinematischen Verfestigungsmodell mit thermomechanischer Kopplung, gelanges durch begrenzte Rechenleistung leider nur, eine qualitative Untersuchung anhand der Lösun-gen auf zweidimensionalen Gittern vorzunehmen, welche jedoch im direkten Vergleich mit dendreidimensionalen Lösungen plausible Ergebnisse lieferten. Insofern betrachten wir die gesetzteZielstellung als erreicht, wobei aufgrund limitierter technischer Möglichkeiten die Umsetzungnicht im vollen Maße möglich war.

52

A Quelltexte

A Quelltexte

A.1 Scherklotz – 3D, isotherm, isotrope Verfestigung

1 from d o l f i n import ∗2 import numpy as np3 from numpy import array4 import matp lo t l i b . pyplot as p l5

6 parameters [ " form_compiler " ] [ " cpp_optimize " ] = True7 parameters [ " form_compiler " ] [ " opt imize " ] = True8

9 " " "−−−−−−−−−−−−newton s o l v e r i n t e r a c t i o n−−−−−−−−−−−−−−" " "10 class i t e r a t e ( NonlinearProblem ) :11 def __init__( s e l f , a , L , bc , exter_B ) :12 NonlinearProblem . __init__( s e l f )13 s e l f . L = L14 s e l f . a = a15 s e l f . bc = bc16 s e l f . ex t e r = exter_B17 def F( s e l f , b , x ) :18 assemble ( s e l f . L , t en so r=b , exter ior_facet_domains=s e l f . ex t e r )19 for cond i t i on in s e l f . bc : cond i t i on . apply (b , x )20 def J ( s e l f , A, x ) :21 assemble ( s e l f . a , t en so r=A, exter ior_facet_domains=s e l f . ex t e r )22 for cond i t i on in s e l f . bc : cond i t i on . apply (A)23

24 " " "−−−−−−−−−−−mesh−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "25 xMin = .026 xMax = .627 xElements = 1228 i , j , k , l , r , m, o , p , s , n , t , q= i nd i c e s (12)29 yMin = .030 yMax = .631 yElements = 1232 zMin = .033 zMax = .634 zElements = 235 mesh = BoxMesh(xMin , yMin , zMin , xMax , yMax , zMax , xElements , yElements , zElements )36

37 " " "−−−−−−−−−f unc t i on spaces−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "38 Coef f = FunctionSpace (mesh , ’CG’ , 1)39 Vector = VectorFunctionSpace (mesh , ’CG’ , 1)40 Tensor = TensorFunctionSpace (mesh , ’CG’ , 1)41

42 " " "−−−−−−−−−−func t ions−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "43 ur e f = Function ( Vector )44 u = Function ( Vector )45 u0 = Function ( Vector )46

47 du = Tria lFunct ion ( Vector )48 delu = TestFunction ( Vector )49

50 " " "−−−−−−−−−− i n i t i a l condi t ions−−−−−−−−−−−−−−−−−−−−−−−−" " "51 class I n i t i a lC ond i t i o n s ( Express ion ) :52 def eva l ( s e l f , vec , x ) :53 vec [ 0 ] = . 0 #u[0]=0

53

A Quelltexte

54 vec [ 1 ] = . 0 #u[1]=055 vec [ 2 ] = . 056 def value_shape ( s e l f ) :57 return ( 3 , )58

59 in iObj = I n i t i a lC ond i t i o n s ( )60 ur e f . i n t e r p o l a t e ( in iObj )61 u . a s s i gn ( u r e f )62

63 " " "−−−−−−−−−−−−boundary condi t ions−−−−−−−−−−−−−−−−−−−−−" " "64

65 top = compile_subdomains ( ’ x [ 1 ] ==Max ’ )66 top .Max= yMax67

68 bot = compile_subdomains ( ’ x [ 1 ] ==Min ’ )69 bot .Min= yMin70

71 i n i tD i s p = Express ion ( ( ’ 0 . 0 ’ , ’ 0 . 0 ’ , ’ 0 . 0 ’ ) )72 bc1 = Dir ichletBC ( Vector , i n i tD i sp , bot )73 shear = Express ion ( ( ’A∗ t ’ , ’ 0 ’ , ’ 0 ’ ) ,A=xMax/350 . , t =0.)74 bc2 = Dir ichletBC ( Vector , shear , top )75 boundaryConditions = [ bc1 , bc2 ]76

77 subdomains = MeshFunction ( " s i z e_t " , mesh , 2)78

79 #−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−80 de l t a = Id en t i t y (3 )81

82 " " "−−−−−−−−−−time−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "83 time =.084 dt = .02585 tMax = 1 .586

87 " " "−−−−−−−−−−mater ia l:−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "88 nu , E = 0 . 3 , 200 . e+0989 E_T = .2∗E90 h=E_T/(1.−E_T/E)91 #h= .1∗E#.4∗E92 lamda = E∗nu/(1.+nu) /(1.−2.∗nu)93 mu = E/(2.+ 2 .∗nu)94

95 Y=Function ( Coe f f )96 Y0 = Constant ( 250 . e+6)97 Y. vec to r ( ) [ : ] = 250 . e+698

99 " " "−−−−−−−−−Strain−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "100 eps0=Function ( Tensor )101 eps = as_tensor ( . 5 ∗ ( u [ i ] . dx ( j )+u [ j ] . dx ( i ) ) , [ i , j ] )102 deps=as_tensor ( eps [ k , l ]− eps0 [ k , l ] , ( k , l ) )103

104 Auf " " "−−−−−−−−−Stress−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "105 a l f a=Function ( Coe f f )106 cauchy0 = Function ( Tensor )107 sdev0=as_tensor ( cauchy0 [ i , j ] −1 ./3 .∗ cauchy0 [ k , k ]∗ de l t a [ i , j ] , [ i , j ] )108 SVM0 = ( 3 . / 2 . ∗ ( sdev0 [ i , k ]∗ sdev0 [ i , k ]+0.00001) ) ∗∗ . 5109 dcauchy = as_tensor (E/(1.+nu) ∗( deps [ i , j ]+nu/(1.−2.∗nu) ∗deps [ k , k ]∗ de l t a [ i , j ]− a l f a ∗ ( 3 .∗E)

/ (3 .∗E+2.∗(1.+nu) ∗h) ∗3 . / 2 .∗ sdev0 [ k , l ]∗ deps [ k , l ] /Y∗ sdev0 [ i , j ] /Y) , ( i , j ) )110 cauchy = as_tensor ( cauchy0 [ i , j ]+dcauchy [ i , j ] , ( i , j ) )

54

A Quelltexte

111 depspl = as_tensor ( a l f a ∗1 ./h ∗3 . / 2 .∗ sdev0 [ i , j ]∗ dcauchy [ i , j ] /Y∗3 . / 2 .∗ sdev0 [ i , j ] /Y, ( i , j ) )112 epsacc=Function ( Coe f f )113

114 " " "−−−−−−−−−−−−−−−−−−−" " "115 Form = cauchy [ i , r ]∗ delu [ i ] . dx ( r ) ∗dx116 Gain = de r i v a t i v e (Form , u , du)117

118 f i l eU = F i l e ( ’ . / Data/U. pvd ’ )119

120 pl . ion ( )121 pl . c l f ( )122 pl . c l a ( )123 pl . x l ab e l ( ’ eps_12 ’ )124 pl . y l ab e l ( ’ sigma_12 ’ )125 pl . t i t l e ( ’ S che rk l o t z ’ )126 pl . g r i d (True )127 p lSt ra in , p l S t r e s s = [ ] , [ ]128 while time <= tMax :129 shear . t=time130 problem = i t e r a t e (Gain , Form , boundaryConditions , subdomains )131 s o l v e r = NewtonSolver ( ’ lu ’ )132 s o l v e r . parameters [ ’ c onve rgence_cr i t e r i on ’ ] = ’ incrementa l ’133 s o l v e r . parameters [ ’ r e l a t i v e_ t o l e r an c e ’ ] = 1 .0 e−3134 s o l v e r . s o l v e ( problem , u . vec to r ( ) )135 cauchyProj=p ro j e c t ( cauchy , Tensor )136 cauchy0 . a s s i gn ( cauchyProj )137 p lS t r a i n . append ( p r o j e c t ( eps , Tensor ) ( (xMax/2 . , 3 . ∗yMax/4 . , zMax/2 . ) ) [ 1 ] )138 p l S t r e s s . append ( cauchyProj ( (xMax/2 . , 3 . ∗yMax/4 . , zMax/2 . ) ) [ 1 ] )139 pl . p l o t ( p lS t ra in , p l S t r e s s , l i n e s t y l e=’− ’ , c o l o r=’ red ’ , l i n ew id th=1, marker=’ ’ ,

markers i ze=11)140 pl . draw ( )141 pl . s a v e f i g ( ’ . / s c h e r k l o t z . eps ’ )142 d i r e c tP r o j = p ro j e c t ( sdev0 [ k , l ]∗ deps [ k , l ] , Coe f f )143 d i r e c t a r r = d i r e c tP r o j . vec to r ( ) . array ( )144 d i r e c t b o o l = np . s i gn ( . 5 ∗ ( d i r e c t a r r +(( d i r e c t a r r ) ∗∗2 . ) ∗∗ . 5 ) )145 epsaccPro j=p ro j e c t ( epsacc+a l f a ∗ ( 2 . / 3 . ∗ depspl [ i , j ]∗ depspl [ i , j ] ) ∗∗ . 5 , Coe f f )146 Y. vec to r ( ) [ : ]= 2 5 0 . e+6+h∗ epsaccPro j . vec to r ( ) . array ( ) [ : ]147 SVM0proj=p ro j e c t (SVM0, Coe f f )148 y i e l dboo l = SVM0proj . vec to r ( ) . array ( ) >= Y. vec to r ( ) . array ( )149 a l f a . vec to r ( ) [ : ] = np . array ( y i e l dboo l ∗ d i r e c tboo l , dtype=in t )150 eps0 . a s s i gn ( p r o j e c t ( eps , Tensor ) )151 epsacc . a s s i gn ( epsaccPro j )152 Yproj=p ro j e c t (Y, Coe f f )153 f i l eU << (u , time )154 time+=dt

A.2 Zugprobe – 3D, isotherm, isotrope Verfestigung



55

A Quelltexte


24 " " "−−−−−−−−−−−mesh−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "25 xMin = .026 xMax = 1 .27 xElements = 2028 i , j , k , l , r , m, o , p , s , n , t , q= i nd i c e s (12)29 yMin = .030 yMax = .231 yElements = 432 zMin = .033 zMax = .134 zElements = 235 mesh = BoxMesh(xMin , yMin , zMin , xMax , yMax , zMax , xElements , yElements , zElements )36

37 " " "−−−−−−−−−f unc t i on spaces−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "38 Coef f = FunctionSpace (mesh , ’CG’ , 1)39 Vector = VectorFunctionSpace (mesh , ’CG’ , 1)40 Tensor = TensorFunctionSpace (mesh , ’CG’ , 1)41



50 " " "−−−−−−−−−− i n i t i a l condi t ions−−−−−−−−−−−−−−−−−−−−−−−−" " "51 class I n i t i a lC ond i t i o n s ( Express ion ) :52 def eva l ( s e l f , vec , x ) :53 vec [ 0 ] = . 0 #u[0]=054 vec [ 1 ] = . 0 #u[1]=055 vec [ 2 ] = . 056 def value_shape ( s e l f ) :57 return ( 3 , )58


63 " " "−−−−−−−−−−−−boundary condi t ions−−−−−−−−−−−−−−−−−−−−−" " "64

65 r i g h t = compile_subdomains ( ’ x [ 0 ] ==Max ’ )66 r i g h t .Max = xMax

56

A Quelltexte

67

68 l e f t = compile_subdomains ( ’ x [ 0 ] ==Min ’ )69 l e f t .Min= xMin70 l e f t .Max = xMax71

72 f r on t = compile_subdomains ( ’ x [ 1 ] ==Max ’ )73 f r on t .Max = yMax74

75 back = compile_subdomains ( ’ x [ 1 ] ==Min ’ )76 back .Min= yMin77

78 i n i tD i s p = Express ion ( ( ’ 0 . 0 ’ , ’ 0 . 0 ’ , ’ 0 . 0 ’ ) )79 bc1 = Dir ichletBC ( Vector , i n i tD i sp , l e f t )80 t e n s i l e = Express ion ( ( ’A∗ t ’ , ’ 0 ’ , ’ 0 ’ ) ,A=xMax/400 . , t =0.)81 bc2 = Dir ichletBC ( Vector , t e n s i l e , r i g h t )82 bc3 = Dir ichletBC ( Vector . sub (1 ) , ’ 0 ’ , r i g h t )83 bc4 = Dir ichletBC ( Vector . sub (2 ) , ’ 0 ’ , r i g h t )84

85 boundaryConditions = [ bc1 , bc2 , bc3 , bc4 ]86


89 #−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−90 de l t a = Id en t i t y (3 )91

92 " " "−−−−−−−−−−time−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "93 time =.094 dt = .00595 tMax = 2 .596

97 " " "−−−−−−−−−−mater ia l:−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "98 nu , E = 0 . 3 , 200 . e+0999 E_T = .2∗E

100 h=E_T/(1.−E_T/E)101 lamda = E∗nu/(1.+nu) /(1.−2.∗nu)102 mu = E/(2.+ 2 .∗nu)103


108 " " "−−−−−−−−−Strain−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "109 eps0=Function ( Tensor )110 eps = as_tensor ( . 5 ∗ ( u [ i ] . dx ( j )+u [ j ] . dx ( i ) ) , [ i , j ] )111 deps=as_tensor ( eps [ k , l ]− eps0 [ k , l ] , ( k , l ) )112

113 " " "−−−−−−−−−Stress−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "114 a l f a=Function ( Coe f f )115 cauchy0 = Function ( Tensor )116 sdev0 = as_tensor ( cauchy0 [ i , j ] −1 ./3 .∗ cauchy0 [ k , k ]∗ de l t a [ i , j ] , [ i , j ] )117 SVM0 = ( 3 . / 2 . ∗ ( sdev0 [ i , k ]∗ sdev0 [ i , k ]+0.00001) ) ∗∗ . 5118 dcauchy = as_tensor (E/(1.+nu) ∗( deps [ i , j ]+nu/(1.−2.∗nu) ∗deps [ k , k ]∗ de l t a [ i , j ]− a l f a ∗ ( 3 .∗E)

/ (3 .∗E+2.∗(1.+nu) ∗h) ∗3 . / 2 .∗ sdev0 [ k , l ]∗ deps [ k , l ] /Y∗ sdev0 [ i , j ] /Y) , ( i , j ) )119 cauchy = as_tensor ( cauchy0 [ i , j ]+dcauchy [ i , j ] , ( i , j ) )120 depspl = as_tensor ( a l f a ∗1 ./h ∗3 . / 2 .∗ sdev0 [ i , j ]∗ dcauchy [ i , j ] /Y∗3 . / 2 .∗ sdev0 [ i , j ] /Y, ( i , j ) )121 epsacc = Function ( Coe f f )122

123 " " "−−−−−−−−−−−−−−−−−−−" " "

57

A Quelltexte

124 Form = cauchy [ i , r ]∗ delu [ i ] . dx ( r ) ∗dx125 Gain = de r i v a t i v e (Form , u , du)126

127 f i l eU = F i l e ( ’ . / Data/U. pvd ’ )128

129 pl . ion ( )130 pl . c l f ( )131 pl . c l a ( )132 pl . x l ab e l ( ’ eps_11 ’ )133 pl . y l ab e l ( ’ sigma_11 ’ )134 pl . t i t l e ( ’ Zugstab ’ )135 pl . g r i d (True )136 p lSt ra in , p l S t r e s s = [ ] , [ ]137 while time <= tMax :138 t e n s i l e . t=time139 problem = i t e r a t e (Gain , Form , boundaryConditions , subdomains )140 s o l v e r = NewtonSolver ( ’ lu ’ )141 s o l v e r . parameters [ ’ c onve rgence_cr i t e r i on ’ ] = ’ incrementa l ’142 s o l v e r . parameters [ ’ r e l a t i v e_ t o l e r an c e ’ ] = 1 .0 e−3143 s o l v e r . s o l v e ( problem , u . vec to r ( ) )144 cauchyProj=p ro j e c t ( cauchy , Tensor )145 cauchy0 . a s s i gn ( cauchyProj )146 p lS t r a i n . append ( p r o j e c t ( eps , Tensor ) ( (xMax/2 . , yMax/2 . , zMax/2 . ) ) [ 0 ] / xMax)147 p l S t r e s s . append ( cauchyProj ( (xMax/2 . ,yMax/2 . , zMax/2 . ) ) [ 0 ] )148 pl . p l o t ( p lS t ra in , p l S t r e s s , l i n e s t y l e=’− ’ , c o l o r=’ red ’ , l i n ew id th=1, marker=’ ’ ,

markers i ze=11)149 pl . draw ( )150 pl . s a v e f i g ( ’ . / webs i t e r e . eps ’ )151 d i r e c tP r o j = p ro j e c t ( sdev0 [ k , l ]∗ deps [ k , l ] , Coe f f )152 d i r e c t a r r = d i r e c tP r o j . vec to r ( ) . array ( )153 d i r e c t b o o l = np . s i gn ( . 5 ∗ ( d i r e c t a r r +(( d i r e c t a r r ) ∗∗2 . ) ∗∗ . 5 ) )154 epsaccPro j=p ro j e c t ( epsacc+a l f a ∗ ( 2 . / 3 . ∗ depspl [ i , j ]∗ depspl [ i , j ] ) ∗∗ . 5 , Coe f f )155 Y. vec to r ( ) [ : ]= 2 5 0 . e+6+h∗ epsaccPro j . vec to r ( ) . array ( ) [ : ]156 SVM0proj=p ro j e c t (SVM0, Coe f f )157 y i e l dboo l = SVM0proj . vec to r ( ) . array ( ) >= Y. vec to r ( ) . array ( )158 a l f a . vec to r ( ) [ : ] = np . array ( y i e l dboo l ∗ d i r e c tboo l , dtype=in t )159 eps0 . a s s i gn ( p r o j e c t ( eps , Tensor ) )160 epsacc . a s s i gn ( epsaccPro j )161 Yproj=p ro j e c t (Y, Coe f f )162 f i l eU << (u , time )163 time+=dt

A.3 Zugprobe – 2D, unter wechselseitiger Belastung – isotherm, isotropeVerfestigung



9 " " "−−−−−−−−−−−−newton s o l v e r i n t e r a c t i o n−−−−−−−−−−−−−−" " "10 class i t e r a t e ( NonlinearProblem ) :

58

A Quelltexte

11 def __init__( s e l f , a , L , bc , exter_B ) :12 NonlinearProblem . __init__( s e l f )13 s e l f . L = L14 s e l f . a = a15 s e l f . bc = bc16 s e l f . ex t e r = exter_B17 def F( s e l f , b , x ) :18 assemble ( s e l f . L , t en so r=b , exter ior_facet_domains=s e l f . ex t e r )19 for cond i t i on in s e l f . bc : cond i t i on . apply (b , x )20 def J ( s e l f , A, x ) :21 assemble ( s e l f . a , t en so r=A, exter ior_facet_domains=s e l f . ex t e r )22 for cond i t i on in s e l f . bc : cond i t i on . apply (A)23

24 " " "−−−−−−−−−−−mesh−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "25 xMin = .026 xMax = 1 .27 xElements = 3028 i , j , k , l , r , m, o , p , s = i nd i c e s (9 )29 yMin = .030 yMax = .231 yElements = 632

33 mesh = RectangleMesh (xMin , yMin , xMax , yMax , xElements , yElements )34

35 " " "−−−−−−−−−f unc t i on spaces−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "36 Coef f = FunctionSpace (mesh , ’CG’ , 1)37 Vector = VectorFunctionSpace (mesh , ’CG’ , 1)38 Tensor = TensorFunctionSpace (mesh , ’CG’ , 1)39



48

49 " " "−−−−−−−−−−−−−−−−−−−−−−−−− i n i t i a l cond−−−−−−−−−−−−−−−−−−−−−−−−−" " "50 class I n i t i a lC ond i t i o n s ( Express ion ) :51 def eva l ( s e l f , vec , x ) :52 vec [ 0 ] = . 053 vec [ 1 ] = . 054 def value_shape ( s e l f ) :55 return ( 2 , )56


61 " " "−−−−−−−−−−−−−−−−−−−−−−−−−boundary cond−−−−−−−−−−−−−−" " "62 r i g h t = compile_subdomains ( ’ x [ 0 ] >= .9∗Max ’ )63 r i g h t .Max = xMax64

65 l e f t = compile_subdomains ( ’ x [ 0 ] <=Min+0.1∗Max ’ )66 l e f t .Min= xMin67 l e f t .Max = xMax68

59

A Quelltexte



75

76 i n i tD i s p = Express ion ( ( ’ 0 . 0 ’ , ’ 0 . 0 ’ ) )77 bc1 = Dir ichletBC ( Vector , i n i tD i sp , l e f t )78 t e n s i l e = Express ion ( ( ’A∗ t ’ , ’ 0 ’ ) ,A=xMax/400 . , t =0.)79 bc2 = Dir ichletBC ( Vector , t e n s i l e , r i g h t )80 bc3 = Dir ichletBC ( Vector . sub (1 ) , ’ 0 ’ , r i g h t )81

82 boundaryConditions = [ bc1 , bc2 , bc3 ]83


86 #−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−87 de l t a = Id en t i t y (2 )88

89 " " "−−−−−−−−−−−−−time−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "90 time =.091 dt = .00592 tMax = 2 .593

94 " " "−−−−−−−−−−mater ia l:−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "95 nu , E = 0 . 3 , 200 . e+0996 E_T = .2∗E97 h=E_T/(1.−E_T/E)98 lamda = E∗nu/(1.+nu) /(1.−2.∗nu)99 mu = E/(2.+ 2 .∗nu)

100


105 " " "−−−−−−−−−Strain−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "106 eps0=Function ( Tensor )107 eps = as_tensor ( . 5 ∗ ( u [ i ] . dx ( j )+u [ j ] . dx ( i ) ) , [ i , j ] )108 deps=as_tensor ( eps [ k , l ]− eps0 [ k , l ] , ( k , l ) )109 epsacc=Function ( Coe f f )110

111 " " "−−−−−−−−−Stress−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "112 a l f a=Function ( Coe f f )113 cauchy0 = Function ( Tensor )114 sdev0=as_tensor ( cauchy0 [ i , j ] −1 ./3 .∗ cauchy0 [ k , k ]∗ de l t a [ i , j ] , [ i , j ] )115 #von−mises e q u i v a l e n t s t r e s s i n c l u d i n g a c o r r e c t i o n f o r 3D −> 2D116 SVM0 = ( 3 . / 2 . ∗ ( sdev0 [ i , k ]∗ sdev0 [ i , k ]+0.00001+( cauchy0 [ i , i ] / 3 . ) ∗∗2 . ) ) ∗∗ . 5117 dcauchy=as_tensor (E/(1.+nu) ∗( deps [ i , j ]+nu/(1.−2.∗nu) ∗deps [ k , k ]∗ de l t a [ i , j ]− a l f a ∗ ( 3 .∗E)

/ (3 .∗E+2.∗(1.+nu) ∗h) ∗3 . / 2 .∗ sdev0 [ k , l ]∗ deps [ k , l ] /Y∗ sdev0 [ i , j ] /Y) , ( i , j ) )118 cauchy=as_tensor ( cauchy0 [ i , j ]+dcauchy [ i , j ] , ( i , j ) )119 depspl=as_tensor ( a l f a ∗1 ./h ∗3 . / 2 .∗ sdev0 [ i , j ]∗ dcauchy [ i , j ] /Y∗3 . / 2 .∗ sdev0 [ i , j ] /Y, ( i , j ) )120

121 " " "−−−−−−−−−−−−−−−−−−−" " "122 Form = cauchy [ i , r ]∗ delu [ i ] . dx ( r ) ∗dx123 Gain = de r i v a t i v e (Form , u , du)124

125 f i l eU = F i l e ( ’ . / Data/U. pvd ’ )

60

A Quelltexte

126

127 pl . ion ( )128 pl . c l f ( )129 pl . c l a ( )130 pl . x l ab e l ( ’ eps_11 ’ )131 pl . y l ab e l ( ’ sigma_11 ’ )132 pl . t i t l e ( ’ Zugstab unter w e c h s e l s e i t i g e r Belastung − i s o . Verf . ’ )133 pl . g r i d (True )134 p lSt ra in , p l S t r e s s = [ ] , [ ]135 while time <= tMax :136 i f time <.5: t e n s i l e . t = time137 e l i f time <1.5: t e n s i l e . t=(1.− time )138 e l i f time <2.0: t e n s i l e . t=(time −2.)139 else : t e n s i l e . t=(time −2.)140 problem = i t e r a t e (Gain , Form , boundaryConditions , subdomains )141 s o l v e r = NewtonSolver ( ’ lu ’ )142 s o l v e r . parameters [ ’ c onve rgence_cr i t e r i on ’ ] = ’ incrementa l ’143 s o l v e r . parameters [ ’ r e l a t i v e_ t o l e r an c e ’ ] = 1 .0 e−3144 s o l v e r . s o l v e ( problem , u . vec to r ( ) )145 cauchyProj=p ro j e c t ( cauchy , Tensor )146 cauchy0 . a s s i gn ( cauchyProj )147 p lS t r a i n . append ( p r o j e c t ( eps , Tensor ) ( (xMax/2 . , yMax/2 . ) ) [ 0 ] )148 p l S t r e s s . append ( cauchyProj ( (xMax/2 . ,yMax/2 . ) ) [ 0 ] )149 pl . p l o t ( p lS t ra in , p l S t r e s s , l i n e s t y l e=’− ’ , c o l o r=’ red ’ , l i n ew id th=1, marker=’ ’ ,

markers i ze=11)150 pl . draw ( )151 pl . s a v e f i g ( ’ . / webs i t e r e . eps ’ )152 d i r e c tP r o j = p ro j e c t ( sdev0 [ k , l ]∗ deps [ k , l ] , Coe f f )153 d i r e c t a r r = d i r e c tP r o j . vec to r ( ) . array ( )154 d i r e c t b o o l = np . s i gn ( . 5 ∗ ( d i r e c t a r r +(( d i r e c t a r r ) ∗∗2 . ) ∗∗ . 5 ) )155 epsaccPro j=p ro j e c t ( epsacc+a l f a ∗ ( 2 . / 3 . ∗ depspl [ i , j ]∗ depspl [ i , j ] ) ∗∗ . 5 , Coe f f )156 Y. vec to r ( ) [ : ]= 2 5 0 . e+6+h∗ epsaccPro j . vec to r ( ) . array ( ) [ : ]157 SVM0proj=p ro j e c t (SVM0, Coe f f )158 y i e l dboo l = SVM0proj . vec to r ( ) . array ( ) >= Y. vec to r ( ) . array ( )159 a l f a . vec to r ( ) [ : ] = np . array ( y i e l dboo l ∗ d i r e c tboo l , dtype=in t )160 eps0 . a s s i gn ( p r o j e c t ( eps , Tensor ) )161 epsacc . a s s i gn ( epsaccPro j )162 Yproj=p ro j e c t (Y, Coe f f )163 f i l eU << (u , time )164 time+=dt

A.4 Zugprobe – 2D, unter wechselseitiger Belastung – isotherm, kinemati-sche Verfestigung



9 " " "−−−−−−−−−−−−newton s o l v e r i n t e r a c t i o n−−−−−−−−−−−−−−" " "10 class i t e r a t e ( NonlinearProblem ) :11 def __init__( s e l f , a , L , bc , exter_B ) :

61

A Quelltexte

12 NonlinearProblem . __init__( s e l f )13 s e l f . L = L14 s e l f . a = a15 s e l f . bc = bc16 s e l f . ex t e r = exter_B17 def F( s e l f , b , x ) :18 assemble ( s e l f . L , t en so r=b , exter ior_facet_domains=s e l f . ex t e r )19 for cond i t i on in s e l f . bc : cond i t i on . apply (b , x )20 def J ( s e l f , A, x ) :21 assemble ( s e l f . a , t en so r=A, exter ior_facet_domains=s e l f . ex t e r )22 for cond i t i on in s e l f . bc : cond i t i on . apply (A)23

24 " " "−−−−−−−−−−−mesh−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "25 xMin = .026 xMax = 1 .27 xElements = 3028 i , j , k , l , r , m, o , p , s , n , t , q= i nd i c e s (12)29 yMin = .030 yMax = .231 yElements = 632


35 " " "−−−−−−−−−f unc t i on spaces−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "36 Coef f = FunctionSpace (mesh , ’CG’ , 1)37 Vector = VectorFunctionSpace (mesh , ’CG’ , 1)38 Tensor = TensorFunctionSpace (mesh , ’CG’ , 1)39

40 " " "−−−−−−−−−−func t ions−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "41 delu = TestFunction ( Vector )42 Du = Tria lFunct ion ( Vector )43

44 u = Function ( Vector )45 u0 = Function ( Vector )46 u00 = Function ( Vector )47

48 " " "−−−−−−−−−−−−−−−−−−−−−−−−− i n i t i a l cond−−−−−−−−−−−−−−−−−−−−−−−−−" " "49 class I n i t i a lC ond i t i o n s ( Express ion ) :50 def eva l ( s e l f , vec , x ) :51 vec [ 0 ] = . 0 #u[0]=052 vec [ 1 ] = . 0 #u[1]=053 def value_shape ( s e l f ) :54 return ( 2 , )55

56 in iObj = I n i t i a lC ond i t i o n s ( )57 u00 . i n t e r p o l a t e ( in iObj )58 u0 . a s s i gn ( u00 )59 u . a s s i gn ( u0 )60


65 l e f t = compile_subdomains ( ’ x [ 0 ] <=Min+0.1∗Max ’ )66 l e f t .Min = xMin67 l e f t .Max = xMax68

69 f r on t = compile_subdomains ( ’ x [ 1 ] ==Max ’ )

62

A Quelltexte

70 f r on t .Max = yMax71


75 i n i tD i s p = Express ion ( ( ’ 0 . 0 ’ , ’ 0 . 0 ’ ) )76 bc1 = Dir ichletBC ( Vector , i n i tD i sp , l e f t )77 t e n s i l e = Express ion ( ( ’A∗ t ’ , ’ 0 ’ ) ,A=xMax/400 . , t =0.)78 bc2 = Dir ichletBC ( Vector , t e n s i l e , r i g h t )79 bc3 = Dir ichletBC ( Vector . sub (1 ) , ’ 0 ’ , r i g h t )80



85 #−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−86 de l t a = Id en t i t y (2 )87

88 " " "−−−−−−−−−−−−−time−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "89 time =.090 dt = .00591 tMax = 2 .592

93 " " "−−−−−−−−−−mater ia l:−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "94 nu , E = 0 . 3 , 200 . e+0995 E_T = .2∗E96 h=E_T/(1.−E_T/E)97 lamda = E∗nu/(1.+nu) /(1.−2.∗nu)98 mu = E/(2.+ 2 .∗nu)99


104 " " "−−−−−−−−−Strain−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "105 eps0=Function ( Tensor )106 eps = as_tensor ( . 5 ∗ ( u [ i ] . dx ( j )+u [ j ] . dx ( i ) ) , ( i , j ) )107 deps=as_tensor ( eps [ k , l ]− eps0 [ k , l ] , ( k , l ) )108

109 " " "−−−−−−−−−Stress−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "110 a l f a=Function ( Coe f f )111 cauchy0 = Function ( Tensor )112 sdev0 = as_tensor ( cauchy0 [ i , j ] −1 ./3 .∗ cauchy0 [ k , k ]∗ de l t a [ i , j ] , [ i , j ] )113 Alph0 = Function ( Tensor )114 depspl0 = Function ( Tensor )115 epsp l0 = Function ( Tensor )116 dcauchy = as_tensor (E/(1.+nu) ∗( deps [ i , j ]+nu/(1.−2.∗nu) ∗deps [ k , k ]∗ de l t a [ i , j ]− a l f a ∗ ( 3 .∗E)

/ (3 .∗E+2.∗(1.+nu) ∗h) ∗ 3 . / 2 . ∗ ( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗deps [ k , l ] /Y0∗( sdev0 [ i , j ]−Alph0 [ i , j] ) /Y0) , ( i , j ) )

117 depspl = as_tensor ( a l f a ∗1 ./h ∗ 3 . / 2 . ∗ ( sdev0 [ i , j ]−Alph0 [ i , j ] ) ∗dcauchy [ i , j ] /Y0∗ 3 . / 2 . ∗ ( sdev0 [ i, j ]−Alph0 [ i , j ] ) /Y0 , ( i , j ) )

118 cauchy = as_tensor ( cauchy0 [ i , j ]+dcauchy [ i , j ] , ( i , j ) )119 dAlph = as_tensor (h∗( depspl [ i , j ]∗ depspl [ i , j ] ∗ 2 . / 3 . ) ∗∗ . 5/Y0∗( cauchy0 [ i , j ]−Alph0 [ i , j ] ) , ( i ,

j ) )120 Alph = as_tensor (Alph0 [ k , l ]+dAlph [ k , l ] , (k , l ) )121 epsp l = as_tensor ( epsp l0 [ i , j ]+ depspl [ i , j ] , ( i , j ) )122 #von−mises e q u i v a l e n t s t r e s s i n c l u d i n g a c o r r e c t i o n f o r 3D −> 2D

63

A Quelltexte

123 f = ( 1 . / 2 . ∗ ( ( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗( sdev0 [ k , l ]−Alph0 [ k , l ] ) +(cauchy0 [ i , i ] / 3 . ) ∗∗2 . ) )−1./3.∗Y∗∗2 .

124

125 " " "−−−−−−−−−−−−−−−−−−−" " "126 Form = cauchy [ i , j ]∗ delu [ j ] . dx ( i ) ∗dx127 Gain = de r i v a t i v e (Form , u , Du)128

129 f i l eU = F i l e ( ’ . /U. pvd ’ )130

131 pl . ion ( )132 pl . c l f ( )133 pl . c l a ( )134 pl . x l ab e l ( ’ eps_11 ’ )135 pl . y l ab e l ( ’ sigma_11 ’ )136 pl . t i t l e ( ’ Zugstab unter w e c h s e l s e i t i g e r Belastung − kin . Verf . ’ )137 pl . g r i d (True )138 p lSt ra in , p l S t r e s s = [ ] , [ ]139 while time <= tMax :140 i f time <.5: t e n s i l e . t = time141 e l i f time <1.5: t e n s i l e . t = (1.− time )142 e l i f time <2.0: t e n s i l e . t = ( time −2.)143 else : t e n s i l e . t = ( time −2.)144 problem = i t e r a t e (Gain , Form , boundaryConditions , subdomains )145 s o l v e r = NewtonSolver ( ’ lu ’ )146 s o l v e r . parameters [ ’ c onve rgence_cr i t e r i on ’ ] = ’ incrementa l ’147 s o l v e r . parameters [ ’ r e l a t i v e_ t o l e r an c e ’ ] = 1 .0 e−3148 s o l v e r . s o l v e ( problem , u . vec to r ( ) )149 cauchyProj=p ro j e c t ( cauchy , Tensor )150 cauchy0 . a s s i gn ( cauchyProj )151 p lS t r a i n . append ( p r o j e c t ( eps , Tensor ) ( (xMax/2 . , yMax/2 . ) ) [ 0 ] )152 p l S t r e s s . append ( cauchyProj ( (xMax/2 . ,yMax/2 . ) ) [ 0 ] )153 pl . f i g u r e (1 )154 pl . p l o t ( p lS t ra in , p l S t r e s s , l i n e s t y l e=’− ’ , c o l o r=’ red ’ , l i n ew id th=1, marker=’ ’ ,

markers i ze=11)155 pl . draw ( )156 pl . f i g u r e (2 )157 pl . p l o t ( tar r , f a r r , l i n e s t y l e=’− ’ , c o l o r=’ red ’ , l i n ew id th=1, marker=’ ’ , markers i ze=11)158 pl . draw ( )159 a lphpro j = p ro j e c t (Alph , Tensor )160 Alph0 . a s s i gn ( a lphpro j )161 f p r o j=p ro j e c t ( f , Coe f f )162 d i r e c tP r o j = p ro j e c t ( ( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗deps [ k , l ] , Coe f f )163 d i r e c t a r r = d i r e c tP r o j . vec to r ( ) . array ( )164 d i r e c t b o o l = . 5∗ ( np . s i gn ( d i r e c t a r r ) +1.)165 y i e l dboo l = f p r o j . vec to r ( ) . array ( ) >= 0 .166 a l f a . vec to r ( ) [ : ] = np . array ( y i e l dboo l ∗ d i r e c tboo l , dtype=in t )167 eps0 . a s s i gn ( p r o j e c t ( eps , Tensor ) )168 epsp l0 . a s s i gn ( p r o j e c t ( epspl , Tensor ) )169 depspl0 . a s s i gn ( p r o j e c t ( depspl , Tensor ) )170 f i l eU << (u , time )171 time+=dt

A.5 Zugprobe – 2D, unter wechselseitiger Belastung mit Temperaturentwick-lung

64

A Quelltexte




24 " " "−−−−−−−−−−−mesh−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "25 xMin = .026 xMax = 1 .27 xElements = 1828 i , j , k , l , r , m, o , p , s , n , t , q= i nd i c e s (12)29 yMin = .030 yMax = .231 yElements = 632


35 " " "−−−−−−−−−f unc t i on spaces−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "36 Coef f = FunctionSpace (mesh , ’CG’ , 1)37 Vector = VectorFunctionSpace (mesh , ’CG’ , 1)38 Tensor = TensorFunctionSpace (mesh , ’CG’ , 1)39 Mixed = MixedFunctionSpace ( [ Vector , Coe f f ] )40

41

42 " " "−−−−−−−−−−func t ions−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "43 vMF = TestFunction (Mixed )44 dMF = Tria lFunct ion (Mixed )45 ( delu , delT ) = s p l i t (vMF)46

47 mixREF = Function (Mixed )48 mix = Function (Mixed )49 mix0 = Function (Mixed )50

51 (u , T) = s p l i t (mix )52 (u0 , T0) = s p l i t (mix0 )53

54

55 T_ini=293.1556 T_ref=T_ini57 " " "−−−−−−−−−−−−−−−−−−−−−−−−− i n i t i a l cond−−−−−−−−−−−−−−−−−−−−−−−−−" " "58 class I n i t i a lC ond i t i o n s ( Express ion ) :

65

A Quelltexte

59 def __init__( s e l f , T_ini ) :60 s e l f . T_ini=T_ini61 def eva l ( s e l f , vec , x ) :62 vec [ 0 ] = . 063 vec [ 1 ] = . 064 vec [ 2 ] = s e l f . T_ini65 def value_shape ( s e l f ) :66 return ( 3 , )67

68 in iObj = I n i t i a lC ond i t i o n s ( T_ini )69 mixREF . i n t e r p o l a t e ( in iObj )70 mix . a s s i gn (mixREF)71 mix0 . a s s i gn (mix )72


77 l e f t = compile_subdomains ( ’ x [ 0 ] <=Min+0.1∗Max ’ )78 l e f t .Min= xMin79 l e f t .Max = xMax80



87 i n i tD i s p = Express ion ( ( ’ 0 . 0 ’ , ’ 0 . 0 ’ ) )88 bc1 = Dir ichletBC (Mixed . sub (0 ) , i n i tD i sp , l e f t )89 t e n s i l e = Express ion ( ( ’A∗ t ’ , ’ 0 ’ ) , A=xMax/400 . , t =0.)90 bc2 = Dir ichletBC (Mixed . sub (0 ) , t e n s i l e , r i g h t )91 bc3 = Dir ichletBC (Mixed . sub (0 ) . sub (1 ) , ’ 0 ’ , r i g h t )92



97 " " "−−−−−−−−−−−−−time−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "98 time =.099 dt = .01

100 tMax = 2 .5101

102 de l t a=Id en t i t y (2 )103

104 " " "−−−−−−−−−−mater ia l:−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "105 nu , E = 0 . 3 , 200 . e+09106 rho = 7800 .107 E_T = .2∗E108 h=E_T/(1.−E_T/E)109 h_th=18.110 lamda = E∗nu/(1.+nu) /(1.−2.∗nu)111 mu = E/(2.+ 2 .∗nu)112 kappa = 57 .113 c_p = 466 .114 C_p = c_p∗ rho115 alpha_T = 11 . e−6116

66

A Quelltexte

117 " " "−−−−−−−−−Strain−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "118 eps0=Function ( Tensor )119 epsp l0=Function ( Tensor )120 epsp l=Function ( Tensor )121 eps = as_tensor ( . 5 ∗ ( u [ i ] . dx ( j )+u [ j ] . dx ( i ) ) , [ i , j ] )122 deps=as_tensor ( eps [ k , l ]− eps0 [ k , l ] , ( k , l ) )123

124 " " "−−−−−−−−−Stress−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "125 Y=Function ( Coe f f )126 Y0 = Constant ( 250 . e+6)127 Y. vec to r ( ) [ : ] = 250 . e+6128

129 a l f a=Function ( Coe f f )130 cauchy0 = Function ( Tensor )131 sdev0=as_tensor ( cauchy0 [ i , j ] −1 ./3 .∗ cauchy0 [ k , k ]∗ de l t a [ i , j ] , [ i , j ] )132

133 Alph0 = Function ( Tensor )134 depspl0 = Function ( Tensor )135 epsp l0 = Function ( Tensor )136

137 dcauchy=as_tensor (E/(1.+nu) ∗( deps [ i , j ]+nu/(1.−2.∗nu) ∗deps [ k , k ]∗ de l t a [ i , j ]− a l f a ∗ ( 3 .∗E)/ (3 .∗E+2.∗(1.+nu) ∗h) ∗( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗deps [ k , l ] /Y0∗ 3 . / 2 . ∗ ( sdev0 [ i , j ]−Alph0 [ i , j] ) /Y0)−( d e l t a [ k , k ]∗ lamda+2.∗mu) ∗alpha_T∗(T−T0) ∗ de l t a [ i , j ] , ( i , j ) )

138

139 depspl=as_tensor ( a l f a ∗1 ./h ∗ 3 . / 2 . ∗ ( sdev0 [ i , j ]−Alph0 [ i , j ] ) ∗dcauchy [ i , j ] /Y0∗ 3 . / 2 . ∗ ( sdev0 [ i , j]−Alph0 [ i , j ] ) /Y0 , ( i , j ) )

140

141 cauchy=as_tensor ( cauchy0 [ i , j ]+dcauchy [ i , j ] , ( i , j ) )142

143 dAlph = as_tensor ( ( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗dcauchy [ k , l ] / ( 2 . / 3 . ∗Y∗∗2 . ) ∗( cauchy0 [ i , j ]−Alph0 [i , j ] ) , ( i , j ) )

144 Alph = as_tensor (Alph0 [ k , l ]+ a l f a ∗dAlph [ k , l ] , (k , l ) )145

146 epsp l = as_tensor ( epsp l0 [ i , j ]+ depspl [ i , j ] , ( i , j ) )147 #von−mises e q u i v a l e n t s t r e s s i n c l u d i n g a c o r r e c t i o n f o r 3D −> 2D148 f =(1 . /2 .∗ ( ( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗( sdev0 [ k , l ]−Alph0 [ k , l ] ) +(cauchy0 [ i , i ] / 3 . ) ∗∗2 . ) ) −1./3.∗Y

∗∗2 .149

150 #−−−−−−−−−−−−−−−−−−151 Form = − dt∗kappa∗T. dx ( i ) ∗( delT/T) . dx ( i ) ∗dx \152 − alpha_T∗delT ∗( lamda∗ de l t a [ r , r ]+2.∗mu) ∗( deps [ k , k]− de l t a [ i , i ]∗ alpha_T∗(T−T0) ) ∗dx \153 − C_p∗delT/T∗(T−T0) ∗dx \154 + delT/T∗2 .∗mu∗( eps [ k , l ]− epsp l0 [ k , l ] ) ∗ depspl [ k , l ]∗ dx \155 + dt∗h_th∗delT/T∗ds \156 + cauchy [ i , r ]∗ delu [ i ] . dx ( r ) ∗dx157

158 Gain = de r i v a t i v e (Form , mix , dMF)159

160 f i l eU = F i l e ( ’ . /U. pvd ’ )161 f i l eT = F i l e ( ’ . /T. pvd ’ )162

163 pl . ion ( )164 pl . f i g u r e (1 , f i g s i z e = [ 7 . 6 , 1 1 . 4 ] )165 pl . c l f ( )166 pl . c l a ( )167 p lSt ra in , p l S t r e s s = [ ] , [ ]168 plT , p l t = [ ] , [ ]169 while time <= tMax :

67

A Quelltexte

170 i f time <.5: t e n s i l e . t = time171 e l i f time <1.5: t e n s i l e . t=(1.− time )172 e l i f time <2.0: t e n s i l e . t=(time −2.)173 else : t e n s i l e . t=(time −2.)174 problem = i t e r a t e (Gain , Form , boundaryConditions , subdomains )175 s o l v e r = NewtonSolver ( ’ lu ’ )176 s o l v e r . parameters [ ’ c onve rgence_cr i t e r i on ’ ] = ’ incrementa l ’177 s o l v e r . parameters [ ’ r e l a t i v e_ t o l e r an c e ’ ] = 1 .0 e−3178 s o l v e r . s o l v e ( problem , mix . vec to r ( ) )179 cauchyProj=p ro j e c t ( cauchy , Tensor )180 cauchy0 . a s s i gn ( cauchyProj )181 p lS t r a i n . append ( p r o j e c t ( eps , Tensor ) ( (xMax/2 . , yMax/2 . ) ) [ 0 ] )182 p l S t r e s s . append ( cauchyProj ( (xMax/2 . ,yMax/2 . ) ) [ 0 ] )183 pl . subplot (211)184 pl . p l o t ( p lS t ra in , p l S t r e s s , l i n e s t y l e=’− ’ , c o l o r=’ red ’ , l i n ew id th=1, marker=’ ’ ,

markers i ze=11)185 pl . subplot (212)186 pl . p l o t ( p l t , plT , l i n e s t y l e=’− ’ , c o l o r=’ red ’ , l i n ew id th=1, marker=’ ’ , markers i ze=11)187 pl . draw ( )188 pl . s a v e f i g ( ’ ’ . j o i n ( [ d i r e c to ry , ’ p l o t . eps ’ ] ) , format=’ eps ’ )189 p l t . append ( time )190 mix0 . a s s i gn (mix )191 plT . append ( p r o j e c t (T, Coe f f ) ( (xMax/2 . ,yMax/2 . ) ) )192 a lphpro j = p ro j e c t (Alph , Tensor )193 Alph0 . a s s i gn ( a lphpro j )194 f p r o j=p ro j e c t ( f , Coe f f )195 d i r e c tP r o j = p ro j e c t ( ( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗deps [ k , l ] , Coe f f )196 d i r e c t a r r = d i r e c tP r o j . vec to r ( ) . array ( )197 d i r e c t b o o l = . 5∗ ( np . s i gn ( d i r e c t a r r ) +1.)198 y i e l dboo l = f p r o j . vec to r ( ) . array ( ) >= 0 .199 a l f a . vec to r ( ) [ : ] = np . array ( y i e l dboo l ∗ d i r e c tboo l , dtype=in t )200 eps0 . a s s i gn ( p r o j e c t ( eps , Tensor ) )201 epsp l0 . a s s i gn ( p r o j e c t ( epspl , Tensor ) )202 depspl0 . a s s i gn ( p r o j e c t ( depspl , Tensor ) )203 f i l eU << (mix . s p l i t ( ) [ 0 ] , time )204 f i l eT << (mix . s p l i t ( ) [ 1 ] , time )205 time+=dt206 print time

A.6 Scherklotz – 2D, unter wechselseitiger Belastung mit Temperaturent-wicklung



9 " " "−−−−−−−−−−−−newton s o l v e r i n t e r a c t i o n−−−−−−−−−−−−−−" " "10 class i t e r a t e ( NonlinearProblem ) :11 def __init__( s e l f , a , L , bc , exter_B ) :12 NonlinearProblem . __init__( s e l f )13 s e l f . L = L

68

A Quelltexte

14 s e l f . a = a15 s e l f . bc = bc16 s e l f . ex t e r = exter_B17 def F( s e l f , b , x ) :18 assemble ( s e l f . L , t en so r=b , exter ior_facet_domains=s e l f . ex t e r )19 for cond i t i on in s e l f . bc : cond i t i on . apply (b , x )20 def J ( s e l f , A, x ) :21 assemble ( s e l f . a , t en so r=A, exter ior_facet_domains=s e l f . ex t e r )22 for cond i t i on in s e l f . bc : cond i t i on . apply (A)23

24 " " "−−−−−−−−−−−mesh−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "25 xMin = .026 xMax = .627 xElements = 1028 i , j , k , l , r , m, o , p , s , n , t , q= i nd i c e s (12)29 yMin = .030 yMax = .631 yElements = 1032


35 " " "−−−−−−−−−f unc t i on spaces−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "36 Coef f = FunctionSpace (mesh , ’CG’ , 1)37 Vector = VectorFunctionSpace (mesh , ’CG’ , 1)38 Tensor = TensorFunctionSpace (mesh , ’CG’ , 1)39 Mixed = MixedFunctionSpace ( [ Vector , Coe f f ] )40

41

42 " " "−−−−−−−−−−func t ions−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "43 vMF = TestFunction (Mixed )44 dMF = Tria lFunct ion (Mixed )45 ( delu , delT ) = s p l i t (vMF)46

47 mixREF = Function (Mixed )48 mix = Function (Mixed )49 mix0 = Function (Mixed )50

51 (u , T) = s p l i t (mix )52 (u0 , T0) = s p l i t (mix0 )53

54

55 T_ini=293.1556 T_ref=T_ini57 " " "−−−−−−−−−−−−−−−−−−−−−−−−− i n i t i a l cond−−−−−−−−−−−−−−−−−−−−−−−−−" " "58 class I n i t i a lC ond i t i o n s ( Express ion ) :59 def __init__( s e l f , T_ini ) :60 s e l f . T_ini=T_ini61 def eva l ( s e l f , vec , x ) :62 vec [ 0 ] = . 063 vec [ 1 ] = . 064 vec [ 2 ] = s e l f . T_ini65 def value_shape ( s e l f ) :66 return ( 3 , )67

68 in iObj = I n i t i a lC ond i t i o n s ( T_ini )69 mixREF . i n t e r p o l a t e ( in iObj )70 mix . a s s i gn (mixREF)71 mix0 . a s s i gn (mix )

69

A Quelltexte

72

73 " " "−−−−−−−−−−−−−−−−−−−−−−−−−boundary cond−−−−−−−−−−−−−−" " "74 top = compile_subdomains ( ’ x [ 1 ] ==Max ’ )75 top .Max = yMax76

77 bot = compile_subdomains ( ’ x [ 1 ] ==Min ’ )78 bot .Min= yMin79

80 i n i tD i s p = Express ion ( ( ’ 0 . 0 ’ , ’ 0 . 0 ’ ) ) #Wie gss81 bc1 = Dir ichletBC (Mixed . sub (0 ) , i n i tD i sp , bot )82 shear = Express ion ( ( ’A∗ t ’ , ’ 0 ’ ) , A=xMax/200 . , t =0.)83 bc2 = Dir ichletBC (Mixed . sub (0 ) , shear , top )84

85 boundaryConditions = [ bc1 , bc2 ]86


89 " " "−−−−−−−−−−−−−time−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "90 time =.091 dt = .0192 tMax = 2 .593

94 de l t a=Id en t i t y (2 )95

96 " " "−−−−−−−−−−mater ia l:−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "97 nu , E = 0 . 3 , 200 . e+0998 rho = 7800 .99 E_T = .2∗E

100 h=E_T/(1.−E_T/E)101 h_th=18.102 lamda = E∗nu/(1.+nu) /(1.−2.∗nu)103 mu = E/(2.+ 2 .∗nu)104 T_ini=293.15105 T_ref=T_ini106 kappa = 57 .107 c_p = 466 .108 C_p = c_p∗ rho109 alpha_T = 11 . e−6110

111 " " "−−−−−−−−−Strain−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "112 eps0=Function ( Tensor )113 epsp l0=Function ( Tensor )114 epsp l=Function ( Tensor )115 eps = as_tensor ( . 5 ∗ ( u [ i ] . dx ( j )+u [ j ] . dx ( i ) ) , [ i , j ] )116 deps=as_tensor ( eps [ k , l ]− eps0 [ k , l ] , ( k , l ) )117

118 " " "−−−−−−−−−Stress−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−" " "119 Y=Function ( Coe f f )120 Y0 = Constant ( 250 . e+6)121 Y. vec to r ( ) [ : ] = 250 . e+6122

123 a l f a=Function ( Coe f f )124 cauchy0 = Function ( Tensor )125

126 sdev0=as_tensor ( cauchy0 [ i , j ] −1 ./3 .∗ cauchy0 [ k , k ]∗ de l t a [ i , j ] , [ i , j ] )127

128 Alph0 = Function ( Tensor )129 depspl0 = Function ( Tensor )

70

A Quelltexte

130 epsp l0 = Function ( Tensor )131

132 dcauchy=as_tensor (E/(1.+nu) ∗( deps [ i , j ]+nu/(1.−2.∗nu) ∗deps [ k , k ]∗ de l t a [ i , j ]− a l f a ∗ ( 3 .∗E)/ (3 .∗E+2.∗(1.+nu) ∗h) ∗( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗deps [ k , l ] /Y0∗ 3 . / 2 . ∗ ( sdev0 [ i , j ]−Alph0 [ i , j] ) /Y0)−( d e l t a [ k , k ]∗ lamda+2.∗mu) ∗alpha_T∗(T−T0) ∗ de l t a [ i , j ] , ( i , j ) )

133

134 depspl=as_tensor ( a l f a ∗1 ./h ∗ 3 . / 2 . ∗ ( sdev0 [ i , j ]−Alph0 [ i , j ] ) ∗dcauchy [ i , j ] /Y0∗ 3 . / 2 . ∗ ( sdev0 [ i , j]−Alph0 [ i , j ] ) /Y0 , ( i , j ) )

135

136

137 cauchy=as_tensor ( cauchy0 [ i , j ]+dcauchy [ i , j ] , ( i , j ) )138

139 dAlph = as_tensor ( ( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗dcauchy [ k , l ] / ( 2 . / 3 . ∗Y∗∗2 . ) ∗( cauchy0 [ i , j ]−Alph0 [i , j ] ) , ( i , j ) )

140 Alph = as_tensor (Alph0 [ k , l ]+ a l f a ∗dAlph [ k , l ] , (k , l ) )141

142 epsp l = as_tensor ( epsp l0 [ i , j ]+ depspl [ i , j ] , ( i , j ) )143 #von−mises e q u i v a l e n t s t r e s s i n c l u d i n g a c o r r e c t i o n f o r 3D −> 2D144 f =(1 . /2 .∗ ( ( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗( sdev0 [ k , l ]−Alph0 [ k , l ] ) +(cauchy0 [ i , i ] / 3 . ) ∗∗2 . ) ) −1./3.∗Y

∗∗2 .145

146 #−−−−−−−−−−−−−−147 Form = − dt∗kappa∗T. dx ( i ) ∗( delT/T) . dx ( i ) ∗dx \148 − alpha_T∗delT ∗( lamda ∗3.+2.∗mu) ∗( deps [ k , k ]−3.∗alpha_T∗(T−T0) ) ∗dx \149 − C_p∗delT/T∗(T−T0) ∗dx \150 + delT/T∗2 .∗mu∗( eps [ k , l ]− epsp l0 [ k , l ] ) ∗ depspl [ k , l ]∗ dx \151 + dt∗h_th∗delT/T∗ds \152 − cauchy [ i , r ]∗ delu [ i ] . dx ( r ) ∗dx153

154 Gain = de r i v a t i v e (Form , mix , dMF)155

156 f i l eU = F i l e ( ’ ’ . j o i n ( [ d i r e c to ry , ’Data/U. pvd ’ ] ) )157 f i l eT = F i l e ( ’ ’ . j o i n ( [ d i r e c to ry , ’Data/T. pvd ’ ] ) )158 f i l eUT = F i l e ( ’ ’ . j o i n ( [ d i r e c to ry , ’Data/UT. pvd ’ ] ) )159 pl . ion ( )160 pl . f i g u r e (1 , f i g s i z e = [ 7 . 6 , 1 1 . 4 ] )161 pl . c l f ( )162 pl . c l a ( )163 p lSt ra in , p l S t r e s s = [ ] , [ ]164 plT , p l t = [ ] , [ ]165 while time <= tMax :166 i f time <.5: shear . t = time167 e l i f time <1.5: shear . t=(1.− time )168 e l i f time <2.0: shear . t=(time −2.)169 else : shear . t=(time −2.)170 problem = i t e r a t e (Gain , Form , boundaryConditions , subdomains )171 s o l v e r = NewtonSolver ( ’ lu ’ )172 s o l v e r . parameters [ ’ c onve rgence_cr i t e r i on ’ ] = ’ incrementa l ’173 s o l v e r . parameters [ ’ r e l a t i v e_ t o l e r an c e ’ ] = 1 .0 e−3174 s o l v e r . s o l v e ( problem , mix . vec to r ( ) )175 cauchyProj=p ro j e c t ( cauchy , Tensor )176 cauchy0 . a s s i gn ( cauchyProj )177 epspro j=p ro j e c t ( eps , Tensor )178 p lS t r a i n . append ( epspro j ( (xMax/2 . , 3 . ∗yMax/4 . ) ) [ 1 ] )179 p l S t r e s s . append ( cauchyProj ( (xMax/2 . , 3 . ∗yMax/4 . ) ) [ 1 ] )180 pl . subplot (211)181 pl . p l o t ( p lS t ra in , p l S t r e s s , l i n e s t y l e=’− ’ , c o l o r=’ red ’ , l i n ew id th=1, marker=’ ’ ,

markers i ze=11)

71

A Quelltexte

182 pl . subplot (212)183 pl . p l o t ( p l t , plT , l i n e s t y l e=’− ’ , c o l o r=’ red ’ , l i n ew id th=1, marker=’ ’ , markers i ze=11)184 pl . draw ( )185 p l t . append ( time )186 mix0 . a s s i gn (mix )187 plT . append ( p r o j e c t (T, Coe f f ) ( (xMax/2 . , 3 . ∗yMax/4 . ) ) )188 a lphpro j = p ro j e c t (Alph , Tensor )189 Alph0 . a s s i gn ( a lphpro j )190 f p r o j=p ro j e c t ( f , Coe f f )191 d i r e c tP r o j = p ro j e c t ( ( sdev0 [ k , l ]−Alph0 [ k , l ] ) ∗deps [ k , l ] , Coe f f )192 d i r e c t a r r = d i r e c tP r o j . vec to r ( ) . array ( )193 d i r e c t b o o l = . 5∗ ( np . s i gn ( d i r e c t a r r ) +1.)194 y i e l dboo l = f p r o j . vec to r ( ) . array ( ) >= 0 .195 a l f a . vec to r ( ) [ : ] = np . array ( y i e l dboo l ∗ d i r e c tboo l , dtype=in t )196 eps0 . a s s i gn ( epspro j )197 epsp l0 . a s s i gn ( p r o j e c t ( epspl , Tensor ) )198 depspl0 . a s s i gn ( p r o j e c t ( depspl , Tensor ) )199 f i l eU << (mix . s p l i t ( ) [ 0 ] , time )200 f i l eT << (mix . s p l i t ( ) [ 1 ] , time )201 time+=dt

72

Literaturverzeichnis


Abali 2011Abali, B.E.: Technical University of Berlin, Institute of Mechanics, Chair of ContinuumsMechanics and Material Theory, Computational Reality. http://www.lkm.tu-berlin.de/

ComputationalReality/, 2011–

Bauschinger 1886Bauschinger, J.: Über die Veränderung der Elastizitätsgrenze und der Festigkeit desEisens und Stahls durch Strecken und Quetschen, durch Erwärmen und Abkühlen unddurch oftmals wiederholte Beanspruchung. In: Mitteilungen des mechanisch-technischenLaboratoriums der Königlich Technischen Hochschule München 13 (1886), S. 1

Bertram 2011Bertram, A.: Elasticity and Plasticity of Large Deformations: An Introduction. 3. Auflage.Berlin, Heidelberg, New York : Springer-Verlag, 2011

Collins u. Houlsby 1997Collins, I.F. ; Houlsby, G.T.: Application of thermomechanical principles to the model-ling of geotechnical materials. In: Proceedings of the Royal Society of London. Series A:Mathematical, Physical and Engineering Sciences 453 (1997), Nr. 1964, S. 1975–2001

Doghri 2000Doghri, I.: Mechanics of Deformable Solids: Linear, Nonlinear, Analytical and Computa-tional Aspects. Springer-Verlag, 2000 (Engineering online library)

Huang u. Khan 1995Huang, S. ; Khan, A.S.: Continuum Theory of Plasticity. New York, London, Sydney,Toronto : John Wiley & Sons, 1995 (Wiley-Interscience publication)

Lévy 1871Lévy, M.: Extrait du Mémoire sur les équations générales des mouvements intérieurs descorps solides ductiles au delà des limites où l’élasticité pourrait les ramener à leur premierétat. In: Journal de mathématiques pures et appliquées 16 (1871), S. 369–372

Lode 1926Lode, W.: Versuche über den Einfluß der mittleren Hauptspannung auf das Fließen derMetalle Eisen, Kupfer und Nickel. In: Zeitschrift für Physik 36 (1926), Nr. 11-12, S. 913–939

Lubliner 1972Lubliner, J.: On the thermodynamic foundations of non-linear solid mechanics. In: Inter-national Journal of Non-Linear Mechanics 7 (1972), Nr. 3, S. 237 – 254

Lubliner 1990Lubliner, J.: Plasticity theory. New York : Macmillan, 1990 (Maxwell Macmillan interna-tional editions in engineering)

73

http://www.lkm.tu-berlin.de/ComputationalReality/

http://www.lkm.tu-berlin.de/ComputationalReality/


Melan 1938Melan, E.: Zur Plastizität des räumlichen Kontinuums. In: Archive of Applied Mechanics9 (1938), S. 116–126

von Mises 1913Mises, R. von: Mechanik der festen Körper im plastisch deformablen Zustand. In: Göttin.Nachr. Math. Phys 1 (1913), S. 582–592

von Mises 1928Mises, R. von: Mechanik der plastischen Formänderung von Kristallen. In: ZAMM -Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematikund Mechanik 8 (1928), Nr. 3, S. 161–185

Müller 2011Müller, W.H.: Streifzüge durch die Kontinuumstheorie. Berlin, Heidelberg, New York :Springer-Verlag, 2011

Odqvist 1933Odqvist, F.K.G.: Die Verfestigung von flußeisenähnlichen Körpern. Ein Beitrag zur Plas-tizitätstheorie. In: ZAMM – Journal of Applied Mathematics and Mechanics / Zeitschriftfür Angewandte Mathematik und Mechanik 13 (1933), Nr. 5, S. 360–363

Ogden 1984Ogden, R.W.: Non-linear elastic deformations. Chichester : Ellis Horwood Limited., 1984(Ellis Horwood series in mathematics and its applications)

Prager 1955Prager, W.: The theory of plasticity: a survey of recent achievements. In: Proceedings ofthe Institution of Mechanical Engineers 169 (1955), Nr. 1, S. 41–57

Prandtl 1924Prandtl, L.: Spannungsverteilung in plastischen Körpern. In: Proceedings of the 1stInternational Congress on Applied Mechanics (1924), S. 43–54

Press u. a. 1986Press, W.H. ; Flannery, B.P. ; Teukolsky, S.A. ; Vetterling, W.T.: NumericalRecipes. 1. Auflage. Cambridge, London, New York, New Rochelle, Melbourne, Sydney :Cambridge University Press, 1986

Reuss 1930Reuss, A.: Berücksichtigung der elastischen Formänderung in der Plastizitätstheorie. In:ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathe-matik und Mechanik 10 (1930), Nr. 3, S. 266–274

de Saint-Venant 1870Saint-Venant, B. de: Memoire sur retablissement des equations differentielles des mou-

74


vement interieurs operes dans les corps solides ductiles au dela des limites oil l’elasticitepourrait les ramener a leur premier etat. In: Comptes Rendus Acad. Sci., Paris 70 (1870),S. 1323–1325

Shield u. Ziegler 1958Shield, R.T. ; Ziegler, H.: On Prager’s hardening rule. In: Zeitschrift für angewandteMathematik und Physik ZAMP 9 (1958), Nr. 3, S. 260–276

Simo u. Hughes 1998Simo, J.C. ; Hughes, T.J.R.: Computational Inelasticity. 7. Auflage. Berlin, Heidelberg,New York : Springer-Verlag, 1998 (Interdisciplinary applied mathematics: Mechanics andmaterials)

Truesdell u. Noll 1965Truesdell, C. ; Noll, W. ; Flügge, S. (Hrsg.): Handbuch der Physik III/3. Bd. Thenonlinear field theories of mechanics. Berlin, Heidelberg, New York : Springer-Verlag, 1965

Žerovnik u. a. 2010Žerovnik, A. ; Kunc, R. ; Prebil, I.: Yield-point phenomenon in constitutive models ofcyclic plasticity. In: Computational Materials Science 49 (2010), Nr. 3, S. 473–482

Wilmański 1998Wilmański, K.: Thermomechanics of Continua. Berlin, Heidelberg, New York : SpringerVerlag, 1998

Wilmański 2008Wilmański, K.: Continuum Thermodynamics, Part 1: Foundations. Singapore : WorldScientific, 2008 (Series on Advances in Mathematics for Applied Sciences)

Ziegler 1959Ziegler, H.: A modification of Prager’s hardening rule. In: Quarterly of Applied Mathe-matics 17 (1959), Nr. 1, S. 55

Ziegler u. Wehrli 1987Ziegler, H. ; Wehrli, Ch.: The Derivation of Constitutive Relations from the free Energyand the Dissipation Function. In: Advances in applied Mechanics. 25 (1987), S. 183–238

75

Klassische Plastizität und ihre thermische Erweiterung ... · 1 Einführung 1 Einführung Die...

Documents

Transcript of Klassische Plastizität und ihre thermische Erweiterung ... · 1 Einführung 1 Einführung Die...