TU Darmstadt / ULB / TUprints

Mathematical modeling and Volume-of-Fluid based simulation of dynamic wetting

Fricke, Mathis (2021):
Mathematical modeling and Volume-of-Fluid based simulation of dynamic wetting. (Publisher's Version)
Darmstadt, Technische Universität Darmstadt,
DOI: 10.12921/tuprints-00014274,
[Ph.D. Thesis]

Available under CC-BY-NC-SA 4.0 International - Creative Commons, Attribution Non-commercial, Share-alike.

Download (22MB) | Preview
Item Type: Ph.D. Thesis
Status: Publisher's Version
Title: Mathematical modeling and Volume-of-Fluid based simulation of dynamic wetting
Language: English

Dynamic wetting phenomena are omnipresent in nature and technology. The legs of the water-strider make use of a sophisticated hierarchical surface structure to achieve superhydrophobicity and to allow the insect to stand and run easily on a water surface. The ability to understand and control dynamic wetting processes is crucial for a variety of industrial and technical processes such as inkjet- or bioprinting or mass transport in microfluidic devices. On the other hand, the moving contact line problem, even in a largely simplified setting, still poses considerable challenges regarding fundamental mathematical modeling as well as numerical methods. The present work addresses both the fundamental modeling and the development of numerical methods based on the geometrical Volume-of-Fluid (VOF) method.

The spreading dynamics of droplets on a bare silicon wafer and on a silicon wafer coated with a polymer brush is studied in cooperation with experimentalists within the Collaborative Research Center (CRC) 1194. An ordinary differential equation describing the spreading dynamics of spherical drops is derived and compared with experimental results. The model is a generalization of a classical model for perfectly wetting drops to the case of partial wetting. Besides these simplified modeling approaches, the main focus of the present work lies on the continuum mechanical description of dynamic wetting.

The moving contact line singularity in the classical hydrodynamic description based on the no-slip boundary condition motivated a lot of research in the past 50 years, aiming at a physically sound model. It has been shown that the Navier slip condition combined with a fixed contact angle leads to a so-called "weak singularity" and it was suspected by Ren and E that the solution may become completely regular for Navier slip combined with a dynamic contact angle (Phys. Fluids, 2007). The central mathematical tool developed in the present work allows to prove that the latter conjecture is false (as long as the slip length is finite). The basic idea is to study the kinematics of wetting in the sharp-interface sharp-contact-line setting independently from the specific continuum mechanical model. An evolution equation for the dynamic contact angle is derived and proven rigorously, assuming that a sufficiently regular velocity field is given on the moving hypersurface with boundary. Thanks to this very general setting, the result is applicable to a large class of continuum mechanical models including the mechanisms of mass transfer between the phases or mass transfer to a surface phase like in the Interface Formation Model. The kinematic result is applied to regular solutions of the "standard model" of dynamic wetting based on the Navier slip condition. It is shown that the system cannot relax to the equilibrium state with a regular solution. Hence, it is concluded that physically sound solutions in the standard model cannot be regular. Moreover, regular solutions to generalizations of the standard model are studied. In particular, it is shown that surface tension gradients at the contact line may give rise to regular solutions.

Furthermore, the compatibility of the boundary conditions at the contact line is studied for the standard model and an adaptation of the Interface Formation Model proposed recently by Lukyanov and Pryer (Langmuir, 2017). It is shown that, depending on the model parameters, the boundary conditions in the model by Lukyanov and Pryer may be compatible at the contact line. In this case, one can even compute explicit expressions for the curvature and the pressure at the moving contact line.

The second part of the present thesis is devoted to numerical methods for dynamic wetting. In order to make the kinematic results easily accessible, an open-source demonstrator code based on a level set representation of the interface is developed and published in an open repository. The current state-of-the-art methods for dynamic wetting based on the geometrical Volume-of-Fluid approach are briefly reviewed. In particular, it is shown that the method to enforce the dynamic contact angle proposed by Afkhami and Bussmann (Int. J. Numer. Methods Fluids, 2008) delivers inconsistent values for the curvature at the contact line. Motivated by the fundamental results on the kinematics of moving contact lines, novel interface reconstruction methods are developed and implemented that allow to reconstruct the free surface close to the domain boundary. In particular, the Boundary ELVIRA method delivers an accurate numerical transport of the contact angle that is consistent with the fundamental kinematics. The latter method greatly improves the accuracy of the VOF method in the presence of contact lines.

Moreover, the numerical approximation of the mean curvature based on the height function technique is studied thoroughly. A rigorous error analysis for the two-dimensional height function method in the presence of data errors is given. In particular, the discrete error amplification is estimated and studied in detail. The latter type of error is rarely discussed in the scientific literature on the topic. But in fact, the impact of the discrete error amplification on the total error can be significant, in particular when disturbances due to transport errors are present in the volume fraction data. The kinematic evolution equation for the mean curvature, which is derived in the first part of this work, serves as a reference solution to validate the numerical transport of the curvature at the contact line. As can be expected from the height function error analysis, the transport error for the curvature is found to be first-order divergent for the Boundary Youngs method and constant for the Boundary ELVIRA method. The latter results clearly show the need for higher-order interface advection methods.

The third part of this work closely investigates two particular wetting flow configurations, namely the capillary rise and the breakup of a liquid bridge on a chemically structured surface. A novel numerical benchmark for wetting flows based on the capillary rise is established with four numerical methods developed within the CRC 1194 at TU Darmstadt. Moreover, a novel adaptation of the Navier slip condition called "staggered slip" is introduced. The goal of the staggered slip condition is to reduce the "numerical slip" inherent to the method. This is achieved by defining the slip length with respect to a virtual boundary that is located in between the physical boundary and location of the face-centered velocity used to transport the volume fraction field. As a result, the discrete viscous dissipation is increased compared to the standard Navier slip implementation. It is shown that the convergence for the capillary rise can be significantly improved if the slip length is not yet resolved. On the other hand, the order of convergence is reduced compared to the standard implementation for a single-phase channel flow example.

The wetting of structured surfaces is studied in joint work with experimentalists in the CRC 1194. The goal is to quantitatively describe the breakup dynamics of a wetting capillary bridge on a structured surface. A major problem for the interpretation of both the experimental and the numerical data arises from the uncertainty in the precise time of pinch-off of the capillary bridge. In order to solve this problem, we introduce a systematic way to analyze the data without the need to determine the pinch-off time. The basic idea, which has been applied before by Li Sprittles (J. Fluid Mech., 2016), is to plot the speed of the breakup process (i.e. the time derivative of the minimum diameter) as a function of the minimum diameter itself. This procedure is well-defined since the minimum diameter is strictly decreasing with time. Indeed, we show that the transformation that maps from the standard representation to the phase space representation is invertible up to a shift in absolute time. With this technique, we are able to study the breakup process in great detail in both the three-dimensional VOF simulation and the experiment. In general, a good agreement is found between experiment and simulation both qualitatively and quantitatively in terms of the time evolution of the minimum diameter. The numerical simulations allowed to identify different regimes in the breakup dynamics that were also found in the experimental data. Remarkably, dynamic surface tension may play a significant role in the breakup dynamics. The agreement between simulation and experiment close to the breakup can be improved by increasing the value of the (constant) surface tension to 90 mN/m. The latter value has been proposed by Hauner et al. (J. Phys. Chem. Lett., 2017) for a freshly created water surface. Moreover, the local rate of interface generation is found to be quite high. However, a fully quantitative assessment of this phenomenon can only be achieved in future work. Finally, we revisit the capillary rise problem in the case of structured surfaces. Interestingly, the surface pattern can be used to construct an energy functional with two stable configurations. The resulting dynamic rise behavior is quite complex with a bifurcation between the two stable configurations at a critical initial rise height.

Alternative Abstract:
Alternative AbstractLanguage
Dynamische Benetzungsphänomene sind in Natur und Technik allgegenwärtig. Die Beine eines Wasserläufers nutzen eine ausgeklügelte hierarchische Oberflächenstruktur, um Superhydrophobie zu erreichen und das Insekt auf einer Wasseroberfläche leicht stehen und laufen zu lassen. Die Fähigkeit, dynamische Benetzungsprozesse zu verstehen und zu steuern, ist entscheidend für eine Vielzahl industrieller und technischer Prozesse wie Bioprinting und Tintenstrahldruck oder Massentransport in Mikrofluidikgeräten. Andererseits birgt das Problem der beweglichen Kontaktlinie selbst in einer stark vereinfachten Formulierung immer noch erhebliche Herausforderungen hinsichtlich der fundamentalen mathematischen Modellierung sowie der numerischen Methoden. Die vorliegende Arbeit befasst sich sowohl mit der grundlegenden Modellierung als auch mit der Entwicklung numerischer Methoden auf der Grundlage der geometrischen VOF (Volume-of-Fluid) Methode. Die Ausbreitungsdynamik von Tropfen auf einem unbehandelten Siliziumwafer und auf einem mit einer Polymerbürste beschichteten Siliziumwafer wird in Zusammenarbeit mit Experimentatoren im Sonderforschungsbreich (SFB) 1194 untersucht. Eine gewöhnliche Differenzialgleichung, die die Ausbreitungsdynamik sphärischer Tropfen beschreibt, wird abgeleitet und mit experimentellen Ergebnissen verglichen. Das Modell ist eine Verallgemeinerung eines klassischen Modells für perfekt benetzende Tropfen auf den Fall der teilweisen Benetzung. Neben diesen vereinfachten Modellierungsansätzen liegt der Schwerpunkt der vorliegenden Arbeit auf der kontinuumsmechanischen Beschreibung der dynamischen Benetzung. Die Singularität an einer sich bewegenden Kontaktlinie in der klassischen hydrodynamischen Beschreibung, basierend auf der Haftbedingung an festen Rändern, hat in den letzten 50 Jahren eine Reihe von Forschungsaktivitäten angeregt, die auf ein physikalisch fundiertes Modell abzielen. Es wurde gezeigt, dass die Navier-Schlupfbedingung in Kombination mit einem festen Kontaktwinkel zu einer sogenannten schwachen Singularität führt, und es wurde vermutet, dass die Lösung für Navier-Schlupf in Kombination mit einem dynamischen Kontaktwinkel vollständig regularisiert wird (Phys. Fluids, 2007). Das in der vorliegenden Arbeit entwickelte zentrale mathematische Werkzeug ermöglicht es zu beweisen, dass die letztere Vermutung falsch ist (solange die Schlupflänge endlich ist). Die Grundidee besteht darin, die Kinematik der Benetzung in der Klasse der sharp-interface sharp-contact-line Modelle unabhängig vom spezifischen kontinuumsmechanischen Modell zu studieren. Eine Evolutionsgleichung für den dynamischen Kontaktwinkel wird abgeleitet und rigoros unter der Annahme bewiesen, dass ein ausreichend reguläres Geschwindigkeitsfeld auf der berandeten Mannigfaltigkeit gegeben ist. Dank dieser sehr allgemeinen Formulierung kann das Resultat auf eine große Klasse von kontinuumsmechanischen Modellen angewendet werden, einschließlich der Mechanismen des Stofftransfers zwischen zwei Phasen oder des Stofftransfers auf eine Oberflächenphase wie im Interface Formation Modell. Das kinematische Resultat wird auf reguläre Lösungen des Standardmodells der dynamischen Benetzung, welches auf der Navier-Schlupfbedingung basiert, angewendet. Es wird gezeigt, dass eine Relaxation des Systems in den Gleichgewichtszustand mit einer regulären Lösung nicht möglich ist. Daher lässt sich die Schlussfolgerung ableiten, dass physikalische Lösungen im Standardmodell nicht regulär sein können. Darüber hinaus werden reguläre Lösungen für Verallgemeinerungen des Standardmodells untersucht. Insbesondere wird gezeigt, dass Oberflächenspannungsgradienten an der Kontaktlinie zu regulären Lösungen führen können. Darüber hinaus wird die Kompatibilität der Randbedingungen an der Kontaktlinie für das Standardmodell und eine kürzlich von Lukyanov und Pryer vorgeschlagene Adaption des Interface Formation Modells (Langmuir, 2017) untersucht. Es wird gezeigt, dass abhängig von den Modellparametern die Randbedingungen im Modell von Lukyanov und Pryer an der Kontaktlinie kompatibel sein können. In diesem Fall lassen sich sogar explizite Ausdrücke für die Krümmung und den Druck an der sich bewegenden Kontaktlinie berechnen. Der zweite Teil der vorliegenden Arbeit befasst sich mit numerischen Methoden zur dynamischen Benetzung. Um die kinematischen Ergebnisse leicht zugänglich zu machen, wird ein, auf der Level Set Methode basierender, Open-Source Demonstratorcode entwickelt und frei zugänglich gemacht. Der aktuelle Stand der Forschung bezüglich geometrischer VOF Verfahren zur Beschreibung von dynamischen Benetzungsprozessen wird kurz zusammengefasst. Insbesondere wird gezeigt, dass das von Afkhami und Bussmann vorgeschlagene Verfahren zur Einstellung des dynamischen Kontaktwinkels (Int. J. Numer. Methods Fluids, 2008) inkonsistente Werte für die Krümmung an der Kontaktlinie liefert. Motiviert durch die grundlegenden Ergebnisse zur Kinematik dynamischer Kontaktlinien werden neuartige Methoden zur Rekonstruktion der freien Grenzflächen nahe der Domänengrenze entwickelt und implementiert. Insbesondere ermöglicht die Boundary ELVIRA-Methode einen genauen numerischen Transport des Kontaktwinkels, der mit der fundamentalen Kinematik übereinstimmt. Das letztere Verfahren verbessert die Genauigkeit des VOF-Verfahrens bei Vorhandensein von dynamischen Kontaktlinien erheblich. Darüber hinaus wird die numerische Approximation der mittleren Krümmung grundlegend untersucht. Eine rigorose Fehleranalyse für die zweidimensionale Höhenfunktionsmethode bei Vorhandensein von Datenfehlern wird durchgeführt. Insbesondere wird die diskrete Fehlerverstärkung abgeschätzt und detailliert untersucht. Die letztere Art von Fehler wird in der wissenschaftlichen Literatur zu diesem Thema selten diskutiert. Tatsächlich kann der Einfluss der diskreten Fehlerverstärkung auf den Gesamtfehler jedoch erheblich sein, insbesondere wenn Störungen aufgrund von Transportfehlern in den Volumenfraktionsdaten vorhanden sind. Die kinematische Evolutionsgleichung für die mittlere Krümmung, die im ersten Teil dieser Arbeit abgeleitet wird, dient als Referenzlösung zur Validierung des numerischen Transports der Krümmung an der Kontaktlinie. Wie aus der Fehleranalyse zu erwarten ist, wird festgestellt, dass der Transportfehler für die Krümmung für die Boundary Youngs-Methode divergent erster Ordnung und für die Boundary ELVIRA-Methode konstant ist. Die Ergebnisse zeigen deutlich den Bedarf an Grenzflächen-Advektionsmethoden höherer Ordnung. Der dritte Teil dieser Arbeit untersucht zwei spezifische Beispiele, nämlich den Anstieg einer Flüssigkeitssäule in einer Kapillare (Capillary Rise) und das Aufbrechen einer Flüssigkeitsbrücke auf einer chemisch strukturierten Oberfläche. Ein vergleichender numerischer Benchmark für dynamische Benetzungsprozesse wird am Beispiel des Capillary Rise entwickelt. Dazu kommen vier völlig verschiedene numerische Methoden zum Einsatz, die im SFB 1194 entwickelt werden. Das Resultat ist eine fundierte Datenbasis zur Validierung künftiger numerischer Verfahren. Darüber hinaus wird eine neuartige Diskretisierung der Navier-Schlupfbedingung eingeführt, die als versetzter Schlupf oder staggered slip bezeichnet wird. Das Ziel der versetzten Schlupfbedingung besteht darin, den der Methode inhärenten numerischen Schlupf zu verringern. Dies wird erreicht, indem die Schlupflänge in Bezug auf einen virtuellen Rand definiert wird, der sich zwischen dem physikalischen Rand und dem Ort der flächenzentrierten Geschwindigkeit befindet, welche zum Transport des Volumenfraktionsfeldes verwendet wird. Infolgedessen ist die diskrete viskose Dissipation im Vergleich zur Standard Navier-Schlupf-Implementierung erhöht. Es wird gezeigt, dass die Konvergenz für den Kapillaranstieg signifikant verbessert werden kann, wenn die Schlupflänge noch nicht aufgelöst ist. Andererseits ist die Konvergenzordnung für eine einphasige Kanalströmung im Vergleich zur Standardimplementierung reduziert. Die Benetzung strukturierter Oberflächen wird in Zusammenarbeit mit Experimentatoren im SFB 1194 untersucht. Ziel ist es, die Aufbruchdynamik einer benetzenden Kapillarbrücke auf einer strukturierten Oberfläche quantitativ zu beschreiben. Ein Hauptproblem für die Interpretation sowohl der experimentellen als auch der numerischen Daten ergibt sich aus der Unsicherheit bezüglich des genauen Zeitpunkts des Aufbruchs der Kapillarbrücke. Um dieses Problem zu lösen, führen wir eine systematische Methode zur Analyse der Daten ein, die keine Bestimmung der Aufbruchzeit erfordert. Die Grundidee, die zuvor schon von Li Sprittles (J. Fluid Mech., 2016) angewendet wurde, besteht darin, die Geschwindigkeit des Prozesses (d.h. die zeitliche Ableitung des minimalen Durchmessers) als Funktion des minimalen Durchmessers selbst zu untersuchen. Dieses Verfahren ist wohldefiniert, da der Mindestdurchmesser monoton abnimmt. In der Tat wird gezeigt, dass die Transformation, die von der Standarddarstellung auf die Phasenraumdarstellung abbildet, invertierbar ist, bis auf eine Verschiebung in der absoluten Zeit. Mit dieser Methode können wir den Aufbruchprozess sowohl in der dreidimensionalen VOF-Simulation als auch im Experiment sehr detailliert untersuchen. Im Allgemeinen wird eine gute Übereinstimmung zwischen Experiment und Simulation sowohl qualitativ als auch quantitativ in Bezug auf die zeitliche Entwicklung des minimalen Durchmessers erreicht. Die numerischen Simulationen ermöglichen es, verschiedene Regime in der Aufbruchdynamik zu identifizieren, die auch in den experimentellen Daten gefunden wurden. Bemerkenswerterweise könnten dynamische Oberflächenspannungseffekte eine wichtige Rolle bei der Aufbruchdynamik spielen. Die Übereinstimmung zwischen Simulation und Experiment kurz vor dem Aufbruch kann verbessert werden, indem der Wert der (konstanten) Oberflächenspannung auf 90 mN/m erhöht wird. Der letztere Wert wurde von Hauner et al. (J. Phys. Chem. Lett., 2017) für eine frisch erzeugte Wasseroberfläche vorgeschlagen. Darüber hinaus wird durch Auswertung der numerischen Daten festgestellt, dass die lokale Grenzflächenerzeugungsrate sehr hoch ist. Eine vollständig quantitative Bewertung dieses Phänomens steht jedoch noch aus. Schließlich betrachten wir noch einmal das Capillary Rise Problem, diesmal für eine strukturierte Oberfläche. Interessanterweise kann die Oberflächenstrukturierung dazu verwendet werden, ein Energiefunktional mit zwei stabilen Konfigurationen zu konstruieren. Das resultierende dynamische Aufstiegsverhalten ist sehr komplex und zeigt eine Bifurkation zwischen den beiden stabilen Gleichgewichtslagen bei einer kritischen initialen Steighöhe.German
Place of Publication: Darmstadt
Collation: xxi, 220 Seiten
Classification DDC: 500 Naturwissenschaften und Mathematik > 510 Mathematik
Divisions: 04 Department of Mathematics > Analysis > Mathematical Modeling and Analysis
Date Deposited: 13 Jan 2021 09:52
Last Modified: 13 Jan 2021 10:11
DOI: 10.12921/tuprints-00014274
URN: urn:nbn:de:tuda-tuprints-142745
Referees: Bothe, Prof. Dr. Dieter and Ulbrich, Prof. Dr. Stefan and Zaleski, Prof. Stéphane
Refereed: 27 October 2020
URI: https://tuprints.ulb.tu-darmstadt.de/id/eprint/14274
Actions (login required)
View Item View Item