auf fungales Wachstum
Mathematische Modellierung und
experimentelle Kultivierung unter Scherströmung
vorgelegt von
M. Sc.
Philipp Jakob Kunz
an der Fakultät III – Prozesswissenschaften der Technischen Universität Berlin zur Erlangung des akademischen Grades
Doktor der Ingenieurwissenschaften – Dr.-Ing. –
genehmigte Dissertation
Promotionsausschuss:
Vorsitzende: Prof. Dr. Steffi Knorn Gutachter: Prof. Dr. Rudibert King Gutachter: Prof. Dr. Rainer Krull
Tag der Wissenschaftlichen Aussprache: 09. Juli 2021
Berlin 2021
Vorwort
Die vorliegende Arbeit entstand anlässlich meiner Tätigkeit als wissenschaftlicher Mitarbeiter am Fachgebiet derMess- und Regelungstechnik, in Kooperation mit dem Fachgebiet für Angewandte und Molekulare Mikrobiologieder Technischen Universität Berlin.
An dieser Stelle möchte ich mich zuallererst bei meinem Mentor, Herrn Prof. Dr. Rudibert King, für seine Zuversicht und sein Vertrauen bedanken, das er mir bei der Betreuung dieser Arbeit entgegengebracht hat. Bei der Bearbeitung ließ er mir viele Freiheiten, sodass diese Arbeit vor allem Ausdruck meiner persönlichen Interessen im Bereich der Ingenieurswissenschaften werden konnte. Ihm und meinen Kollegen im Fachgebiet habe ich unzählige lustige Erinnerungen und - vor allem - Denkanstöße zu verdanken.
Mein Dank gilt ebenfalls Herrn Prof. Dr. Rainer Krull, der sich bereiterklärt hat, als Zweitgut- achter dieser Arbeit zu fungieren und mit dem ich viele interessante Gespräche im RahmenAG Filamentöse Pilze in den vergangenen Jahren führen konnte. Ebenfalls danke ich Frau Prof. Dr.
Steffi Knorn für ihre Unterstützung als Vorsitzende des Prüfungsausschusses.
Mein besonderer Dank gilt auch Frau Prof. Dr. Vera Meyer, die mich in ihrem Fachgebiet An- gewandte und Molekulare Mikrobiologie willkommen geheißen, jederzeit fachlich unterstützt und mir viel Vertrauen entgegengebracht hat. Besonderen Dank auch an ihr gesamtes Team - ins- besondere an Herrn Lars Barthel - für die freundschaftliche Arbeitsatmosphäre, viele wertvolle Anregungen und seine stete Hilfsbereitschaft, die wesentlich zum Gelingen dieser Arbeit beige- tragen haben.
Unter allen Kollegen am Fachgebiet der Mess- und Regelungstechnik möchte mich zuerst bei Herrn Lennart Alker bedanken, dass er mich in den letzten Jahren mit großem Einsatz im Labor unterstützt hat. Unter den zahlreichen Freunden und Kollegen, an die ich während der Bearbeitung dieses Themas mit zahllosen Korrekturwünschen herangetreten bin, möchte ich Herrn Dr. Kars- ten Pufahl, Herrn Karl Neuhäuser, Frau Sarah Wendland, Herrn Andreas Lange, Frau Dr. Flavia Neddermeyer und Herrn Dominik Krämer für ihre Hilfe ebenfalls herzlichst danken!
Danke an die Deutsche Forschungsgemeinschaft für die Förderung dieses Projekts im Schwer- punktprogramm SPP 1934 zur Dispersitäts-, Struktur- und Phasenänderungen von Proteinen und biologischen Agglomeraten in biotechnologischen Prozessen.
Mein letzter und herzlichster Dank gilt aber Herrn Moritz Decker, für seine uneingeschränkte, liebevolle und vielseitige Unterstützung während meines Studiums, ohne die diese Arbeit so nicht möglich gewesen wäre.
Berlin, 13. Juli 2021
P. Kunz
Zusammenfassung
Der filamentöse PilzAspergillus niger zeichnet sich durch eine hohe Sekretionskapazität aus und hat besondere Bedeutung für die industrielle Herstellung von Zitronensäure sowie verschiedener Enzyme wie Amylasen, Pektinasen und Proteasen. Auf intrazellulärer Ebene ist die Produktsekreti- on das Ergebnis von Transportvorgängen, bei denen gebildete Proteine in sekretorischen Vesikeln verkapselt, zur Zellmembran transportiert und, schließlich, durch Fusion mit selbiger in den extra- zellulären Raum entlassen werden. Außerdem werden sekretorische Vesikel als Cargo-Vehikel für Zellwandmaterial genutzt, wodurch neue Zellwand aufgebaut und existierende Zellwand verstärkt werden können. Sekretorische Vesikel haben damit sowohl Einfluss auf das Wachstum als auch auf die Sekretion von intrazellulären Metaboliten und haben somit eine Doppelfunktion.
Die vorliegende Arbeit untersucht den Vesikeltransport auf zwei methodisch unterschiedlichen Wegen: Zunächst werden, im ersten Hauptteil, die sekretorischen Vesikel inAspergillus nigerma- thematisch über das Zellwachstum abgebildet. Mit diesem Ziel wird ein mechanistisches Modell aufgebaut, in dem sekretorische Vesikel entlang von Hyphen transportiert und zum Zwecke des Wachstums an der Hyphenspitze verbraucht werden. Das Modell wird beispielhaft an konfokale Messdaten mehrerer junger Keimlinge in Flüssigkultur angepasst. Die Verteilung der sekretori- schen Vesikel im Pilz wird durch ein fluoreszierendes Markerprotein intrazellulär sichtbar gemacht und zusammen mit der Gesamtlänge der Hyphe durch bildanalytische Methoden extrahiert. Die angepassten Längen- und Vesikelverläufe belegen die generelle Eignung des Modells, die Dyna- mik der Etablierung einer hohen Vesikeldichte an den Spitzen junger Keimschläuche inAspergillus niger zu beschreiben. Die Annahme, dass die effektive Transportgeschwindigkeit der sekretori- schen Vesikel zur Spitze hin abnimmt, führte zu einer Transportdynamik, welche die gemessene Akkumulation sekretorischer Vesikel im Spitzenbereich nachstellte.
Im zweiten Hauptteil dieser Arbeit wird der Einfluss des hydrodynamischen Stresses auf die sekretorischen Vesikel untersucht. Hierfür wird eine neuartige Wachstumskammer entwickelt, umAspergillus niger unter Einwirkung zeitvariabler Wandschubspannung zu kultivieren. In der Kammer kommt es hinter einer rückwärts-gewandten Stufe zur Aufscherung der Grenzschicht der Nährflüssigkeit und zur Bildung von Wirbeln. Diese führen zu zeitlich veränderlichen Wandschub- spannungen am Boden der Kammer und treten in Interaktion mit dort wachsendem Pilzmyzel. In Experimenten werden die Auswirkungen dieser Strömungsumgebung auf die Morphologie, die Verteilung und die dynamische Entwicklung der Biomasse, der sekretorischen Vesikel und eines Zielenzyms (Glucoamylase) vermessen und bildanalytisch ausgewertet. Die Wandschubspannung am Boden der Kammer mit rückwärts-gewandter Stufe wird durch CFD-Simulationen als Pseudo- Messwerte prädiziert und für vier Volumenströme des Nährmediums mit den mikroskopischen Aufnahmen verglichen. Obwohl die Quantität der gebildeten Biomasse vom Volumenstrom un- beeinflusst bleibt, bilden sich verschiedenartige morphologische Strukturen nahe der Stufe und im mittleren Bereich der Kammer aus. Während die gemessene Menge der sekretorischen Vesikel bei erhöhter Strömungsbelastung abnahm, erwiesen sich sowohl eine niedrige als auch eine hohe Stömungsbelastung als nachteilig für die Produktion von Glucoamylase.
Abstract
The filamentous fungi Aspergillus niger is well-known for its high production capacity. Its high industrial importance is due to the production of citric acid and several enzymes, such as amylases, pectinases, and proteases. On an intracellular level, product secretion is the result of transport pro- cesses which encapsulate produced proteins in secretory vesicles, move them to the cell membrane and fuse them in order to unload its contents into extracellular space. Furthermore, secretory vesic- les serve as cargo vehicles for cell-wall-related material, which is needed to build up new biomass at the fungal tips or to further strengthen the existing cell wall. Thus, secretory vesicles have an impact on both growth and the secretion of intracellular metabolites and, therefore, have a double function in the metabolism.
This thesis contributes to the mentioned field of research by investigating secretory vesicle transport in two, methodically different ways: In the first part of this thesis, secretory vesicles ofAspergillus niger and cell growth are modeled mathematically. For this purpose, a mechanistic model was constructed that describes the intracellular transport of secretory vesicles through the hyphae and their consumption for the purpose of tip growth. The model was fit to data of several young germlings, which were measured in three-dimensional space by means of confo- cal laser-scanning microscopy. The distribution of eGFP-tagged secretory vesicles was extracted through image analysis. The parameter estimation result proved the overall ability of the model to describe the dynamics of the establishment of high vesicle density at the tips of young germ tubes inAspergillus niger. Our assumption that the effective transport velocity decreases towards the tip led to transport dynamics which recreates the meassured accumulation of secretory vesicles in the tip area.
In the second part of this thesis, the influences of hydrodynamic stress on secretory vesicles, the biomass and secretion are investigated. For this purpose,Aspergillus nigerwas cultivated in a self-developed growth chamber. The flow in the chamber develops an unstable shear layer after flowing over a backward-facing step and leads to the development and transport of vortices, leading to changing wall shear stress at the bottom of the chamber. In experiments, the effects of the flow environment on the morphology, the temporal development and distribution of the biomass, the secretory vesicles, and a target enzyme (glucoamylase) were then measured and evaluated by image analysis. The wall shear stress at the bottom of the chamber was predicted by CFD simulations as pseudo-measured values and compared with the microscopic images for four different flow rates of the culture medium. Even though the quantity of the formed biomass remained uninfluenced by the flow rate, different morphological structures form near the step and in the middle area of the chamber. While the measured amount of secretory vesicles decreased with increased flow stress, both low and high flow stress proved detrimental to glucoamylase production.
Inhaltsverzeichnis
1 Einleitung 1
2 Numerische Grundlagen 5
2.1 Allgemeine Kontinuitätsgleichung . . . 5
2.2 Lösung partieller, eindimensionaler Differentialgleichungen . . . 6
2.2.1 Partielle Differentialgleichungen . . . 6
2.2.2 Allgemeine Linienmethode . . . 6
2.2.3 Finite Volumina . . . 7
2.2.4 Gewichtete, essenziell nicht-oszillierende Schemata (WENO-LF) . . . 9
2.2.5 Zeitdiskretisierung: Expliziter RK-45 . . . 12
2.3 Computergestützte Fluidsimulation . . . 15
2.3.1 Wandschubspannung, Definition . . . 15
2.3.2 Verwendete Wirbelmodelle . . . 17
2.3.3 Zeitlich gemittelte Kenngrößen . . . 18
2.3.4 Örtlich gemittelte Kenngrößen . . . 19
2.3.5 Software zur Fluidsimulation . . . 19
2.3.6 Physikalische Kenngrößen des Nährmediums . . . 20
2.4 Parameteridentifikation und -unsicherheit . . . 21
2.4.1 Optimierungsbasierte Ansätze . . . 21
2.4.2 Markov-Ketten-Monte-Carlo (MCMC) . . . 22
3 Bilddatengewinnung und -auswertung 27 3.1 Bildanalytische Grundlagen . . . 27
3.1.1 Hellfeldmikroskopie . . . 28
3.1.2 Konfokale Laser-Scanning-Mikroskopie . . . 31
3.2 Verwendete Reporterstämme . . . 32
3.3 Aufnahme-Einstellungen am CLSM . . . 34
3.3.1 Konfiguration 1: Einzelhyphen . . . 34
3.3.2 Konfiguration 2: Übersichtsaufnahmen . . . 35
3.4 Programme zur Bildanalyse . . . 36
3.4.1 Verwendete Standardoperationen der Bildanalyse . . . 36
3.4.2 Programm zurmorphologischen 3-D-Rekonstruktion von Filamenten 38 3.4.3 Programme zur Quantifizierung der Biomasse und der Fluoreszenz . . . 40
4 Versuchsaufbau 43 4.1 Auslegung und Fertigung eines Prototyps der speziellen Wachstumskammer . . . . 43
4.2 Beschreibung des Pumpsystems . . . 47
4.3 Vor- und Nachbereitung zur Kultivierung unter Strömung . . . 50
5 Mikroskalen-Modell des Wachstums von Einzelhyphen 51 5.1 Stand der Forschung . . . 51
5.1.1 Morphologische Entwicklungsphasen vonA. niger: Von der Spore zum Myzel . . . 51
5.1.2 Mikro-Skalen-Modellbildung: Stand der Forschung . . . 57
5.3 Herleitung eines Vesikel-basierten Wachstumsmodells . . . 61
5.4 Versuchsbedingungen und Stammdaten . . . 69
5.5 Parameteridentifikation . . . 69
5.6 Ergebnisse . . . 72
5.7 Diskussion und Fazit . . . 79
5.7.1 Diskussion der Ergebnisse . . . 79
5.7.2 Diskussion der Methoden . . . 81
5.7.3 Fazit und Ausblick . . . 84
6 Entwicklung einer speziellen Wachstumskammer und Kultivierung vonA. niger 87 6.1 Stand der Forschung . . . 88
6.1.1 Zellstress während der Fermentation . . . 88
6.1.2 Durchflusskammer mit rückwärts-gewandter Stufe: Stand der Forschung . . 91
6.2 Strömungssimulation . . . 96
6.2.1 Innenraumgeometrie . . . 96
6.2.2 Rechengitter . . . 97
6.2.3 Rand- und Anfangswerte . . . 101
6.2.4 Ergebnisse der Simulation (Stufenzone, LES) . . . 106
6.3 Strömungsversuche in der speziellen Wachstumskammer . . . 115
6.3.1 Stammdaten und Probenvorbereitung . . . 115
6.3.2 Vorversuch zur Widerstandsfähigkeit adhärenter Sporen gegen Strömung . 115 6.3.3 Inkubation, Strömungsverlauf und Bildaufnahme . . . 116
6.4 Ergebnisse . . . 119
6.4.1 Übersichtsaufnahmen und Metastrukturen . . . 119
6.4.2 Dreidimensionale Verteilung . . . 121
6.4.3 Biomasse . . . 123
6.4.4 Sekretorische Vesikel (SV-eGFP) . . . 125
6.4.5 Glucoamylase (GlaA-dTomato) . . . 127
6.5 Diskussion und Fazit . . . 128
6.5.1 Diskussion der Ergebnisse . . . 128
6.5.2 Diskussion der Methoden . . . 132
6.5.3 Fazit und Ausblick . . . 136
7 Zusammenfassung dieser Arbeit 139 A Materialanhang 151 A1 Koordinatentransformation . . . 151
A2 WENO-Schema: WENO35-LF . . . 151
A3 Zeitlich gemittelte Wandschubspannungen, Ergebnisse der CFD-Simulation . . . . 153
A4 Konstruktionszeichnungen . . . 159
A4.1 PPFC mit BFS . . . 159
A4.2 Verworfene Entwürfe . . . 163
Symbolverzeichnis
Abkürzungen / Akronyme Abk. Beschreibung
AR Erstreckung (aspect ratio)
BFS rückwärts-gewandte Stufe (backward-facing step)
CFD Computer-gestützte Fluiddynamik (computational fluid dynamics) CFL Courant-Friedrichs-Lewy-Zahl
CLSM konfokales Laser-Scanning Mikroskops (confocal laser-scanning microscope) DNS direkte numerische Simulation
eGFP verbessertes, grün-fluoreszierendes Protein (enhanced green fluorescent protein) ER relative Querschnittserweiterung (expansion ratio)
GlaA Glucoamylase
LES Large-Eddy-Simulation
OPSL optically pumped semiconductor laser
PPFC parallele Plattenströmungskammer (parallel plate flow chamber) PT parallel tempering
RANS Reynolds-gemittelte Simulation (reynolds-averaged numerical simulation) SNARE solubale N-ethylmalemide sensitive factor attachment protein receptor SV sekretorische Vesikel
WALE Wand-adaptive lokale Wirbelviskosität (wall-adapting local eddy-viscosity) WENO gewichtete, essentiell nicht oszillierende Schemata (weighted essentially non-
osciallatory scheme)
WSS Wandschubspannung (wall shear stress) Lateinische Buchstaben, klein
Buchstabe Beschreibung Einheit Kap./Abs.
𝒾,𝒿,𝓀 Laufindizes
m,n,o Laufindizes für diskrete Bilddatenstapel
𝑤˜(1,2,...) Gewichtungsfaktor, WENO 2.2.4
mSkelett Vektor aus𝑥-Positionen der Mittellinie (Skelett) 3.4.2
nSkelett Vektor aus𝑦-Positionen der Mittellinie (Skelett) 3.4.2
oSkelett Vektor aus𝑧-Positionen der Mittellinie (Skelett) 3.4.2
𝑛 Normalenvektor 2.1
𝑢 Geschwindigkeit, vektoriell m s−1 2.3.1
𝑑(1,2,...) Lineargewicht, WENO 2.2.4
Φ
𝑓IC Faktor für Anfangsbedingungen 5.6
𝑓K Körperkräfte, allgemein N 2.3.1
𝑓LB Abklingfaktor des linken Randwertes min−1 5.3
𝑓sl Instabilitätsfrequenz s−1 6.2.3
ℎ freie Höhe des Kanals μm 6.1.2
ℎRK Zeitschrittweise bei Runge-Kutta-Integration s 2.2.5
𝑗𝑖 Flussdichte kg s−1μm−25.3
𝑘C maximale Geschwindigkeit der Spitzenerweiterung min−1 5.3
𝑙thres Schwellenwert bei Binarisierung, Bildanalyse 3.4.2
𝑚 Masse kg 5.3
𝑚V Masse sekretorischer Vesikel kg 5.3
𝑚MT Masse Mikrotubuli kg 5.3
𝑝 Druck Pa 2.3.1
𝑞∗ Kandidaten-Verteilung, MCMC 2.4.2
𝑟 Hyphenradius μm 5
𝑠 Laufkoordinate entlang der Mittellinie 3.4.2
𝑠 Stufenhöhe μm 6.1.2
𝑡 Zeit s 2.1
𝑡0 Anfangszeitpunkt s 2.2.1
𝑡end Endzeitpunkt s 2.2.1
𝑢(𝑥) Geschwindigkeitsfeld im 1-D-Hyphenmodell μm min−1 5.3
𝑢∞ Geschwindigkeit der Kernströmung, CFD m s−1 6.1.2
𝑢max maximale Geschwindigkeit in der Hyphe μm min−1 5.3 𝑢min,∞ minimale Geschwindigkeit an der Hyphenspitze μm min−1 5.3
𝑣 normierte Konzentration sekretorischer Vesikel 5.3
𝑤(1,2,...) normierter Gewichtungsfaktor, WENO 2.2.4
𝑥 erste kartesische Koordinate μm 6.1.2
𝑥r Wiederanlegelänge bzw. -linie μm 6.1.2
𝑦 zweite kartesische Koordinate μm 6.1.2
𝑧 dritte kartesische Koordinate μm 6.1.2
𝓅 Wahrscheinlichkeit, Bayes 2.4
Lateinische Buchstaben, groß
𝑉 Volumenstrom mL min−1 6.2.3
A Entscheidungskriterium 2.4.2
F Fluoreszenz-Kanal 6.3
G grüner Fluoreszenz-Kanal 6.3
M Modell, allgemein, MCMC 2.4.2
O Konvergenzordnung 2.4.2
P A-posteriori-Verteilung 2.4.2
R roter Fluoreszenz-Kanal 6.3
S großer Rechenstern 2.2.4
T Transmissions-Kanal 6.3
V Konzentration sekretorische Vesikel kg m−3 5.3
Y Messwerte, allgemein, MCMC 2.4.2
I 3-D-Bilddaten-Matrix der FormI[𝑚, 𝑛, 𝑜, 𝑡𝑖] 3.1
St Strouhal-Zahl 6.2.3
M Bildbreite, diskret Pxl 3
N Bildhöhe, diskret Pxl 3
O Gesamttiefe 3-D-Bilddatenstapel, diskret Vxl 3
𝐶1,2 Mittellinie in zwei Raumdimensionen μm 3.4.2
𝐴 Querschnittfläche µm2 5.3
𝐷 Diffusionskonstante μm2min−1 5.3
𝐸{·} Erwartungswert 2.4.2
𝐹(𝜙) Fluss U𝜙 s−1 2.2.3
𝐻 Gesamthöhe des Kanals μm 6.1.2
𝐾C empirische Sättigungskonstante 5.3
𝐿 Hyphenlänge μm 5.3
𝐿0 Anfangslänge des simulierten Keimlings μm 5.5
𝐿ref Referenzlänge für𝑅𝑒 μm 6.1.2
𝑁C Anzahl Knoten inΔ𝑥C 2.2.3
𝑁it Anzahl an Iterationen bzw. Anzahl der akzeptierten Parameter-Sets, MCMC
5.6
𝑁rot Drehzahl rad min−1 6.1.1
𝑁T Anzahl Läufer in Parallel Tempering, MCMC 2.4.2
𝑁𝑡 Anzahl Zeitschritte bzw. Frames (Bildanalyse), diskret 3
𝑁𝑥 Anzahl Knoten in 1-D-Finite-Volumen-Gitter 2.2.3
𝑁grau Anzahl der Grauwerte im Digitalbild 42
𝑁mult Anzahl der Startwert-Vektoren, Multi-Startwerte- Algorithmus
2.4.1
𝑃(1,2,...) Teil-Polynome, WENO 2.2.4
𝑄(𝑥) Groß-Polynom, WENO 2.2.4
𝑆𝜙 Quellterm für𝜙, allgemein 2.1
𝑋 Bildbreite μm 3
𝑌 Bildhöhe μm 3
𝑌𝐿 Längenausdehnungsfaktor μm−2 5
𝑍 Gesamttiefe 3-D-Bilddatenstapel μm 3
Re Reynolds-Zahl 6.1.2
Griechische Buchstaben, klein
𝛼 Einstellparameter für Parallel Tempering 2
𝛽RK Sicherheitsfaktor bei Runge-Kutta-Integration 2.2.5
𝛽(1,2,...) lokaler Glattheits-Indikator 2.2.4
𝛿 Dicke der laminaren Grenzschicht μm 6.1.2
𝛿∗ Verdrängungsdicke m 6.2.3
𝜖 Null-Vermeidungs-Faktor 2.2.4
𝜖Gitter relativer Fehler der Lösung im Rechengitter mit 𝑁KnotenKnoten
6.2.2
𝜂 tranformierte𝑥-Koordinate – 5.3
𝜆 Wellenlänge nm 3.2
𝜇 dynamische Viskosität Pa s 6.1.2
𝜇C Konsumterm m−3min−1 5.3
𝜇P Produktionsterm m−3min−1 5.3
𝜇𝑡 Wirbelviskosität Pa s 2.3.1
𝜙 intensive Größe, allgemein UΦm−3 2.1
𝜌 Dichte kg m−3 2.3.6
𝜎MCMC Parameterunsicherheit, errechnet aus MCMC-Kette 2.4.2 𝜏w,𝑖 erste Komponente des Wandschubspannungsvektors Pa 2.3.1 𝜏w,𝑗 zweite Komponente des Wandschubspannungsvektors Pa 2.3.1
𝜏𝑖𝑖 𝑖-te Normalenspannung von𝜏 Pa 2.3.1
𝜏𝑖 𝑗 Schubspannungskomponente von𝜏 Pa 2.3.1
𝜃P Produktionsrate min−1 5.3
𝜇 Mittelwertvektor für𝜃-dimensionale Normalverteilung, Parallel Tempering
2
𝜏w Wandschubspannung Pa 2.3.1
𝜃 Parametersatz, vektoriell 2.4
𝜏 Stresstensor 2.3.1
Griechische Buchstaben, groß
Δ𝑡sim Zeitschrittweite in der LES s 6.2.3
Δ𝑡tot,𝒾 Totzeit-Term des𝒾-ten Keimlings min 5.2
Δ𝑥 Abstand von Gitterknoten bzw. Zellbreite μm 2.2.3
Δ𝑥C Verbrauchszone von SV hinter der Spitze μm 5.3
Δ𝑥P Produktionszone von SV hinter der Spitze μm 5.3
Δ𝑥𝑢 Geschwindigkeits-Abstiegs-Zone μm 5.3
Δ𝑧 Abstand, axial μm 5
Φ extensive Größe, allgemein UΦ 2.1
Σ Kovarianzmatrix für𝜃-dimensionale Normalverteilung, Parallel Tempering
2
Ξ Zielfunktion, allgemein 5.5
Zahlenmengen, Operatoren, Verteilungen
ℝ reeller Zahlenraum 2.2
N Normalverteilung 2.4.2
U Gleichverteilung, stetig 2.4.2
∇ Nabla-Operator 2.3.1
Indizes und Modifikatoren
(¯¯·) zeitliche Varianz 2.3.3
(·)¯ zeitliche Mittlung 2.3.3
(·)ˆ geschätzter Wert 5.5
(·) transformierte Koordinate, front-fixing-Methode 5.3
(·)∗ Kandidat von (·), z. B.𝜃∗ 2.4.2
(·)−1 invertierte Intensität, Bildanalyse 3.4.3
(·)F feinste Lösung, CFD 6.2.2
(·)in zufließende Größe 5.3
(·)L linker Randwert 2.2.1
(·)out abfließende Größe 5.3
(·)R rechter Randwert 2.2.1
(·)𝑖 erste, d. h.𝑥-Komponente des Vektors 6.1.2
(·)𝑗 zweite, d. h.𝑦-Komponente des Vektors 6.1.2
(·)𝑘 dritte, d. h.𝑧-Komponente des Vektors 6.1.2
(·)𝑧 Mittlung über Slices in Stack in𝑧-Richtung, Bildanalyse 3.4.3
(·)max Maximalwert 5.3
(·)min Minimalwert 5.3
(·) zeitliche Ableitung 5.3
(·)0 Anfangswert, z. B.𝑡0 2.2.1
· räumlicher Mittelwert 2.3.4
· räumliche Varianz 2.3.4
(·)˜ transformierte Koordinate, Bildanalyse 3.4.2
1 Einleitung
Myzelförmig wachsende Mikroorganismen spielen eine entscheidende Rolle in der biotechnolo- gischen Industrie für die Produktion verschiedener Enzyme, Antibiotika, organischer Säuren und Aromastoffe. Pilze (lat.: Fungi) aus der Familie der Aspergilliproduzieren u. a. Enzyme, die zur Verarbeitung von Stärke in der Back-, Brauerei- sowie der Getränkeindustrie eingesetzt werden.
Der schwarze GießkannenschimmelAspergillus nigerweist die höchste Sekretionskapazität inner- halb dieser Gattung auf, worin sich sein besonderer industrieller Nutzen begründet. Neben seinen herausragenden Eigenschaften bei der Kompostierung von Abfällen wirdAspergillus niger groß- technisch für die Herstellung von Zitronensäure, Glukonsäure und Glucoamylase genutzt [3].
Die industrielle Fermentation von A. niger wurde über Jahrzehnte hinweg in Bezug auf ihre Produktivität optimiert [116]. Während filamentöse Pilze in der Natur oft auf festen Substraten gedeihen, werden sie in der Industrie meistens in Flüssigkultur in Bioreaktoren gezüchtet. Diese industrielle Art der Kultivierung führt – je nach gewählten Prozessbedingungen – zu dispergier- tem Myzel oder zu kugelförmigen Mikrokolonien, auch Pellets genannt. Der Effekt von morpholo- gischen Veränderungen auf die Produktionsleistung ist inA. niger unbestritten groß, allerdings wurde bislang keine klare Beziehung zwischen Morphologie und Produktivität gefunden [54].
Prozessparameter, die die Morphogenese filamentöser Organismen in Reaktorfermentationen be- einflussen, sind vielfältig und reichen von der Nährstoffzusammensetzung und der Temperatur des Mediums, der Konzentration des Gelöstsauerstoffs bis zur Reaktorgeometrie – und damit zur Fluiddynamik. In Kulturbrühen mit komplexen, makroskopischen Myzelien ändern sich organis- menbedingt die rheologischen Eigenschaften des Mediums, was den Masse- und den Gastransport zusätzlich erschwert. Einige Studien finden zwar Relationen zwischen Prozessbedingungen und Produktausbeuten, aber die Ursachen dieser Zusammenhänge, d. h. die Auswirkungen der genann- ten Prozessparameter auf zellphysiologischer Ebene, bleiben weitestgehend unbekannt [151, 153, 121]. Da inAspergillusdie Produktion und die Sekretion von Metaboliten von den intrazellulären Transportvorgängen abhängen, kann ein tieferes Verständnis der metabolischen Zusammenhän- ge zwischen Nährstoffaufnahme, Zellwachstum und Sekretion diese Lücke schließen und wichtige Impulse für das Design zukünftiger Kultivierungsprozesse liefern.
Als essenzielle Organelle spielen sekretorische Vesikel (Abk.: SV) eine entscheidende Doppelrol- le, sowohl für den intrazellulären Transport, die Sekretion als auch für das Wachstum: Zum einen transportieren SV Enzyme und wichtige Metabolite zu den Hyphenspitzen des Pilzes, wo diese mit der Zellmembran fusionieren und den Cargo an die Umgebung abgeben (Exozytose). Hierbei wird umgebendes organisches Material zersetzt und für die anschließende Aufnahme in den Pilz (Endo- zytose) konditioniert. Zum andern liefern SV auch Grundbestandteile zum Aufbau neuer Zellwand und ermöglichen so das apikale Wachstum in filamentöser Form zu sog. Hyphen. Die durch SV gelieferten Baustoffe dienen so dem Erhalt und der Verdichtung der Zellwand und damit der Si- cherung der Integrität der Zelle. Durch die moderne konfokale Laser-Scanning-Mikroskopie (Abk.:
CLSM) und die gezielte, genetische Veränderung können SV (und andere Organelle) mit Fluores- zenzproteinen markiert und ihre dreidimensionale Verteilung im Pilz sichtbar gemacht werden
𝑦 𝑥
Spore
Hyphenspitze
Abbildung 1.1:Aufnahme einer jungen Hyphe nach der Keimung. Grün fluoreszierende Marker- Moleküle an den sekretorischen Vesikeln weisen auf Akkumulation in der Spore und an der Hyphenspitze hin. Die dreidimensionale Aufnahme wurde auf die 𝑥𝑦-Ebene durch Maximum- Projektion abgebildet.
(siehe Abbildung 1.1).
In dieser Doppelrolle haben SV unmittelbaren Einfluss auf Fermentationen in der industriel- len Praxis. In mechanisch gerührten Reaktorsystemen wirdA. niger lokalen Kräften ausgesetzt, welche durch Fluid-Wirbel-Interaktionen, durch Interaktionen mit der Reaktorwand, dem Rühr- erblatt oder anderem Pilzmyzel entstehen. Der Rührvorgang stellt die Sauerstoffversorgung der Pilzbrühe sicher, übt allerdings auch erhöhten hydrodynamischen Stress auf Myzelien vonA. niger aus. Morphologische Veränderungen als Folge scherreicher Kultivierung reichen von einer Ver- dickung oder Verhärtung der Zellwände bis zu veränderter Verzweigung des Myzels [117, 115].
Diese Veränderungen erhöhen – aus Sicht des Pilzes – die Widerstandsfähigkeit gegenüber dieser externen Krafteinwirkung, können aber auch rein mechanisch bedingt durch Bruch-Aggregations- Interaktionen entstehen. Die veränderte morphologische Zusammensetzung der Pilzbrühe hat oft eine Abnahme der Produktionsleistung zur Folge [47]. Eine Verschiebung der Beladung der SV zugunsten von Zellwandbestandteilen bei äußerer Krafteinwirkung könnte die Transportkapazi- tät von Produkten mindern und damit ursächlich für dieses Phänomen sein. Diese wird hier als Hypothese im Rahmen eines Forschungsprojektes im Schwerpunktprogramm 1934 der Deutschen Forschungsgemeinschaft untersucht.
Die beschriebene Doppelrolle unterstreicht den Umstand, dass SV entscheidend für das Ver- ständnis der Sektretionsleistung von filamentösen Pilzen sind. Sie untermauert außerdem die Notwendigkeit, die Rolle der SV nicht losgelöst vom Kontext realer Fermentationen zu betrachten.
Gleichzeitig stellt sich allerdings die praktische Frage, wie SV-Transport- und Wachstumsvorgänge im Labormaßstab vermessen, und, wie der Einfluss von definierter Scherströmung auf das Myzel- wachstum quantifiziert werden können.
Um beiden Problemstellungen zu begegnen, wurde diese Arbeit in zwei, voneinander unabhän- gige Hauptkapitel unterteilt: Zunächst wird SV-Transport im ersten Hauptkapitel unter synthe- tischen Umgebungsbedingungen (d. h. ohne hydrodynamische Belastung) auf einer Mikroskala modelliert. Anschließend, im zweiten Hauptteil, wird die zeitliche Entwicklung und Verteilung der SV unter kultivierungsähnlichen Bedingungen, d. h. unter Scherstress auf einer Makroska- la untersucht. Für die Scherstressuntersuchungen wurde eine durchflossene Wachstumskammer
hergestellt, in der sowohl laminare als auch transiente Strömung mit und ohne Wirbelbildung er- zeugt werden kann. In beiden Hauptteilen findet jeweils ein Abgleich zwischen Messungen und Simulationen statt.
Gliederung dieser Arbeit
Zur Vorbereitung der beiden Hauptkapitel werden inKapitel 2die notwendigen Voraussetzun- gen für die mathematische Modellbildung, die CFD-Simulation und die Parameteridentifikation geschaffen.Kapitel 3fasst anschließend die verwendeten mikroskopischen Methoden (Hellfeld- Mikroskopie und CLSM) und die dazugehörigen optischen Grundlagen zusammen. Das Kapitel beinhaltet zudem zwei Programme zur Bildanalyse, mit denen aus vieldimensionalen CLSM- Messdaten Informationen für die beiden Hauptkapitel extrahiert werden können. InKapitel 4wird der Aufbau des Experiments zur Kultivierung vonA. niger unter Strömung erörtert. Hierfür wird zunächst die Dimensionierung der entwickelten Wachstumskammer mit rückwärts-gewandter Stufe, dann das Druck-getriebene Pumpsystem beschrieben. Die für die Probenvorbereitung, den Anschluss der Kammer und die Reinigung des Pumpsystems nötigen Schritte werden außerdem erklärt.
Imersten Hauptteil,Kapitel 5wird anschließend die elementare Rolle der SV an der Entwick- lung vonA. niger auf derMikroskalauntersucht. Hierfür wird zunächst der Stand der Forschung zum SV-Transport inA. niger und zur Mikroskalen-Modellbildung erfasst und, im Anschluss, ein mechanistisches Modell zur Beschreibung des Wachstums einzelner Keimschläuche formuliert.
Für die Validierung des Modells werden dreidimensionale Bilddaten vonA. niger-Keimlingen mit einem CLSM aufgenommen sowie die Verteilung der SV im Pilz durch bildanalytische Metho- den zugänglich gemacht. Dieser Teil der Arbeit entwickelt das in der vorausgegangenen Veröf- fentlichung von Kunz et al. [79] beschriebene Modell weiter, umfasst eine größere Anzahl von Keimlingen und dient damit dem Verständnis des Wachstums und der Transportkinetik sekretori- scher Vesikel. Es lässt aber zunächst den in Kultivierungen vorkommenden Scherstress außer Acht.
Imzweiten Hauptteil,Kapitel 6dieser Arbeit wird auf den kultivierungsbedingten Scherstress eingegangen. Zu diesem Zweck wurde eine durchflossene Kammer hergestellt, in der Hyphenpilze im Labormaßstab unter Scherströmung herangezogen werden können. Das Kapitel beginnt mit ei- ner Übersicht über den Stand der Forschung zu Zellstress in Fermentationen und weiht Leser*innen in die fluiddynamischen Hintergründe dieser speziellen Strömungsgeometrie ein. Die Fluidströ- mung und die entstehende zeitliche Variation der örtlichen Wandschubspannung wurden über CFD-Simulationen quantifiziert. In Experimenten wurde die Reaktion des am Kammerboden wach- senden PilzesA. niger auf den extrazellulären Scherstress anschließend für vier verschiedene Be- lastungsniveaus ermittelt. Die CLSM-Messungen ermöglichten Einblicke in die Morphologie, in die Verteilung und die zeitliche Entwicklung der Biomasse, der SV und des Zielenzyms Glucoamylase.
Die dynamischen Verläufe werden der prädizierten Wandschubspannung aus CFD-Simulationen gegenübergestellt. Darauf folgt eine Diskussion der Ergebnisse, der verwendeten Methoden und ein Ausblick auf den zukünftigen Nutzen der Wachstumskammer im Forschungsprojekt.
2 Numerische Grundlagen
In dieser Arbeit kommen unterschiedliche numerische Methoden zur Anwendung, deren Grund- lagen in diesem Kapitel beschrieben werden. Hierbei werden Flüsse von Masse und Impuls durch räumliche Bilanzgebiete bilanziert, weshalb, zunächst, der Fluss einer allgemeinen, extensiven Grö- ße in Unterkapitel 2.1 mathematisch formuliert wird. Zur Lösung der im ersten Hauptkapitel die- ser Arbeit aufgestellten partiellen Differentialgleichung wurde ein numerisches Verfahren pro- grammiert, das sich aus mehreren Einzelbauteilen zusammensetzt. Unterkapitel 2.2 fasst hierzu die Grundbausteine zur Lösung eindimensionaler, partieller Differentialgleichungen mittels finiter Volumenapproximation zusammen und beschreibt das Verfahren zur Ortsdiskretisierung. Für die numerische Fluidsimulation von Medium durch die entwickelte Wachstumskammer im zweiten Hauptkapitel werden anschließend die Grundlagen der CFD-Simulation erläutert und die nötigen statistischen Kenngrößen zu wandnaher Scherströmung definiert (Unterkapitel 2.3). Dieses Kapitel schließt mit der Beschreibung der verwendeten Methoden zur Parameteridentifikation in Unter- kapitel 2.4 ab.
2.1 Allgemeine Kontinuitätsgleichung
Der Fluss einer extensiven Größe,Φ, durch ein Bilanzvolumen,𝑉, kann durch eine volumenbezo- gene, d.h. intensive, hinreichend glatte Zustandsgröße,𝜙, mit𝜙 :ℝ×ℝ+0 →ℝund der allgemeinen Einheit [𝜙] =UΦm−3beschrieben werden. Damit ist𝜙 üblicherweise eine Dichte, also z. B. eine Massen-, Impuls- oder Energiedichte. Die zeitliche Änderung dieser Größe kann in differenzieller Form
𝜕
𝜕𝑡𝜙(𝑥, 𝑡) + ∇𝑓(𝜙(𝑥, 𝑡))=𝑆𝜙 (2.1) mit∇𝑓 =
𝑑𝑓𝑑𝑥, . . .
bzw. in integraler Form 𝑑
𝑑𝑡
∫
𝑉(𝑡)𝜙·𝑑𝑉 +
∫
𝑉(𝑡)∇𝑓(𝜙(𝑥, 𝑡)) ·𝑑𝑉 =
∫
𝑉(𝑡)𝑆𝜙 ·𝑑𝑉 (2.2)
angegeben werden.𝑓 bezeichnet eine gegebene Flussfunktion, während𝑆𝜙 Quell- und Senkterme des Systems beinhaltet. Gl. 2.2 vereinfacht sich durch den Gauß’schen Integralsatz zu
𝑑 𝑑𝑡
∫
𝑉(𝑡)𝜙·𝑑𝑉 +
∮
𝐴𝑓(𝜙(𝑥, 𝑡)) ·𝑛·𝑑𝐴=
∫
𝑉(𝑡)𝑆𝜙 ·𝑑𝑉 (2.3) d. h. die zeitliche Veränderung von𝜙pro Kontrollvolumen kann durch die Integration der in das Kontrollvolumen eintretenden Flüsse durch die Randflächen𝐴und der Quell- bzw. Senkterme𝑆𝜙
beschrieben werden.
2.2 Lösung partieller, eindimensionaler Differentialgleichungen
Die numerische Integration von partiellen Differentialgleichungen erfolgt in dieser Arbeit nach dem Lösungsansatz der allgemeinen Linienmethode (engl.: method of lines), der Orts- und Zeit- diskretisierung getrennt voneinander behandelt. Im folgenden Abschnitt werden eindimensionale partielle Differentialgleichungen zunächst mathematisch erfasst (Abschnitt 2.2.1). Daraufhin wird das methodische Vorgehen der Linienmethode zur Übersicht erläutert und die Ortsdiskretisierung sowie der Transport über Zellen finiter Volumina beschrieben (Abschnitt 2.2.3). Als Verfahren zur Approximation der Variablen auf den Zellengrenzen wird das WENO-Verfahren ausführlich be- handelt (gewichtete, essenziell nicht-oszillierende Schemata, Abschnitt 2.2.4).
2.2.1 Partielle Differentialgleichungen
Im Rahmen der Modellierung vieler Phänomene aus Natur- und Ingenieurwissenschaften stellt der Umgang mit zeitabhängigen, partiellen Differentialgleichungen eine herausfordernde Aufgabe dar.
Aus komplexer werdenden, zunehmend realitätsnahen Modellen ergeben sich schnell hochgradig nichtlineare, gekoppelte Systeme, in denen eine Differentialgleichung von mehreren unabhängigen Variablen abhängt. Da das in dieser Arbeit behandelte mathematische Modell dieser Problemklasse angehört, wird hier die Notation kurz eingeführt.
Diese Arbeit behandelt nichtlineare, partielle Differentialgleichungen in einer Raumdimension bis zu zweiter Ordnung (in𝑥) der allgemeinen Form
𝐺(𝑥, 𝑡, . . . , 𝜙,𝜕𝜙/𝜕𝑥,𝜕2𝜙/𝜕𝑥2,𝜕𝜙/𝜕𝑡) =0, (2.4) beziehungsweise explizit zur Zeitableitung
𝜕𝜙/𝜕𝑡 =𝑔(𝑥, 𝑡, 𝜙,𝜕𝜙/𝜕𝑥,𝜕2𝜙/𝜕𝑥2), (2.5) im Zeitintervall
𝑡 ∈Ω𝑡 :=[𝑡0, 𝑡𝑁𝑡] ⊂ℝ1 (2.6)
sowie im Ortsgebiet
𝑥 ∈Ω𝑥 :=[𝑥L, 𝑥R] ⊂ℝ1. (2.7)
Das Problem gilt dann als korrekt gestellt ("well posed"), wenn seine Lösung𝜙 existiert, ein- deutig ist und kontinuierliche Rand- und Anfangsbedingungen definiert sind. Im in dieser Arbeit betrachteten Fall eines Konvektions-dominierten Transportmodells ist das System mit den An- fangsbedingungen
𝜙(𝑥, 𝑡0) =𝜙0(𝑥), 𝑡 =𝑡0 (2.8)
und beidseitigen örtlichen Randbedingungen vollständig definiert.𝜙(𝑥, 𝑡) heißt Lösung, wenn𝜙 Gleichung 2.4 in allen Punkten𝑥 ∈Ω𝑥 und für alle Zeiten𝑡 ∈Ω𝑡 erfüllt.
2.2.2 Allgemeine Linienmethode
Im Kern basiert die Idee der Linienmethode auf der Diskretisierung aller örtlichen Variablen, um ein System partieller (Abschnitt 2.2.1) in ein System gewöhnlicher Differentialgleichungen zu transfor-
2.2 Lösung partieller, eindimensionaler Differentialgleichungen
mieren. Das ursprüngliche Randwertproblem wird damit in ein Anfangswertproblem erster Ord- nung überführt, welches mit Ein- oder Mehrschritt-Zeit-Integratoren numerisch gelöst werden kann. Abbildung 2.1 stellt das Grundprinzip der Linienmethode dar.
semi-diskrete DGL
𝜕𝜙(𝑥𝒾, 𝑡)
𝜕𝑡 =𝑔(𝑥, 𝑡, 𝜙(𝑥)) partielle DGL 𝜕𝜙(𝑥, 𝑡)𝜕𝑡 =
𝑔(𝑥, 𝑡, 𝜙, 𝜕𝜙𝜕𝑥, 𝜕𝜕𝑥2𝜙2)
diskrete Lösung 𝜙(𝑥𝒾, 𝑡) = 𝜙𝒾(𝑡)
Gittergenerierung; Diskretisierung
zeitliche Integration
Abbildung 2.1:Grundprinzip der Linienmethode
Bei der Linienmethode werden klassischerweise folgende Schritte durchlaufen: Beginnend mit der Gittergenerierung (engl.: gridding) wird das Lösungsgebiet in diskrete Gitterpunkte aufgeteilt.
Danach wird das Gleichungssystem auf dem Rechengitter örtlich diskretisiert und damit in ein semi-diskretes System von𝑁𝑥 Differentialgleichungen
𝜕𝜙(𝑥𝒾, 𝑡)
𝜕𝑡 =𝑔𝒾(𝑥, 𝑡, 𝜙(𝑥, 𝑡)) mit𝒾=1,2, . . . , 𝑁𝑥 (2.9) mit𝑥 = [𝑥1, . . . , 𝑥𝒾, . . . , 𝑥𝑁𝑥]𝑇 überführt. Die in dieser Arbeit verwendete Methode der finiten Vo- lumina wird in Abschnitt 2.2.3 vorgestellt. Die Berechnung nach der finiten Volumen-Methode erfordert hierbei die Approximation der ortsverteilten Variablen an den Zellgrenzen durch das WENO-Interpolationsschema, welches in Abschnitt 2.2.4 eingeführt und dessen Anwendung auf Erhaltungsgleichungen demonstriert wird. Im letzten Schritt erfolgt die Lösung des Differential- gleichungssystems mit einem geeigneten numerischen Integrator. In dieser Arbeit findet ein expli- zites Runge-Kutta-Schema, siehe Abschnitt 2.2.5, Verwendung.
2.2.3 Finite Volumina
Erhaltungsgleichungen werden üblicherweise über numerische Verfahren angenähert, welche die Lösung auf einem diskreten örtlichen Rechengitter approximieren. In dieser Arbeit wird die Me- thode der Zellen-zentrierten finiten Volumina1verwendet, welche aus der integralen Formulierung der physikalischen Erhaltungsgleichungen, Gl. 2.2 resultiert. Das eindimensionale, örtliche Bilanz- gebiet gem. Gl. 2.7 der Länge,𝐿=𝑥R−𝑥L, wird hierbei in eine endliche Anzahl kleinerer Teilgebiete – sog. Volumina oder Zellen – aufgeteilt, wie in Abbildung 2.2 dargestellt. Die Bilanzierung einer extensiven, als global glatt angenommenen Größe𝜙 erfolgt pro Zelle𝒾am Ortspunkt des jeweili- gen Zellmittelknotens𝑥𝑖
𝜙𝒾(𝑡) 𝜙(𝑥𝒾, 𝑡), mit𝒾=1,2, . . . , 𝑁𝑥 . (2.10)
1Im Gegensatz zu Knoten-zentrierten finiten Volumina.
L R
1 2 . . . 𝒾 . . . 𝑁𝑥−1 𝑁𝑥
Zellmittelknoten𝑥𝒾, mit𝒾=
𝐿 Δ𝑥
1
2 1+12 𝒾−12 𝒾+12 𝑁𝑥−12 𝑁𝑥+12
Zellgrenzen𝑥𝒾, mit𝒾=
Abbildung 2.2:Rechengitter des Zellen-zentrierten finiten Volumen-Verfahrens.
Das Rechengitter sei äquidistant, d. h. es giltΔ𝑥 =𝑥𝒾+1−𝑥𝒾 =konst., mit
𝑥1<𝑥2<· · · <𝑥𝑁𝑥−1 <𝑥𝑁𝑥 . (2.11) Alle Zellen werden durch Zellränder begrenzt, die explizit nicht Teil der angrenzenden Zellen sind, d. h. für die zunächst keine der Werte auf den angrenzenden Zellmittelknoten gelten. Die Zellengrenzen werden in Halbschritten als Index geführt, d. h., die Zellengrenze zwischen den Zellmittelknoten𝑥𝒾 und𝑥𝒾+1wird als𝑥𝒾+1
2 bezeichnet. Für𝜙 gilt an den Zellengrenzen entspre- chend
𝜙𝒾(𝑡) 𝜙(𝑥𝒾, 𝑡), 𝒾 = 1 2,3
2, . . . , 𝑁𝑥 +1
2 . (2.12)
Der Nettofluss eines Kontrollvolumens, d. h.𝜙in−𝜙out, ist gleich der Summe der Integrale der Flächen-spezifischen Flüsse𝑓 über alle vorhandenen Seitenflächen𝑆𝑘
∮
𝐴𝑓(𝜙(𝑥, 𝑡)) ·𝑛·𝑑𝐴=
𝑘
∮
𝑆𝑘
𝑓(𝜙(𝑥, 𝑡)) ·𝑛·𝑑𝐴 . (2.13) Im eindimensionalen Fall sind – mit Ausnahme der Randzellen – für𝑥𝒾zwei Grenzflächen in𝑥𝒾−1/2
und in𝑥𝒾+1/2vorhanden. Die Flächen-spezifischen Flüsse𝑓 durch die beiden Randflächen der𝒾-ten Zelle vereinfachen sich in diesem Fall zu
∮
𝐴𝒾±1/2
𝑓(𝜙𝒾) ·𝑛·𝑑𝐴=𝑓(𝜙𝒾±1/2) ·𝐴𝒾±1/2=±𝐹(𝜙𝒾±1/2), (2.14) mit dem flächenspezifischen Fluss𝑓(𝜙) bzw. dem Fluss𝐹(𝜙). Für das Volumenelement𝒾verein- facht sich die Erhaltungsgleichung Gl. 2.3, S. 5 damit zu
𝑑 𝑑𝑡
∫
𝑉𝜙𝒾·𝑑𝑉 +
𝐹(𝜙𝒾+1/2) −𝐹(𝜙𝒾−1/2) =
∫
𝑉𝑆𝜙·𝑑𝑉 . (2.15)
Die Flüsse auf den Rändern,𝐹(𝜙𝒾±1/2), werden durch die sie umgebenden Knotenwerte mittels ver- schiedener Diskretisierungsverfahren als sog. numerische Flüsse ˆ𝐹(𝜙𝒾±1/2) rekonstruiert. Das zu diesem Zweck in dieser Arbeit verwendete WENO-Verfahren wird nachfolgend genauer erläutert.
2.2 Lösung partieller, eindimensionaler Differentialgleichungen
2.2.4 Gewichtete, essenziell nicht-oszillierende Schemata (WENO-LF)
Folgende Ausführungen zu gewichteten, essenziell nicht-oszillierenden Rechenschemata (Abk.
WENO2) basieren auf den Arbeiten von Shu [135] und dienen der Herleitung des Interpolations- schemas.
Es ist bekannt, dass die klassischen, linearen Diskretisierungsverfahren im Falle steiler Gradi- enten zu Über-, Unterschwingungen oder, allgemein, zu nicht-physikalischem numerischem Ver- halten neigen [41]. Um diese Mängel zu kompensieren, wurden zahlreiche, erweiterte Verfahren entwickelt. Bis in die 80er-Jahre kamen zur Begrenzung von Oszillationen nahe Unstetigkeitsstel- len zwei Methoden zur Anwendung. In der Ersten sorgte ein künstlicher, sorgfältig und problem- spezifisch ausgelegter Viskositätsterm für eine ausreichende Dämpfung, um Oszillationen als Folge von steilen Gradienten zu unterdrücken. Für die zweite Methode wurden Begrenzungsfunktionen eingeführt, welche unstetige Bereiche mit nicht-oszillierenden Verfahren niedriger Ordnung ap- proximierten mit dem Ziel, durch die Wahl der lokalen Approximationsmethode entsprechende Oszillationen zu unterbinden. Für Details zu diesen Methoden sei auf LeVeque [84] verweisen.
L R
1 2 . . . 𝒾 . . . 𝑁𝑥−1 𝑁𝑥
𝒾+12
Polynom:𝑄(𝑥)mit𝑥 ∈ S
𝑃(1)(𝑥)mit𝑥 ∈ S1
𝑃(2)(𝑥)mit𝑥 ∈ S2
𝑃(3)(𝑥) mit𝑥 ∈ S3
Flussrichtung
𝐿
Geister-Knoten Geister-Knoten
Abbildung 2.3: Rechengitter und -stern für Interpolation 5. Ordnung im Ortsgitterpunkt 𝑥𝒾+12 (Kreuz). blau: Zusammengesetzter Rechenstern, 5. Ordnung. rot: Teil-Rechensterne 𝑃(1)(𝑥), 𝑃(2)(𝑥), 𝑃(3)(𝑥), mit Bias in Flussrichtung, jew. 3. Ordnung. Gestrichelt: sog. Geister- knoten außerhalb der Rechendomäne. Grau: (theoretische) Zellgrenzen als Äquivalent zur FV- Approximation. Bei umgekehrter Flussrichtung sind die Rechensterne entsprechend um𝑥𝒾+12 zu spiegeln.
WENO-Schemata entwickelten sich aus den sog. ENO Verfahren (engl. Akronym: essentially non-oscillatory schemes), welche – durch Harten et al. [55] begründet – erstmalig eine von Git- tergrößen unabhängige, nicht oszillierende Interpolation bei gleichbleibend hoher Approximati- onsordnung ermöglichten. ENO- und WENO-Verfahren nutzen auf die Glattheit der numerischen Lösung angepasste, adaptive Rechensterne, um eine konstant hohe Ordnung und Oszillationsfrei- heit (besonders) in Gegenwart von Unstetigkeiten sicherzustellen.
Die Werte an den Zellgrenzen der Zelle𝒾, z. B. also in𝑥𝒾+12 (Abbildung 2.3,rotesKreuz) sollen aus den bekannten Knotenwerten approximiert werden. Folgende Ausführungen beziehen sich zu- nächst auf den Fall, dass ein positiver Fluss (in Abbildung 2.3 vonlinksnachrechts) vorliegt. Für die Adaption auf einen negativen Nettoflusses sind nachfolgende Ausführungen um den Knoten- wert𝑥𝒾+1
2 zu spiegeln. Zur Approximation des Knotenwertes in𝑥𝒾+1
2 können einfache Polynome𝑄
2engl. Akronym:Weighted Essentially Non-Oscillatory scheme
konstruiert werden, welche𝑟+1 Rechenknoten links bzw.𝑟Rechenknoten rechts der Zellengrenze 𝑥𝒾+12 einschließen (blauin Abbildung 2.3). Der Schätzfehler des gesuchten Knotenwertes nimmt in diesem Fall mit𝑘-ter Ordnung bei Reduktion der GitterweiteΔ𝑥ab:
𝜙ˆ𝒾+1
2 𝑄(𝜙𝒾−𝑟, . . . , 𝜙𝒾+𝑟), 𝒾 =1,2, . . . , 𝑁𝑥 , (2.16) 𝜙ˆ𝒾+1
2 −𝜙𝒾+1
2 =O(Δ𝑥𝑘), 𝑘=2𝑟+1. (2.17)
Die Gesamtspanne der bei der Berechnung abgedeckten Gitterpunkte wird alsgroßerRechenstern bezeichnet
S=
𝑥𝒾−𝑟, . . . , 𝑥𝒾, 𝑥𝒾+1
2, 𝑥𝒾+1, . . . , 𝑥𝒾+𝑟
. (2.18)
Zusätzlich können innerhalb der Spannbreite des großen Sterns(𝑘+1)/2 weitere Polynome,𝑃(𝑠)(𝑥) mit𝑠 =0, . . . ,(𝑘+1)/2−1, konstruiert werden, welche𝑟 +1 Knotenwerte innerhalb vonkleinen Rechensternen S𝑘 ⊆ S (rot in Abbildung 2.3) um die Zielposition 𝑥𝒾+1
2 einschließen und eine verminderte Ordnung von𝑘+1/2=𝑟+1 aufweisen
𝜙ˆ𝒾+(𝑠)1
2 𝑃(𝑠)(𝜙𝒾+𝑠−𝑟, . . . , 𝜙𝒾+𝑠), 𝒾=1,2, . . . , 𝑁𝑥 , und 𝑠 =0, . . . ,(𝑘+1)/2−1 (2.19) 𝜙ˆ𝒾+(𝑠)1
2 −𝜙𝒾+12 =O(Δ𝑥𝑘+12 ) .
Die Kernidee des ENO-Verfahrens besteht darin, zur Interpolation in𝒾+ 12eine der möglichen po- lynomialen Approximationen ˆ𝜙𝒾+(𝑠)1
2
zu wählen, abhängig davon, wie lokal glatt die Knotenwerte sind bzw. ob im großen Rechenstern steile Gradienten vorhanden sind. Der Algorithmus bevor- zugt damit RechensterneS𝑘, welche diese Diskontinuitäten nicht beinhalten. Dies sorgt global für eine Approximationsgüte mindestens(𝑘+1)/2-ter Ordnung, während durch die Wahl des kleinen Rechensterns um den Zielknoten Oszillationen vermieden werden. Dem WENO-Verfahren, als des- sen Erweiterung, liegt zusätzlich die Beobachtung zugrunde, dass sich die Approximation höchster Ordnung ˆ𝜙𝒾+12 als konvexe Linearkombination
𝜙ˆ𝒾+12 =
(𝑘+1)/2−1 𝑠=0
𝑤𝑠𝜙ˆ𝒾+(𝑠)1 2
; 𝑤𝑠 ≥0;
(𝑘+1)/2−1 𝑠=0
𝑤𝑠 =1 (2.20)
der Kleinapproximationen ˆ𝜙𝒾+(𝑠)1 2
darstellen lässt. Zudem gilt
𝑤𝑠 = 𝑤˜𝑠 (𝑘+1)/2−1
𝑠=0 𝑤˜𝑠 mit 𝑤˜𝑠 = 𝑑𝑠
(𝜖+𝛽𝑠)2,
(𝑘+1)/2−1 𝑠=0
𝑑𝑠 =1. (2.21) 𝜖wird zur Vermeidung der Nennernull typischerweise zu𝜖 =10−6gesetzt und hat keinen Einfluss auf die Genauigkeit. Die Gewichte𝑤𝑠 konvergieren bei einer glatten Funktion gegen die Linear- gewichte𝑑𝑠, welche nicht von den Funktionswerten, sondern von dem gewählten Interpolations- polynomen𝑃(𝑠) abhängen. Im Fall einer Diskontinuität konvergieren die Gewichte𝑤𝑠gegen 0 .
2.2 Lösung partieller, eindimensionaler Differentialgleichungen
Der lokale Glattheitsindikator𝛽𝑠 errechnet sich über 𝛽𝑠 =𝑘
𝑙=1
Δ𝑥2𝑙−1
∫ 𝑥
𝒾+1 2 𝑥𝒾−1
2
𝑑𝑙
𝑑𝑥𝑙𝑃(𝑠)(𝑥) 2
𝑑𝑥 (2.22)
als die Summe aller skalierten𝐿2-Normen aller𝑙 Ableitungen des Interpolationspolynoms𝑃(𝑠)(𝑥) im Intervall[𝑥𝒾−1
2, 𝑥𝒾+1
2]. Diese können analytisch aus den Interpolationspolynomen (Gleichung 2.19) bestimmt werden [135].
Die in dieser Arbeit genutzten Rechensterne sind dritter bzw. insgesamt fünfter Ordnung (wie in Abb. 2.3 dargestellt), wobei das resultierende Schema als „WENO-35” abgekürzt wird. Die ver- wendeten Lagrange-Polynome 𝑄 bzw. 𝑃(0), . . . , 𝑃(2), die expliziten Formeln für 𝛽+/−0 , . . . , 𝛽2+/−
und ˜𝑤0,0, . . . ,𝑤˜2,0 und die Parameter 𝑤0, . . . ,𝑤2 bzw. 𝑑0, . . . ,𝑑2 dieser Polynome sind in An- hang A2, S. 151 zu finden.
Berechnung des Transports von Erhaltungsgrößen
Da WENO-Methoden typischerweise zur Lösung partieller Differentialgleichungen, insbesondere physikalischer, hyberbolischer Kontinuitätsgleichungen in differenzieller Form eingesetzt werden, soll an dieser Stelle die Anwendung auf ein Transportsystem folgen. Um mit gegebener Methode eindimensionale Erhaltungsgleichungen aus Zellwerten approximieren zu können, muss eine (im Wortsinne) konservative, d. h. erhaltende Schätzung der Ableitungen einer beliebigen Größe aus diskreten Gitterpunkten berechnet werden. Für die Herleitung wird an dieser Stelle die differen- tielle Form der Erhaltungsgleichung 2.1 mit𝑆𝜙 = 0 genutzt. Eine an den Rechenknoten bekannte Größe𝜙(𝑥)werde damit durch
𝜕𝜙(𝑥, 𝑡)
𝜕𝑡 + 𝜕𝑓(𝜙(𝑥, 𝑡))
𝜕𝑥 =0 (2.23)
über die gegebene Flussfunktion𝑓 transportiert. Hierfür solle die Ortsableitung𝜕𝑓(𝜙(𝑥, 𝑡))/𝜕𝑥aus der Differenz der sog. numerischen Flüsse ˆ𝑓𝒾±1
2 𝑓(𝜙ˆ𝒾±1
2) mit einer Genauigkeit k-ter Ordnung approximiert werden
1 Δ𝑥
𝑓ˆ𝒾+1
2 −𝑓ˆ𝒾−1
2
= 𝜕𝑓(𝜙𝒾, 𝑥𝒾)
𝜕𝑥 + O(Δ𝑥𝑘), 𝒾 =1,2, . . . , 𝑁𝑥, (2.24) wobei die numerischen Flüsse konsistent mit den physikalischen sein müssen. Mit Gl. 2.3 folgt daraus:
𝑑𝜙𝒾(𝑡) 𝑑𝑡 ≈ − 1
Δ𝑥
𝑓ˆ𝒾+12 −𝑓ˆ𝒾−12
. (2.25)
Die Berechnung der numerischen Flüsse ˆ𝑓𝒾±1
2 kann über das WENO-Verfahren gemäß Gl. 2.19, 2.20 und 2.21 aus den umgebenden Knotenwerten
𝑓ˆ𝒾±12 =𝑓
(𝑘+1)/2−1
𝑠=0
𝑤𝑠𝜙ˆ𝒾±(𝑠)1 2
=𝑓
(𝑘+1)/2−1
𝑠=0
𝑤𝑠𝑃(𝑠)(𝜙𝒾+𝑠−𝑟, . . . , 𝜙𝒾+𝑠)
, 𝑖 =1,2, . . . , 𝑁𝑥 , (2.26)
oder – im Falle einer linearen Flussfunktion 𝑓(𝜙) – direkt aus den an den Knotenwerten errech- neten Flüssen𝑓𝒾+12 =𝑓(𝜙𝒾+12)
𝑓ˆ𝒾±12 =
(𝑘+1)/2−1 𝑠=0
𝑤𝑠𝑓ˆ𝒾±(𝑠)1
2 =
(𝑘+1)/2−1 𝑠=0
𝑤𝑠𝑃(𝑠)(𝑓𝒾+𝑠−𝑟, . . . , 𝑓𝒾+𝑠), 𝒾=1,2, . . . , 𝑁𝑥 (2.27) mit den nach Gleichung 2.21 errechneten Gewichten erfolgen.
Stabilität kann für advektionsdominante Probleme nur für Rechensterne garantiert werden, die aus überwiegend Knotenwerten stromauf, d. h. mit einem Bias entgegen der Strömungsrichtung, konstruiert werden [41]. Dies wurde bereits bei der Konstruktion der WENO-Polynome berück- sichtigt, siehe Abbildung 2.3. Die Rekonstruktion eines Knotenwertes in𝑥𝒾+1
2 auf Grundlage von ausschließlich Knotenwerten stromab dieses Knotens (im Falle eines reinen Advektionstransports) wäre – auch nach allgemeinem Verständnis – in hohen Maße unphysikalisch.
Um mit dem WENO-Verfahren auch bei einer Umkehr der Strömungsrichtung stabil die Werte an den Zellgrenzen schätzen zu können, muss die Schätzung des numerischen Flusses ˆ𝑓𝒾+1
2 mit einem Rechenstern erfolgen, der exakt um den Punkt𝑥𝒾+12 gespiegelt ist. Dies wird typischerweise durch einen Teilungsansatz (engl.: flux splitting) realisiert, der ˆ𝑓𝒾+12 in links- bzw. rechtsseitig wirkende Flussanteile (Superindex−bzw.+) aufteilt
𝑓ˆ𝒾±1
2(𝜙) =𝑓ˆ𝒾±+1
2(𝜙) +𝑓ˆ𝒾±−1
2(𝜙), mit 𝜕𝑓ˆ+
𝜕𝜙 ≥0 und 𝜕𝑓ˆ−
𝜕𝜙 ≤ 0. (2.28)
Für den bekanntesten, in dieser Arbeit verwendeten Ansatz nach Lax - Friedrich [82] gilt 𝑓ˆ𝒾±±1
2(𝜙) = 1 2
𝑓ˆ𝒾±1
2(𝜙) ±𝛼𝜙
mit 𝛼=max
𝜙 |𝑓(𝜙)| =konst. (2.29) Die Abkürzung des verwendeten Teilungsansatzes wird üblicherweise als Suffix (-LF) an den Na- men des Diskretisierungsverfahrens angehängt (hier: WENO-LF).
2.2.5 Zeitdiskretisierung: Expliziter RK-45
Die Stabilität von WENO-Methoden in Kombination mit expliziten Integratoren wurde numerisch untersucht [93]. Hierbei erwiesen sich Runge-Kutta-Methoden höherer Ordnung als am besten geeignet für den Nutzen in Kombination mit WENO-Schemata. Die zeitliche Integration erfolgt in dieser Arbeit daher mit einem expliziten Verfahren vierter Ordnung nach Fehlberg [40], welches als Runge-Kutta-Fehlberg-Verfahren bekannt wurde. Bei diesem Verfahren handelt es sich um ein Einschritt-Verfahren mit adaptiver Schrittweiten-Steuerung, das in jedem Integrationsschritt zwei verschiedene Approximationen der Lösung berechnet und vergleicht. Die Integrationsstufen des Verfahrens werden hierbei rekursiv berechnet.
