Institut für Mechanik Technische Universität Berlin

Kontinuumsmechanik und MaterialtheorieProf. W. H. Müller, M.Sc. B. E. Abali

Simulationstools und ihre Anwendung:Simulation der Temperaturverteilung beim Laserschweißen von

dünnen Blechen

Wintersemester 2010/11

5. April 2011

Gruppe A

Teilnehmer Studiengang Matrikelnr. E-Mail1 Liang Cheng PI 213878 chengliang@hotmail.de2 Paul Lofink PI 309284 pavidlo@gmail.com3 Felix Maucher PI 205661 felix.maucher@googlemail.com

I

Inhaltsverzeichnis

1 Laserschweißen 1

1.1 Der Erste Hauptsatz der Thermodynamik . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Das FOURIER’sche Gesetz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Die innere Energie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2 Analytische Lösung 3

2.1 Allgemein . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Eindimensionale Lösung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3 FEM (Finite-Elemente-Methode) 8

3.1 Zeitliche Diskretisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83.2 Variationelle Formulierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

4 Simulation mit FEniCS 9

4.1 Was ist FEniCS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94.2 Vorbereitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94.3 Das Programm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

5 Ergebnisse 12

A Programmcode 15

Abbildungsverzeichnis

1.1 Laserschweißen - Anwendung in Automobilindustrie (N-TV, 2009) . . . . . . . . . . . . . . 11.2 Kontrollvolumen mit Gitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22.1 Kontrollvolumen Ω mit Rand ∂Ω . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Homogenisierung der Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.3 Randbedingungen bei einem allgemeinem Kontrollvolumen . . . . . . . . . . . . . . . . . . . 44.1 schematisches Rechengebiet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

4.2 Laserverteilung Laser = Pe− x2+y2

r2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115.1 Simulationsbeispiel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125.2 Temperaturverteilung bei P= 3MW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

Tabellenverzeichnis

4.1 Materialkennwerte (Woo u. Cho, 1999, S. 704) . . . . . . . . . . . . . . . . . . . . . . . . . 10

II

Nomenklatur

ρ Dichte

κ Wärmeleitfähigkeit

σi j Spannungstensor

cp,cv isobare, isochore Wärmekapazität

Ekin,Ein kinetische Energie, innere Energie

h Konvektionskoeffizient

L Laserleistung

Q Wärme

qi spez. Wärmeflussvektor

r spez. Wärmestrahlung

T Temperatur

t Zeit

u spez. innere Energie

v spez. Volumen

W mechanische Arbeit

1 Laserschweißen

1 Laserschweißen

Laserschweißen ist heutzutage eine weit verbreitete Schweißtechnik, die vor allem in der Automobil-, Luftfahrt-und Schifffahrtindustrie eingesetzt wird. Die Vorteile sind ein sehr kurze Schweißzeit, hohe Genauigkeit, hoherArbeitsabstand, der es ermöglicht auch an schwer erreichbaren Stellen zu schweißen und der hohe Grad anAutomatisierbarkeit, der eine einfache und effiziente Fließbandproduktion erlaubt. Die gängigsten Lasertypensind Nd:YAG (Neodym dotierte Yttrium Aluminium Granat Laser) und CO2-Laser, die je nach Werkstückund Gegebenheit eingesetzt werden. Beim Laserschweißen werden zwei Werkstücke nah an einander gebracht.Genau in der Fügestelle wird die Optik des Lasers fokussiert. Damit steigt im sogenannten Brennfleck die Tem-peratur extrem schnell über die Schmelztemperatur, sodass die zu fügenden Werkstücke lokal auf schmelzenund beim Wiederabkühlen und Erstarren aus der Schmelze eine Verbindung entsteht.1

Abbildung 1.1: Laserschweißen - Anwendung in Automobilindustrie (N-TV, 2009)

In der vorliegenden Arbeit wird eine numerische Herangehensweise für die Simulation mit dem Open SourceTool FEniCS vorgestellt. Es wird die Temperaturverteilung bei einem Schweißprozess mit Hilfe der Finiten-Elemente-Methode (FEM) modelliert und berechnet. Dabei gehen wir von dünnen Blechen aus, sodass dieStrahlung des Lasers in guter Näherung instantan in einer finiten Umgebung ins Körperinnere eingebrachtwird. Somit kann sie als volumetrische Größe aufgefasst werden.

1.1 Der Erste Hauptsatz der Thermodynamik

Der erste Hauptsatz der Thermodynamik (Energieerhaltungssatz) besagt, dass die zeitliche Änderung der in-neren Energie und kinetische Energie des gesamten Systems der Summe der Leistung der äußeren Kräfte undder zu/abgeführten Wärme (Vernachlässigung anderen Energien wie chemische Energie etc.) entspricht.

Ein+ Ekin = W +Q (1.1)

Da die zu schweißenden Platten nicht bewegt werden sollen, ist die Änderung der kinetischen Energie Null.Als Kontrollvolumen wählen wir die gesamte Platte und formulieren die globale Energiebilanz.2

1Wikipedia (2011); Manke (2010, S. 62-72)2vgl. Müller (2006, S. 31-44)

1

1 Laserschweißen

—

Abbildung 1.2: Kontrollvolumen mit Gitter

ˆ

V

ρ udV

︸ ︷︷ ︸

=Ein

=

ˆ

V

σ ji

∂vi∂x j

dV

︸ ︷︷ ︸

=W

+

(

−˛

∂VqinidA+

ˆ

V

ρrdV

)

︸ ︷︷ ︸

=Q

(1.2)

Wobei:´

Vρ udV Änderung der inneren Energie

´

Vσ ji

∂vi∂x j

dV Leistung der Oberflächenspannungen

−¸

∂V qinidA Wärmefluss´

VρrdV Strahlungsleistung

Da das System nicht mechanisch beansprucht wird, kann die Leistung der Oberflächenkräfte vernachlässigtwerden.

ˆ

V

ρ udV =−˛

∂VqinidA+

ˆ

V

ρrdV (1.3)

1.2 Das FOURIER’sche Gesetz

Der spez. Wärmefluss kann nach Fourier beschrieben werden. Dabei ist der spez. Wärmefluss proportional zumGradienten der Temperatur.

qi =−κ∂T

∂xi(1.4)

Die dabei auftretende Proportionalitätskonstante ist die Wärmeleitfähigkeit κ .3

1.3 Die innere Energie

Die spez. innere Energie u ist eine Funktion von spez. Volumen und Temperatur.

u= u(v,T ) ⇒ du=

(∂u

∂T

)

v

dT +

(∂u

∂v

)

T

dv (1.5)

Unter der Annahme eines inkommpressiblen Stoffes oder unter Vernachlässigung einer Volumenänderung gilt.

u= u(T ) ⇒ du=

(∂u

∂T

)

v

dT (1.6)

Die partielle Ableitung der spez. inneren Energie nach der Temperatur bei konstantem Volumen wird als spez.isochore Wärmekapazität bezeichnet cv. Die spez. isobare4 Wärmekapazität cp und cv sind bei festen Körpern in

3vgl. Stephan u. a. (2007, S. 357ff.)4bei konstantem Druck

2

2 Analytische Lösung

guter Näherung gleich, sodass ab jetzt nur noch cp ≈ cv = c geschrieben wird.5 Diese wird für den betrachtetenTemperaturbereich als konstant angenommen.

du= cdT ⇒ u= c∂T

∂ t(1.7)

Nun können die gefundenen Beziehungen in die globale Energiebilanz (Gleichung 1.3) eingesetzt werden.Daraus ergibt sich eine Gleichung die nur noch die Temperatur enthält.

ˆ

V

ρc∂T

∂ tdV =

˛

∂Vκ

∂T

∂xinidA+

ˆ

V

ρrdV (1.8)

Der Laser wird im Material fokussiert, dadurch wird die Laserleistung direkt ins Innere eingebracht. Im vor-liegenden Modell wird dies durch den Strahlungsterm realisiert r = L(xi, t).

ˆ

V

ρc∂T

∂ tdV =

˛

∂Vκ

∂T

∂xinidA+

ˆ

V

ρL(xi, t)dV (1.9)

2 Analytische Lösung

Zunächst wollen wir uns die analytischen Lösung eines einfachen 1-D Problems anschauen: ein homogenerStab in dem sich die Wärme nur in eine Richtung ausbreitet (siehe Abbildung 2.1).

T

x

Ω ∂Ω

Abbildung 2.1: Kontrollvolumen Ω mit Rand ∂Ω

2.1 Allgemein

Dazu wird die globale Energiebilanz zunächst in eine lokale Bilanz umgeschrieben. Um dies zu bewerkstelligenwird das Flächenintegral mit dem Gaußsche Integralsatz in ein Volumenintegral umgewandelt.

ˆ

V

ρc∂T

∂ tdV =

ˆ

V

∂

∂xi

(

κ∂T

∂xi

)

dV +

ˆ

V

ρL(xi, t)dV (2.1)

Die Wärmeleitfähigkeit wird als konstant angenommen. Für stetige Felder kann das Volumenintegral an-schließend auf einen Punkt zusammengezogen werden. Hieraus folgt die lokale Bilanz oder auchWärmeleitungs-gleichung.

ρc∂T

∂ t−κ

∂ ²T

∂xi∂xi= ρL, für Ω (2.2)

Randbedinungen für ∂Ω (2.3)

Dies ist eine partielle Differentialgleichung (pDGL) im Gebiet Ω mit dem Rand ∂Ω. Dabei treten zwei Prob-leme auf.

1. Inhomogene pDGL (rechte Seite ρL(t))

(a) Lösung ist eine Linearkombination aus homogenen Lösungen und partikulärer Lösung (siehe Ab-schnitt 2.2)

5vgl. Stephan u. a. (2007, S. 104f.)

3

2 Analytische Lösung

2. Inhomogene Randbedingungen

(a) Lösung: Substitution bei vorgegebenen Randtemperatur T = T −αT

T

T

T

Tx

Abbildung 2.2: Homogenisierung der Randbedingungen

(b) Lösung: Substitution bei vorgegebener Randwärmefluss T = T −β ∂T∂xi

niDabei können folgenden Randbedingungstypen auftreten:

i. Dirichlet: α 6= 0, β = 0

ii. Neumann: α = 0, β 6= 0

iii. Robin: α 6= 0, β 6= 0

Als Beispiel werden ein Dirichlet- und ein Neumannrand untersucht.

ΓD

ΓN

Abbildung 2.3: Randbedingungen bei einem allgemeinem Kontrollvolumen

∂Ω = ΓD∪ΓN (2.4)

Nun ist das Problem vollständig definiert.Die Wärmeleitungsgleichung

ρc∂T

∂ t−κ

∂ ²T

∂xi∂xi= ρL (2.5)

wird nun substituiert mit:

T = T − αT∣∣ΓD

(2.6)

T = T − β∂T

∂xini

∣∣∣∣ΓN

(2.7)

Damit folgt:

ρc∂ T

∂ t

∣∣∣∣Ω

− κ∂ 2T

∂xi∂xi

∣∣∣∣Ω

= −ακ∂ 2T

∂xi∂xi

∣∣∣∣ΓD

− βκ∂ 3T

∂xi∂xi∂x jn j

∣∣∣∣ΓN

+ ρcα∂T

∂ t

∣∣∣∣ΓD

︸ ︷︷ ︸

=0

+ ρcβ∂ 2T

∂ t∂xini

∣∣∣∣ΓN

︸ ︷︷ ︸

=0

+ ρL|Ω

(2.8)Die beiden letzten Terme fallen weg, weil die Randbedingungen sich zeitlich nicht ändern und damit auch dieSubstitutionsfunktion T nicht von der Zeit abhängt. Analoges gilt für den Neumannrand.

4

2 Analytische Lösung

• T = T (x)∣∣ΓD

= const.t ⇒ ∂T∂ t

∣∣∣ΓD

= 0

• ∂T∂xi

ni =∂T∂xi

ni

∣∣∣ΓN

= const.t ⇒ ∂ 2T∂ t∂xi

ni

∣∣∣ΓN

= 0

ρc∂ T

∂ t

∣∣∣∣Ω

− κ∂ 2T

∂xi∂xi

∣∣∣∣Ω

= −ακ∂ 2T

∂xi∂xi

∣∣∣∣ΓD

︸ ︷︷ ︸

=0

− βκ∂ 3T

∂x j∂xi∂xin j

∣∣∣∣ΓN

︸ ︷︷ ︸

=0

+ ρL|Ω (2.9)

Die beiden letzten Terme fallen weg, wenn die Substitutionsfunktion linear ist, was im vorliegenden Fall völligausreicht.

ρc∂ T

∂ t−κ

∂ 2T

∂xi∂xi= ρL in Ω (2.10)

2.2 Eindimensionale Lösung

Als einfachstes Beispiel wird ein Stab der Länge l analysiert (siehe Abbildung 1.2). Für diesen gilt:

ρc∂T

∂ t−κ

∂ 2T

∂x∂x= ρL(x, t) , auf Ω = (0, l) (2.11)

Dabei werden beispielhaft folgende Randbedingungen angenommen:

T (x= 0) = T1, T (x= l) = T2 (2.12)

für 1D kann man die inhomogenen Randbedingungen mit folgender einfachen Substitution umgehen:

T − T = T ⇒ T = T + T , mit T = T1(1−x

l)+T2

x

l(2.13)

Damit reduziert sich das Problem nur noch auf T , welches nun aber mit homogenen Randbedingungen gesuchtist.Die Anfangsbedingung ist:

T (x, t = 0) = ϑ0 (x) ⇒ T = ϑ0− T = ϑ1 (x) (2.14)

Mit der Substitution aus Gleichung 2.13 folgt für die beschreibende pDGL analog zum Abschnitt 2.1:

ρc∂ T

∂ t−κ

∂ 2T

∂x∂x= ρL−ρc

∂ T

∂ t︸ ︷︷ ︸

=0

+ κ∂ 2T

∂x∂x︸ ︷︷ ︸

=0, T=lineare Fkt.

(2.15)

Die Lösung dieser linearen partiellen Differentialgleichung lässt sich additiv aufspalten in die homogene Lö-sung Th und partikulären Lösung Tp.

T = Th+ Tp (2.16)

Zunächst die Lösung der homogenen pDGL:

ρc∂ Th∂ t

−κ∂ 2Th

∂x∂x= 0 (2.17)

Dies führt auf den Bernoulli-Separationsansatz:

Th = χ (x)θ (t) (2.18)

Einsetzen in die pDGL liefert:

ρ cχ θ −κ χ ′′θ = 0 ⇔ ρc

κ

θ

θ=

χ ′′

χ= const. =−λ ∈R (2.19)

5

2 Analytische Lösung

Da die linke Seite nur von t abhängt und die rechte Seite nur von x abhängt, können beide Seiten höchstensnoch konstant sein. Dadurch kann man die partielle DGL in zwei gewöhnliche DGLen überführen.

θ +κλ

ρcθ = 0 mit dem Ansatz: θ (t) = Ae

− κλρc t (2.20)

χ ′′+λ χ = 0 mit dem Ansatz: χ (x) = Bsin(√

λx)

+C cos(√

λx)

(2.21)

Mit den nun homogenen Randbedingungen:

T (x= 0) = T (x= l) = 0 ⇒ χ (0)θ (t)!= 0

!= χ (l)θ (t) (2.22)

Da diese Bedingung für alle Zeiten t gelten sollen und es mit Sicherheit eine Zeit gibt wo θ (t) 6= 0 gilt, könnendiese Gleichungen nur durch χ erfüllt werden:

χ (x= 0) = 0, χ (x= l) = 0 (2.23)

⇒ χ (0) =C!= 0 (2.24)

χ (l) = Bsin(√

λ l)

!= 0 ⇒

√

λkl = πk ⇒ λk =π2

l2k2 (2.25)

Die allgemeine Lösung kann sich aus beliebigen Bestandteilen zusammensetzten und kann deshalb als Summerdargestellt werden.

χ (x) =∞

∑k=0

Bk sin(√

λkx)

=∞

∑k=0

Bk sin

(kπ

lx

)

(2.26)

θ (t) =∞

∑k=0

Ake− κλk

ρc t (2.27)

Damit ist die vollständige homogene Lösung:

Th (x, t) =∞

∑k=0

AkBk︸︷︷︸

Dk

sin

(kπ

lx

)

e− κλk

ρc t =∞

∑k=0

Dk sin

(kπ

lx

)

e− κλk

ρc t (2.28)

Um die Lösung an die Anfangsbedingung

T (x,0) = Th (x,0)+ Tp (x,0) = ϑ1 (x) (2.29)

anzupassen kann man ohne Beschränkung der Allgemeinheit, folgenden Ansatz wählen:

Th (x,0) = ϑ1 (x) und Tp (x,0) = 0 (2.30)

Wenn man nun die Funktion ϑ1 (x) in eine Fourierreihe entwickelt, dann kann man durch einen Koeffizienten-vergleich die Konstanten Dk bestimmen.

ϑ1 (x) =∞

∑k=0

ck sin

(kπ

lx

)

mit ck =2

l

ˆ l

0ϑ1 (x) sin

(kπ

lx

)

dx (2.31)

Und damit gilt:

Dk =2

l

ˆ l

0ϑ1 (x) sin

(kπ

lx

)

dx (2.32)

Für die partikuläre Lösung kann man ebenfalls einen Sinusansatz wählen.

Tp (x, t) =∞

∑n=0

θn (t)sin(nπ

lx)

(2.33)

Wobei für die Koeffizientenθn (t = 0) = 0 (2.34)

6

2 Analytische Lösung

gilt, um den Ansatz in Gleichung 2.30 zu erfüllen.Der Ansatz für die partikuläre Lösung wird nun in die pDGL eingesetzt.

∞

∑n=0

[ρcθn (t)+κλnθn (t)

]sin(nπ

lx)

= ρL(x, t) (2.35)

Nun wird die rechte Seite in eine Fourier-Reihe entwickelt.

L(x, t) =∞

∑n=0

Ln (t)sin(nπ

lx)

mit Ln (t) =2

l

lˆ

0

L(x, t) sin(nπ

lx)

dx (2.36)

⇒∞

∑n=0

[ρcθn (t)+κλnθn (t)−ρLn (t)

]

︸ ︷︷ ︸

fn(t)

sin(nπ

lx)

= 0 (2.37)

Um diese Gleichung zu lösen wird auf beiden Seiten mit sin(mπlx)multipliziert und dann über die gesamte

Länge integriert.lˆ

0

∞

∑n=0

fn (t)sin(nπ

lx)

sin(mπ

lx)

dx= 0 (2.38)

Mit den Orthogonalitätseigenschaften der Sinusfunktionenˆ l

0sin(nπ

lx)

sin(mπ

lx)

dx=l

2δnm für n,m = 1,2 . . . (2.39)

lässt sich die Gleichung vereinfachen.∞

∑n=1

fn (t)l

2δnm = 0 und damit fm (t) = 0 (2.40)

Dadurch entsteht ein System aus linearen DGLen:

θm (t)+κλm

ρcθm (t) =

1

cLm (t) (2.41)

Diese DGLen können gelöst werden, in dem auf beiden Seiten mit eκλmρc t multipliziert wird.

(

θm (t)+κλm

ρcθm (t)

)

eκλmρc t

︸ ︷︷ ︸

ddt

(

θm(t)eκλmρc t

)

=1

cLm (t)e

κλmρc t (2.42)

Eine Integration über t liefert dann.

θm (t)eκλmρc t = θm (0)

︸ ︷︷ ︸

=0, Gl. 2.34

+1

c

ˆ t

0Lm (s)e

κλmρc s

ds (2.43)

⇒ θm (t) =1

c

ˆ t

0Lm (s)e

κλmρc (s−t)

ds (2.44)

Mit dem oben gemachten Ansatz (Gleichung 2.30) lässt sich nun die Gesamtlösung hinschreiben.

T (x, t) = T1

(

1− x

l

)

+T2x

l︸ ︷︷ ︸

T

+∞

∑k=0

Dk sin

(πk

lx

)

e− κπ2k2

ρcl2t

︸ ︷︷ ︸

Th

+∞

∑n=0

θn (t)sin(nπ

lx)

︸ ︷︷ ︸

Tp

(2.45)

mit:

• Dk =2l

´ l

0 ϑ1 (τ)sin(kπl

τ)dτ = 2

l

´ l

0

(ϑ0 (τ)−T1(1− τ

l)−T2

τl

)sin(kπl

τ)dτ

• θn (t) =2cl

´ t

0

´ l

0 L(ξ ,s)sin(nπl

ξ)dξ e

κλnρc (s−t)

ds

Diese Lösung erfüllt die partielle Differentialgleichung und alle Rand- und Anfangsbedingungen. Man erkennthieraus wie umfangreich und aufwendig sich die analytische Lösung selbst für ein eindimensionales Problemgestaltet, daher soll nun im folgenden eine numerische Herangehensweise vorgestellt werden.

7

3 FEM (Finite-Elemente-Methode)

3 FEM (Finite-Elemente-Methode)

Damit die numerisch mit der Finiten-Elemente-Methode untersucht werden kann, muss es zunächst noch einwenig präpariert werden. Ausgangspunkt ist nun wieder die globale Energiebilanz:

ˆ

V

∂T

∂ tdV =

ˆ

V

κ

ρc︸︷︷︸

d

∂ 2T

∂xi∂xidV +

ˆ

V

L

c︸︷︷︸

f

dV (3.1)

Als natürlichste Randbedingungen werden gemischte oder Robin’sche Randbedingungen angenommen.

qini = h(T −T0) auf ∂V (3.2)

Dies soll zu allen Zeiten t gelten und stellt den Fluss an der Oberfläche dar, der aus der Temperaturdifferenzzwischen den Randelementen des Bleches ∂V und der Umgebungstemperatur T0 resultiert.6 Als Anfangsbedin-gung wird für das gesamte Blech die Umgebungstemperatur angenommen.

T (t = 0) = T0 auf V (3.3)

Damit ist das Problem vollständig definiert.

3.1 Zeitliche Diskretisierung

Die Zeitlich Ableitung wird durch einen rückwärtigen Differenzenquotienten diskretisiert:7

∂T

∂ t=

T (t)−T (t−∆t)

∆t+O (∆t) =

T k−T k−1

∆t+O (∆t) (3.4)

Dabei bezeichnet k immer den aktuellen Zeitschritt.ˆ

V

T k−T k−1

∆tdV −

ˆ

V

d∂ 2T k

∂xi∂xidV −

ˆ

V

f kdV = 0 (3.5)

3.2 Variationelle Formulierung

Bei der variationelle Formulierung multipliziert man die obige Gleichung unter dem Integral mit einer Test-funktion δT .

ˆ

V

T k−T k−1

∆tδTdV −

ˆ

V

d∂ 2T k

∂xi∂xiδTdV −

ˆ

V

f kδTdV = 0 (3.6)

Wenn diese Gleichung für alle Testfunktionen δT erfüllt ist, dann ist T k die exakte Lösung für den aktuellenZeitschritt.8 Wenn wir jedoch nur eine Näherungslösung suchen dann genügt es nur über einen endlich Menge,bzw. besser gesagt über einen endlich dimensionalen Funktionenraum von Testfunktionen zu ’varieren’. DieMenge über die variiert wird ist der Sobolev-Raum H1 (V ).

H1(V ) =

v

∣∣∣∣∣∣

ˆ

V

v2+

(n

∑i=1

∂v

∂xi

)2

dV < ∞

(3.7)

Dieser Raum wird nach der Galerkin-Methode diskretisiert, also durch einen finiten Funktionenraum Vh ⊂H1 (V ) ersetzt.9

Vh =v ∈C0 (V )

∣∣v ist ein lineares Polynom auf jedem Element

(3.8)

6vgl. Woo u. Cho (1999, S. 697)7vgl. Plato (2006, S. 228ff.)8vgl. Zulehner (2008, S. 9f.)9vgl. Sete (2010, S. 58ff.)

8

4 Simulation mit FEniCS

Dabei istC0 (V ) die Menge der stetigen Funktionen im VolumenV .10 Nun sind wir fast schon so weit das ganzezu implementieren, doch vorher muss noch der zweite Term in der obige Gleichung leicht modifiziert werden.Dazu nutzen wird die Produktregel.

∂

∂xi

(∂T k

∂xiδT

)

=∂ 2T k

∂xi∂xiδT +

∂T k

∂xi

∂δT

∂xi⇒ ∂ 2T k

∂xi∂xiδT =

∂

∂xi

(∂T k

∂xiδT

)

− ∂T k

∂xi

∂δT

∂xi(3.9)

Damit folgt für das zweite Integral in Gleichung 3.6:

ˆ

V

d∂ 2T k

∂xi∂xiδTdV =

ˆ

V

∂

∂xi

(

d∂T k

∂xiδT

)

dV −ˆ

V

d∂T k

∂xi

∂δT

∂xidV =

˛

∂Vd ni

∂T k

∂xiδTdA−

ˆ

V

d∂T k

∂xi

∂δT

∂xidV

(3.10)Für das Oberflächenintegral kann die Randbedingung eingesetzt werden.

qini =−κ∂T k

∂xini = h

(

T k−T0

)

⇒ ∂T k

∂xini =− h

κ

(

T k−T0

)

(3.11)

Damit kann nun Gleichung 3.6 wie folgt umgeschrieben werden:

ˆ

V

T k−T k−1

∆tδTdV +

˛

∂V

h

ρc

(

T k−T0

)

δTdA+

ˆ

V

d∂T k

∂xi

∂δT

∂xidV −

ˆ

V

f kδTdV = 0 (3.12)

4 Simulation mit FEniCS

4.1 Was ist FEniCS

FEniCS ist eine Sammlung von Freie Software Paketen (open source), die mit C++ oder Python programmiertwerden. Ziel ist es damit partielle Differentialgleichungen mit der Finiten-Elemente-Methode und Lineare Al-gebra Problem mit den entsprechenden Lösern automatisch zu simulieren und zu rechnen und damit dem Inge-nieur ein mächtiges und freies Instrument an die Hand zu geben.11 Entstanden ist es aus einem Forschungspro-jekt der University of Chicago und der Chalmers University of Technology.12

4.2 Vorbereitung

Um Gleichung 3.12 in FEniCS einzubinden, schreiben wir alle bekannten Größen auf die rechte Seite und alleunbekannten auf die linke. Wobei wir uns von einem Zeitschritt zum nächsten hangeln werden und deshalbT k−1 ebenfalls als bekannt voraussetzten.ˆ

V

T kδT dV +

ˆ

V

∆t d∂T k

∂xi

∂δT

∂xidV +

˛

∂V

h∆t

ρcT kδT dA=

ˆ

V

T k−1δTdV +

˛

∂V

h∆t

ρcT0δT dA+

ˆ

V

∆t f kδT dV

(4.1)

4.3 Das Programm

Nun schauen wir uns den Aufbau des Programms im Detail an. Der komplette Code befindet sich im AnhangA.Die ersten beiden Zeilen importieren die für diese Berechnungen notwendigen Klassen. Diese sind unter an-derem alle DOLFIN Klassen zum Generieren von Gittern, Definieren der entsprechenden Funktionenräumeund natürlich die entsprechenden algebraischen Löser.13

from d o l f i n impor t *from numpy impor t a r r a y

10siehe auch Zulehner (2008, S. 35ff.)11vgl. Fenics (2010, S. 123)12Wikipedia (2010)13vgl. Langtangen (2009, S. 9ff.)

9

4 Simulation mit FEniCS

Als nächstes werden Materialparameter definiert. Als Beispiel haben wir Werte für Baustahl verwendet.

Eigenschaft Dichte ρ Wärmekapazität c Wärmeleitfähigkeit κ Konvektionskoeffizient h

Wert 7860 kg

m3 624 JkgK

30 WmK

18 Wm2K

Tabelle 4.1: Materialkennwerte (Woo u. Cho, 1999, S. 704)

Diese tauchen auch im Programm auf.

rho =7860.0 # D ich t e i n kg /m^3c =624.0 # Wae rmekapaz i t a e t i n J / ( kg K )kappa =30.1 # Wae rm e l e i t f a e h i g k e i t i n W / (m K)h=18.0 # spez . K o n v e k t i o n s k o e f f i z i e n td=kappa / rho / c

Nun folgt die Dimensionierung der Platte, die in der vorliegenden Berechnung als quadratische angesetzt wurde(siehe Abbildung 1.2).

l =0 .1 # Laenge und B r e i t e de r q u a d r a t i s c h e n P l a t t e i n md i ck e =0.001 # Dicke de r P l a t t e m

Als nächstes werden die Lasereigenschaften dem Programm mitgeteilt.

P=2000000 # L a s e r l e i s t u n g i n W / kgr a d i u s =0.002 # Radius des L a s e r s i n mspeed =0.01 # g e s c hw i n d i g k e i t i n y r i c h t u n g i n m/ s

Nachdem nun die einzelnen Größen dem Programm mitgeteilt wurden wird im nächsten Schritt das Gittergeneriert.

mesh = Box ( 0 . 0 , 0 . 0 , 0 . 0 , l , l , d i cke , 1 0 , 1 0 , 5 )V = Func t i onSpace ( mesh , " Lagrange " , 1 )

Mit dieser Zeile wird ein quaderförmiges Gitter erzeugt mit den Eckpunkten A(0,0,0) und B(l, l,dicke) (sieheAbbildung 4.1). Das Gitter ist in Abbildung 1.2 zu sehen. Außerdem werden die Seiten in 10 x 10 x 5 Abschnitteeingeteilt auf denen dann die entsprechenden dreieckigen finiten Elemente definiert werden.

A

B

Abbildung 4.1: schematisches Rechengebiet

Anschließend wird auf diesem Gitter der endlichdimensionale Funktionenraum definiert der aus linearen Poly-nomen vom Grad kleiner gleich 1 sein soll, was in der Literatur auch als Lagrangeelemente vom Grad 1 oderlineare Dreieckselemente bekannt ist.14 Das entspricht gerade dem diskreten Funktionenraum Vh aus unsererAufgabenstellung. Als nächstes wird der zeitliche Rahmen definiert.

t =0 .0T=10.0d t =0 .1

Die Zeit soll also bei t0 = 0s beginnen und mit der Schrittweite dt = 0.1s stetig ansteigen bis sie die EndzeitT = 10s erreicht.

# Def ine i n i t i a l c o n d i t i o nT0=300.0u0 = Func t i o n (V)u0 = Exp r e s s i o n ( " 3 0 0 . 0 " )

14vgl. Sete (2010, S. 61f.)

10

4 Simulation mit FEniCS

Hier wurde die Anfangsbedingung implementiert, dabei wurde zunächst T0 der Wert 300.0K zugewiesen unddieser dann einer Feldfunktion u0 zugewiesen. Zuvor wurde noch u0 als solch eine Feldfunktion definiert, diegerade im Funktionenraum Vh lebt.

# Def ine v a r i a t i o n a l problemu = T r i a l F u n c t i o n (V)v = Te s t F unc t i o n (V)

Nun geht es ans Eingemachte, u stellt die gesuchte Temperaturverteilung T (x,y,z, t) und wird deshalb Trial-Function im Funktionenraum Vh definiert. v beschreibt die Testfunktion über die variiert oder besser gesagtgetestet wird, was in unserem Problem δT entspricht. Theoretisch hätte man die Namen der Größen auchentsprechend unseren Gegebenheiten anpassen können, aber so wird es in der Literatur häufig verwendet unddeshalb haben wir es auch dabei belassen. Im nächsten Schritt folgt die linke Seite von Gleichung 4.1.

l h s = u*v*dx + d* d t *Dx( u , i )*Dx( v , i )* dx + h* d t / rho / c*u*v* ds

Dabei werden die Volumenintegral durch dx repräsentiert und die Integrale über die Oberflächen bzw. Rän-

der durch ds. Der Gradient wird durch die Funktion D(·, i), was im Prinzip ∂ (·)∂xi

entspricht. Man beachte die

unverkennbare Ähnlichkeit zwischen den beiden Formulierung. Dies stellt wahrscheinlich auch den größtenVorteil von FEniCS für den Ingenieur dar, dass man die hergeleiteten Feldgleichungen fast eins zu eins inFEniCS einprogrammieren kann.

f i l e = F i l e ( " run8 / po i s s on . pvd " )f i l e << u

In diesem Schritt wird eine Datei erstellt mit der Endung .pvd. In diese Datei werden später für die einzelnenZeitschritte mit der zweiten Anweisung die gerade gefundenen Temperaturwerte für jeden Punkt x, y, z hinein-geschrieben.

wh i l e t <=T :Lase r = Exp r e s s i o n ( " P*exp (−(pow ( x [0]−0.5* l * ( 1 + 0 . 5 * . . .

s i n ( 6 . 2 8 / T* t ) ) , 2 ) + pow ( x [1]−0.01* t , 2 ) ) / 0 . 0 0 0 0 4 ) " )Lase r . r a d i u s = r a d i u sLase r . P=PLase r . T=TLase r . l = lLa se r . t = t

r h s = ( Lase r / c* d t + u0 )* v*dx + h* d t / rho / c*T0*v*ds

Nun wird eine while-Schleife gestartet die solange wiederholt wird bis die Zeit t die Endzeit T erreicht. In derSchleife wird als erstes die Laserfunktion mit der Expression-Anweisung definiert. Die Laserfunktion wird alseine Gauß’sche Verteilung der Leistung modelliert.

Abbildung 4.2: Laserverteilung Laser = Pe− x2+y2

r2

Die folgenden Zeilen geben der Expression-Anweisung zu verstehen welche Variablen mit den entsprechendenZeichen in der Anweisung verknüpft werden sollen. In der letzten Zeile wird die rechte Seite von Gleichung

11

5 Ergebnisse

4.1 implementiert. Das muss innerhalb der while-Schleife stehen, da sich der Laser bewegt und deshalb einezeitlich Abhängigkeit besitzt. Als Ergebnis des vorherigen Zeitschrittes T k−1 steht repräsentativ u0, das injedem Zeitschritt mit der aktuellen Lösung überschrieben wird.

problem = Va r i a t i o n a l P r o b l em ( lh s , r h s )u = problem . s o l v e ( )

In diesen Zeilen wird das Variationsproblem rechnerfreundlich formuliert in dem man ihm nur mitteilt was dielinke und was die rechte Seite der Gleichung ist. Den Rest erledigt FEniCS von selbst und löst in der nächstenZeile das Variationsproblem auch schon.

t += d tu0=u

Am Schluss der Schleife wird noch die Zeit inkrementiert und die aktuelle Lösung überschrieben damit diesedann im nächsten Zeitschritt als vorherige Lösung zur Verfügung steht. Am Schluss befindet sich noch einePlotanweisung.

# P l o t s o l u t i o np l o t ( u , i n t e r a c t i v e =True )

5 Ergebnisse

Im endgültigen Programm befinden sich 4 while-Schleifen, die 4 Teilprozesse darstellen. Zunächst fährt derLaser in einer S-Kurve von oben nach unten über die Platte, danach wird die Platte eine Zeit lang ruhen gelassendamit man auch die Konvektion der Wärme und damit auch der Temperatur beobachten kann. Und anschließendwird die Platte noch von links nach rechts in einer S-Kurve durchfahren. Damit soll ein Schweißprozess vonvier gekrümmte Platten simulierte werden, die sehr nah an einander gelegt werden.

Abbildung 5.1: Simulationsbeispiel

Wenn der Laser über die Fügestellen fährt schmelzen die entsprechenden Bereiche der Bleche und verbindensich beim Wiederabkühlen. Der Schmelzvorgang kann mit dem vorgestellten Modell zwar nicht simuliert wer-den, aber man kann z.B. durch Tests herausfinden welche Leistung des Lasers nötig ist um die Schmelztemper-atur vom Werkstück zu erreichen. Die Schmelztemperatur von Stahl beträgt ca. 1800K.Die Ergebnisse können mit dem Postprocessingtool ParaView dargestellt und auch modifiziert werden.

12

5 Ergebnisse

Abbildung 5.2: Temperaturverteilung bei P= 3MW

Man erkennt, dass bei einer eingestrahlten Leistung von P = 3 MWkg

die Schmelztemperatur von Stahl ungefährerreicht wird. Dabei verteilt sich diese Leistung auf ein kleines Zylinderelement mit dem Radius r = 2mm undeiner Höhe von h = 1mm, was der Dicke des Bleches entspricht. Daraus resultiert ein Volumen von Vhohl ≈13mm3 bzw. einer Masse von mhohl = ρ ·Vhohl ≈ 9,88 · 10−5 kg. Das entspricht einer Stahlungsleistung vonp= P ·mhohl ≈ 300W . Solche Leistungen sind durchaus realisierbar mit üblichen CO2-Lasern.15

15vgl. Woo u. Cho (1999, S. 706ff.)

13

Literatur

Literatur

[Fenics 2010] FENICS: FEniCS. http://www.feni sproje t.org/. Version: November 2010

[Langtangen 2009] LANGTANGEN, H. P.: A FEniCS Tutorial. Fenics, 2009

[Manke 2010] MANKE, Nico: Schweißen im Stahlbau in der historischen Entwicklung. Hochschule Neubran-denburg, 2010

[Müller 2006] MÜLLER, Wolfgang H.: Grundlagen der Kontinuumstheorie I / Tensoranalysis. TechnischeUniversität Berlin, Institut für Mechanik, 2006/2007

[N-TV 2009] http://www.n-tv.de/img/41/417827/O_1000_680_680_lasers hweissen.jpg[Plato 2006] PLATO, Robert: Numerische Mathematik kompakt: Grundlagenwissen für Studium und Praxis. 3.,aktualisierte und verbesserte Auflage. Vieweg, 2006

[Sete 2010] SETE, Olivier: Vorlesungsnotizen zur Veranstaltung ’Numerik II für Ingenieure’ von Prof. Dr. Jörg

Liesen. TU-Berlin, 2010

[Stephan u. a. 2007] STEPHAN, Peter ; SCHABER, Karlheinz ; MAYINGER, Franz: Thermodynamik Grundlagenund technische Anwendung. Bd. 1: Einstoffsysteme. Springer-Lehrbuch, 2007

[Wikipedia 2010] WIKIPEDIA: FEniCS Project. http://en.wikipedia.org/wiki/FEniCS_Proje t.Version: Dezember 2010

[Wikipedia 2011] WIKIPEDIA: Schweißen. http://de.wikipedia.org/wiki/S hwei%C3%9Fen.Version:März 2011

[Woo u. Cho 1999] WOO, H G. ; CHO, H S.: Three-dimensional temperature distribution in laser surface hard-ening processes. In: Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering

Manufacture 213 Part B (1999), S. 695–712

[Zulehner 2008] ZULEHNER, Walter: Numerische Mathematik: Eine Einführung anhand von Differentialgle-

ichungsproblemen. Bd. 1:Stationäre Probleme. Birkhäuser, 2008

14

A Programmcode

A Programmcode

from d o l f i n impor t *from numpy impor t a r r a y

rho =7860.0 # D ich t e i n kg /m^3c =624.0 # Wae rmekapaz i t a e t i n J / ( kg K )kappa =30.1 # Wae rm e l e i t f a e h i g k e i t i n W / (m K)h=18.0 # spez . K o n v e k t i o n s k o e f f i z i e n t

d=kappa / rho / c

l =0 .1 # Laenge und B r e i t e de r q u a d r a t i s c h e n P l a t t e i n md i ck e =0.001 # Dicke de r P l a t t e m

P=3000000 # L a s e r l e i s t u n g be i ca . W / mr a d i u s =0.002 # Radius des L a s e r s i n mspeed =0.01 # g e s c hw i n d i g k e i t i n y r i c h t u n g i n m/ s

# C r e a t e mesh and d e f i n e f u n c t i o n spacemesh = Box ( 0 . 0 , 0 . 0 , 0 . 0 , l , l , d i cke , 1 0 , 1 0 , 5 )V = Func t i onSpace ( mesh , " Lagrange " , 1 )

t =0 .0T=10.0d t =0 .1

# Def ine i n i t i a l c o n d i t i o nT0=300.0u0 = Func t i o n (V)u0 = Exp r e s s i o n ( " 3 0 0 . 0 " )

# Def ine v a r i a t i o n a l problemu = T r i a l F u n c t i o n (V)v = Te s t F unc t i o n (V)

l h s = u*v*dx + d* d t *Dx( u , i )*Dx( v , i )* dx + h* d t / rho / c*u*v* ds

u = Func t i o n (V)f i l e = F i l e ( " run8 / po i s s on . pvd " )

wh i l e t <=T :Lase r = Exp r e s s i o n ( " P*exp (−(pow ( x [0]−0.5* l * (1+0 . 5* . . .

s i n ( 6 . 2 8 / T* t ) ) , 2 ) + pow ( x [1]−0.01* t , 2 ) ) / 0 . 0 0 0 0 4 ) " )Lase r . r a d i u s = r a d i u sLase r . P=PLase r . T=TLase r . l = lLa se r . t = tr h s = ( Lase r / c* d t + u0 )* v*dx + h* d t / rho / c*T0*v*dsproblem = Va r i a t i o n a l P r o b l em ( lh s , r h s )p r i n t " Time " , t

15

A Programmcode

u = problem . s o l v e ( )t += d tu0=uf i l e << u

wh i l e t <=2*T :Lase r = Exp r e s s i o n ( " 0 . 0 " )r h s = ( Lase r / c* d t + u0 )* v*dx + h* d t / rho / c*T0*v*dsproblem = Va r i a t i o n a l P r o b l em ( lh s , r h s )p r i n t " Time " , tu = problem . s o l v e ( )t += d tu0=uf i l e << u

wh i l e t <=3*T :Lase r = Exp r e s s i o n ( " P*exp (−(pow ( x [1]−0.5* l * . . .

(1−0.5* s i n ( 6 . 2 8 / T*( t−2*T ) ) ) , 2 ) + pow ( x [0] −0 .01*( t−2*T ) , 2 ) ) / 0 . 0 0 0 0 4 ) " )Lase r . r a d i u s = r a d i u sLase r . P=PLase r . T=TLase r . l = lLa se r . t = tr h s = ( Lase r / c* d t + u0 )* v*dx − d* d t * a l ph a *T0*v*dsproblem = Va r i a t i o n a l P r o b l em ( lh s , r h s )p r i n t " Time " , tu = problem . s o l v e ( )t += d tu0=uf i l e << u

wh i l e t <=4*T :Lase r = Exp r e s s i o n ( " 0 . 0 " )r h s = ( Lase r / c* d t + u0 )* v*dx − d* d t * a l ph a *T0*v*dsproblem = Va r i a t i o n a l P r o b l em ( lh s , r h s )p r i n t " Time " , tu = problem . s o l v e ( )t += d tu0=uf i l e << u

# P l o t s o l u t i o np l o t ( u , i n t e r a c t i v e =True )

16