
MATHEMATICAL MODELING AND

VOLUME-OF-FLUID BASED SIMULATION

OF DYNAMIC WETTING

Vom Fachbereich Mathematik
der Technischen Universität Darmstadt

zur Erlangung des Grades eines
Doktors der Naturwissenschaften (Dr. rer. nat.)

genehmigte Dissertation von

M. SC. MATHIS FRICKE

aus Alsfeld

Referent: Prof. Dr. Dieter BOTHE, TU Darmstadt

Korreferent: Prof. Dr. Stefan ULBRICH, TU Darmstadt

Korreferent: Prof. Stéphane ZALESKI, Sorbonne Université, Paris

Darmstadt 2020
D 17


Fricke, Mathis: Mathematical modeling and Volume-of-Fluid based
simulation of dynamic wetting,
Darmstadt, Technische Universität Darmstadt

Tag der mündlichen Prüfung: 27.10.2020
Jahr der Veröffentlichung der Dissertation of TUprints: 2021

URN: urn:nbn:de:tuda-tuprints-142745

This work is published under the Creative Commons License:
Attribution – NonCommercial – ShareAlike (CC BY-NC-SA 4.0)
https://creativecommons.org/licenses/by-nc-sa/4.0/

urn:nbn:de:tuda-tuprints-142745
https://creativecommons.org/licenses/by-nc-sa/4.0/


Für meine Familie

MAGDALENE, ANTON EMIL UND JAKOB ELIAS.


Acknowledgments
First of all, I would like to express my deep gratitude to my supervisor Prof. Dieter Bothe for supporting me as a
person and as a young researcher. He gave me the opportunity to join a unique research environment in the Math-
ematical Modeling and Analysis (MMA) group at TU Darmstadt. I am very grateful for his continuous support of
my work and lots of deep scientific discussions.

I very much appreciated being part of the MMA team over the last years. In particular, I thank my colleagues
Johannes Kromer and Dirk Gründing for proof-reading parts of the manuscript of the present thesis and a lot of
enlightening discussions. Furthermore, I very much enjoyed working closely with Maximilian Hartmann to match
numerical simulations and experiments in dynamic wetting. The numerous discussions with Dirk Gründing in our
daily work within CRC 1194 were always useful to make the project a success. Moreover, I am very lucky to have
Mattias Köhne (HHU Düsseldorf) as a close collaborator in my project. He greatly supported me in my first steps
as a researcher and I am looking forward to continue our joint work in the future.

Referees for the Thesis: I am very grateful to Prof. Dieter Bothe, Prof. Stefan Ulbrich, and Prof. Stéphane
Zaleski for acting as referees for the present thesis.

Close collaboration within CRC 1194: I greatly benefited from being part of a great team of researches within
CRC 1194. In particular, I would like to thank Maximilian Hartmann, Dirk Gründing, Tomislav Marić, Mar-
tin Smuda, Thomas Antritter, Daniel Rettenmaier, Beatrice Fickel, Elisabeth Diehl, Florian Kummer, Holger
Marschall, Steffen Hardt, Ilia Roisman, Tatiana Gambaryan-Roisman, Tobias Pfeiffer, Edder J. Garcı́a, Nico van
der Vegt, Günter Auernhammer and Peyman Rostami for the close collaboration within CRC 1194.

Organization of CRC 1194: I would like to express my gratefulness to Prof. Peter Stephan, Prof. Dieter Bothe,
Dr. Benjamin Lambie, Monika Medina España and Sebastian Keuth for their daily efforts in running the collab-
orative research center and thereby creating an excellent environment for young researches like myself to follow
their passion for science.

Supervised students within CRC 1194: I am happy to thank Fabian Gabel, João Victor Machado Lima, Fu
Xing Long and Aleksandar Vučković for their efforts as undergraduate assistants in subproject B01 of CRC 1194.
In particular, Aleksandar Vučković contributed greatly to the project with the implementation of a level set method
in C++. I am grateful for his persistent effort in supporting the project.

Scientific discussions: During the last years I very much appreciated scientific exchange and discussions with
Prof. Shahriar Afkhami, Dr. Terry Blake, Dr. Elmar Bonaccurso, Prof. Longquan Chen, Prof. Joel De Coninck,
Prof. Serafim Kalliadasis, Dr. Dirk Peschka, Prof. Weiqing Ren, Prof. Yulii Shikhmurzaev, Prof. Jacco Snoeijer
and Prof. James Sprittles.

Financial support: I kindly acknowledge the financial support by the German Research Foundation (DFG)
within the Collaborative Research Centre 1194 “Interaction of Transport and Wetting Processes” – Project-ID
265191195 , subproject B01.

HPC infrastructure: Calculations for this research were conducted on the Lichtenberg high performance com-
puter of the TU Darmstadt.


Danksagung
Schließlich möchte ich meiner Familie meine tiefe Dankbarkeit für all die Unterstützung in den letzten Jahren
ausdrücken. Ohne Euch wäre der erfolgreiche Abschluss meiner Promotion niemals möglich gewesen. Besonders
meiner Frau Magdalene und meinen Söhnen Anton Emil und Jakob Elias möchte ich danken für den geduldigen
Verzicht auf viele gemeinsame Stunden in der letzten Zeit und all die Unterstützung. Auch meinen Eltern Susanne
und Michael und meinem Bruder Fabian gilt mein besonderer Dank für Eure Begleitung und Unterstützung vom
Kindesalter über die Schule bis zum Studium und weit darüber hinaus.


Summary
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 experi-
mental 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 con-
dition 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 kine-
matic 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

vii

https://www.doi.org/10.1063/1.2646754
http://dx.doi.org/10.1021/acs.langmuir.7b02409
https://www.doi.org/10.1002/fld.1651
https://www.doi.org/10.1002/fld.1651


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 wet-
ting 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
90mN/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 quan-
titative 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.

viii

https://www.doi.org/10.1017/jfm.2016.276
http://dx.doi.org/10.1021/acs.jpclett.7b00267


Zusammenfassung
Dynamische Benetzungsphänomene sind in Natur und Technik allgegenwärtig. Die Beine eines Wasserläufers nut-
zen eine ausgeklügelte hierarchische Oberflächenstruktur, um Superhydrophobie zu erreichen und das Insekt auf
einer Wasseroberfläche leicht stehen und laufen zu lassen. Die Fähigkeit, dynamische Benetzungsprozesse zu ver-
stehen und zu steuern, ist entscheidend für eine Vielzahl industrieller und technischer Prozesse wie Bioprinting und
Tintenstrahldruck oder Massentransport in Mikrofluidikgeräten. Andererseits birgt das Problem der beweglichen
Kontaktlinie selbst in einer stark vereinfachten Formulierung immer noch erhebliche Herausforderungen hinsicht-
lich 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 Poly-
merbürste beschichteten Siliziumwafer wird in Zusammenarbeit mit Experimentatoren im Sonderforschungsbreich
(SFB) 1194 untersucht. Eine gewöhnliche Differenzialgleichung, die die Ausbreitungsdynamik sphärischer Trop-
fen beschreibt, wird abgeleitet und mit experimentellen Ergebnissen verglichen. Das Modell ist eine Verallgemei-
nerung eines klassischen Modells für perfekt benetzende Tropfen auf den Fall der teilweisen Benetzung. Neben
diesen vereinfachten Modellierungsansätzen liegt der Schwerpunkt der vorliegenden Arbeit auf der kontinuums-
mechanischen Beschreibung der dynamischen Benetzung.

Die Singularität an einer sich bewegenden Kontaktlinie in der klassischen hydrodynamischen Beschreibung, basie-
rend auf der Haftbedingung an festen Rändern, hat in den letzten 50 Jahren eine Reihe von Forschungsaktivitä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 Singularität“ führt, und es
wurde vermutet, dass die Lösung für Navier-Schlupf in Kombination mit einem dynamischen Kontaktwinkel
vollständig regularisiert wird (Ren and E, Phys. Fluids, 2007). Das in der vorliegenden Arbeit entwickelte zen-
trale mathematische Werkzeug ermöglicht es zu beweisen, dass die letztere Vermutung falsch ist (solange die
Schlupflänge endlich ist). Die Grundidee besteht darin, die Kinematik der Benetzung in der Klasse der sharp-
interface sharp-contact-line Modelle unabhängig vom spezifischen kontinuumsmechanischen Modell zu studieren.
Eine Evolutionsgleichung für den dynamischen Kontaktwinkel wird abgeleitet und rigoros unter der Annahme
bewiesen, dass ein ausreichend regulä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 Oberflächenphase wie im Interface Formation Modell. Das kinematische Resultat wird auf
reguläre Lö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
regulären Lösung nicht möglich ist. Daher lässt sich die Schlussfolgerung ableiten, dass physikalische Lösungen
im Standardmodell nicht regulär sein können. Darüber hinaus werden reguläre Lösungen für Verallgemeinerungen
des Standardmodells untersucht. Insbesondere wird gezeigt, dass Oberflächenspannungsgradienten an der Kon-
taktlinie zu regulären Lösungen führen können.
Darüber hinaus wird die Kompatibilität der Randbedingungen an der Kontaktlinie für das Standardmodell und ei-
ne kürzlich von Lukyanov und Pryer vorgeschlagene Adaption des Interface Formation Modells (Langmuir, 2017)
untersucht. Es wird gezeigt, dass abhängig von den Modellparametern die Randbedingungen im Modell von Lu-
kyanov und Pryer an der Kontaktlinie kompatibel sein können. In diesem Fall lassen sich sogar explizite Ausdrücke
für die Krü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 zugänglich zu machen, wird ein, auf der Level Set Methode basierender, Open-
Source Demonstratorcode entwickelt und frei zugänglich gemacht. Der aktuelle Stand der Forschung bezüglich
geometrischer VOF Verfahren zur Beschreibung von dynamischen Benetzungsprozessen wird kurz zusammenge-
fasst. Insbesondere wird gezeigt, dass das von Afkhami und Bussmann vorgeschlagene Verfahren zur Einstellung
des dynamischen Kontaktwinkels (Int. J. Numer. Methods Fluids, 2008) inkonsistente Werte für die Krümmung an
der Kontaktlinie liefert. Motiviert durch die grundlegenden Ergebnisse zur Kinematik dynamischer Kontaktlinien
werden neuartige Methoden zur Rekonstruktion der freien Grenzflächen nahe der Domänengrenze entwickelt und

ix

https://www.doi.org/10.1063/1.2646754
http://dx.doi.org/10.1021/acs.langmuir.7b02409
https://www.doi.org/10.1002/fld.1651


implementiert. Insbesondere ermöglicht die Boundary ELVIRA-Methode einen genauen numerischen Transport
des Kontaktwinkels, der mit der fundamentalen Kinematik übereinstimmt. Das letztere Verfahren verbessert die
Genauigkeit des VOF-Verfahrens bei Vorhandensein von dynamischen Kontaktlinien erheblich. Darüber hinaus
wird die numerische Approximation der mittleren Krümmung grundlegend untersucht. Eine rigorose Fehlerana-
lyse für die zweidimensionale Höhenfunktionsmethode bei Vorhandensein von Datenfehlern wird durchgeführt.
Insbesondere wird die diskrete Fehlerverstärkung abgeschätzt und detailliert untersucht. Die letztere Art von Feh-
ler wird in der wissenschaftlichen Literatur zu diesem Thema selten diskutiert. Tatsächlich kann der Einfluss der
diskreten Fehlerverstärkung auf den Gesamtfehler jedoch erheblich sein, insbesondere wenn Störungen aufgrund
von Transportfehlern in den Volumenfraktionsdaten vorhanden sind. Die kinematische Evolutionsgleichung für die
mittlere Krümmung, die im ersten Teil dieser Arbeit abgeleitet wird, dient als Referenzlösung zur Validierung des
numerischen Transports der Krümmung an der Kontaktlinie. Wie aus der Fehleranalyse zu erwarten ist, wird fest-
gestellt, dass der Transportfehler für die Krümmung für die Boundary Youngs-Methode divergent erster Ordnung
und für die Boundary ELVIRA-Methode konstant ist. Die Ergebnisse zeigen deutlich den Bedarf an Grenzflächen-
Advektionsmethoden höherer Ordnung.

Der dritte Teil dieser Arbeit untersucht zwei spezifische Beispiele, nämlich den Anstieg einer Flüssigkeitssäule
in einer Kapillare (”Capillary Rise“) und das Aufbrechen einer Flüssigkeitsbrücke auf einer chemisch strukturier-
ten Oberfläche. Ein vergleichender numerischer Benchmark für dynamische Benetzungsprozesse wird am Beispiel
des Capillary Rise entwickelt. Dazu kommen vier völlig verschiedene numerische Methoden zum Einsatz, die im
SFB 1194 entwickelt werden. Das Resultat ist eine fundierte Datenbasis zur Validierung künftiger numerischer
Verfahren. Darüber hinaus wird eine neuartige Diskretisierung der Navier-Schlupfbedingung eingeführt, die als

”versetzter Schlupf“ oder ”staggered slip“ bezeichnet wird. Das Ziel der versetzten Schlupfbedingung besteht dar-
in, den der Methode inhärenten ”numerischen Schlupf“ zu verringern. Dies wird erreicht, indem die Schlupflänge
in Bezug auf einen virtuellen Rand definiert wird, der sich zwischen dem physikalischen Rand und dem Ort der
flächenzentrierten Geschwindigkeit befindet, welche zum Transport des Volumenfraktionsfeldes verwendet wird.
Infolgedessen ist die diskrete viskose Dissipation im Vergleich zur Standard Navier-Schlupf-Implementierung
erhöht. Es wird gezeigt, dass die Konvergenz für den Kapillaranstieg signifikant verbessert werden kann, wenn die
Schlupflänge noch nicht aufgelöst ist. Andererseits ist die Konvergenzordnung für eine einphasige Kanalströmung
im Vergleich zur Standardimplementierung reduziert.

Die Benetzung strukturierter Oberflächen wird in Zusammenarbeit mit Experimentatoren im SFB 1194 untersucht.
Ziel ist es, die Aufbruchdynamik einer benetzenden Kapillarbrücke auf einer strukturierten Oberfläche quantitativ
zu beschreiben. Ein Hauptproblem für die Interpretation sowohl der experimentellen als auch der numerischen
Daten ergibt sich aus der Unsicherheit bezüglich des genauen Zeitpunkts des Aufbruchs der Kapillarbrücke. Um
dieses Problem zu lösen, fü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 Durchmes-
sers) 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 Standard-
darstellung auf die Phasenraumdarstellung abbildet, invertierbar ist, bis auf eine Verschiebung in der absoluten
Zeit. Mit dieser Methode können wir den Aufbruchprozess sowohl in der dreidimensionalen VOF-Simulation
als auch im Experiment sehr detailliert untersuchen. Im Allgemeinen wird eine gute Übereinstimmung zwischen
Experiment und Simulation sowohl qualitativ als auch quantitativ in Bezug auf die zeitliche Entwicklung des
minimalen Durchmessers erreicht. Die numerischen Simulationen ermöglichen es, verschiedene Regime in der
Aufbruchdynamik zu identifizieren, die auch in den experimentellen Daten gefunden wurden. Bemerkenswerter-
weise könnten dynamische Oberflächenspannungseffekte eine wichtige Rolle bei der Aufbruchdynamik spielen.
Die Übereinstimmung zwischen Simulation und Experiment kurz vor dem Aufbruch kann verbessert werden, in-
dem der Wert der (konstanten) Oberflächenspannung auf 90mN/m erhöht wird. Der letztere Wert wurde von
Hauner et al. (J. Phys. Chem. Lett., 2017) für eine frisch erzeugte Wasseroberfläche vorgeschlagen. Darüber hin-
aus wird durch Auswertung der numerischen Daten festgestellt, dass die lokale Grenzflächenerzeugungsrate sehr
hoch ist. Eine vollständig quantitative Bewertung dieses Phänomens steht jedoch noch aus. Schließlich betrachten
wir noch einmal das Capillary Rise Problem, diesmal für eine strukturierte Oberfläche. Interessanterweise kann
die Oberflächenstrukturierung dazu verwendet werden, ein Energiefunktional mit zwei stabilen Konfigurationen

x

https://www.doi.org/10.1017/jfm.2016.276
http://dx.doi.org/10.1021/acs.jpclett.7b00267


zu konstruieren. Das resultierende dynamische Aufstiegsverhalten ist sehr komplex und zeigt eine Bifurkation
zwischen den beiden stabilen Gleichgewichtslagen bei einer kritischen initialen Steighöhe.

xi


Contents
1. Fundamental concepts of wetting 1

1.1. Capillarity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2. Static Wetting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3. Variational derivation of the Young–Dupré equation in two dimensions . . . . . . . . . . . . . . . 7
1.4. Dynamic Wetting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

I. Mathematical Modeling & Qualitative Analysis 15

2. A geometrical model for spreading drops in the partial wetting regime 17
2.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2. Spreading dynamics of spherical drops . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3. Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.4. Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3. Sharp interface models 31
3.1. Notation and mathematical setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2. Entropy production . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3. The standard model for moving contact lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4. A note on the generalized Navier boundary condition . . . . . . . . . . . . . . . . . . . . . . . . 37
3.5. Interface Formation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

4. Kinematics of moving contact lines 39
4.1. Kinematic transport of the contact angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2. Contact angle evolution in the standard model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.3. Remarks on more general models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.4. Kinematic transport of the mean curvature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5. Boundary conditions for dynamic wetting - A mathematical analysis 65
5.1. Compatibility conditions for partial differential equations . . . . . . . . . . . . . . . . . . . . . . 66
5.2. Compatibility analysis for the standard Navier slip model . . . . . . . . . . . . . . . . . . . . . . 66
5.3. Compatibility analysis for a model by Lukyanov and Pryer . . . . . . . . . . . . . . . . . . . . . 69
5.4. Comparison with results by Lukyanov and Pryer . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.5. Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

II. Numerical Methods 77

6. Numerical methods for dynamic wetting: An overview 79
6.1. Overview of numerical methods for multiphase flows . . . . . . . . . . . . . . . . . . . . . . . . 79
6.2. Contact line advection using the level set method . . . . . . . . . . . . . . . . . . . . . . . . . . 82

7. The geometrical Volume-of-Fluid method 85
7.1. The VOF solver FS3D: An overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
7.2. Interface advection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7.3. Interface reconstruction methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
7.4. Numerical modeling of surface tension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

xiii


Contents

7.5. State of the art: Wetting with geometrical VOF methods . . . . . . . . . . . . . . . . . . . . . . . 93

8. Contact line advection using the geometrical VOF method 99
8.1. An analytical reference solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
8.2. Interface reconstruction close to the boundary in two dimensions . . . . . . . . . . . . . . . . . . 102
8.3. Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
8.4. Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

9. Curvature from Height Functions: A mathematical error analysis 113
9.1. Analytical error estimate in two dimensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
9.2. Computational examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
9.3. A note on higher-order methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

10. Numerical transport of the mean curvature 125
10.1. Reference solutions for the curvature transport in two dimensions . . . . . . . . . . . . . . . . . . 126
10.2. Volume-of-Fluid-based curvature transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

III. Applications 137

11. The capillary rise problem 139
11.1. Static capillary rise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
11.2. Continuum mechanical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
11.3. Dynamics of capillary rise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
11.4. Staggered slip boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149

12. Wetting of structured surfaces 155
12.1. Static wetting of an inhomogeneous substrate in two dimensions . . . . . . . . . . . . . . . . . . 155
12.2. Breakup dynamics of droplets on structured surfaces . . . . . . . . . . . . . . . . . . . . . . . . 161
12.3. A note on dynamic surface tension effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
12.4. Structured capillary rise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177

IV. Conclusion and Outlook 181

13. Concluding remarks 183

V. Appendix 185

A. Proof of the entropy production theorem 187

B. Kinematics of moving contact lines: Some mathematical details 191

C. Some geometrical relations 195

D. On the implementation of the Boundary ELVIRA method 197

Bibliography 201

Curriculum Vitae 215

xiv


List of Figures

Figure 1.1 Kinematics of wetting with no-slip. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

Figure 1.2 Water strider resting on a water surface and microscopic images of the animals leg. . . . . . . . 3

Figure 1.3 Surface tension as a force per unit length. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

Figure 1.4 Young diagram and equilibrium contact angle. . . . . . . . . . . . . . . . . . . . . . . . . . . 4

Figure 1.5 Sketch of the capillary rise problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

Figure 1.6 Wetting of a patterned surface with hydrophilic and hydrophobic stripes. . . . . . . . . . . . . 6

Figure 1.7 Height function representation of a sessile wetting droplet in two dimensions. . . . . . . . . . 7

Figure 1.8 Sketch of the channel flow problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

Figure 2.1 Experimental data for the spreading of a droplet. . . . . . . . . . . . . . . . . . . . . . . . . 17

Figure 2.2 Deposition of a pure water droplet on a silicon wafer. . . . . . . . . . . . . . . . . . . . . . 19

Figure 2.3 Kistler’s empirical function and the Cox-Voinov law for θeq = 0. . . . . . . . . . . . . . . . 20

Figure 2.4 Some basic geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

Figure 2.5 Experimental setup for the spreading experiments. . . . . . . . . . . . . . . . . . . . . . . . 23

Figure 2.6 Spreading of a 75%-glycerol droplet on a silicon wafer. . . . . . . . . . . . . . . . . . . . . 24

Figure 2.7 Experimental data for the droplet height compared to hcap(Lexp,Vexp). . . . . . . . . . . . . . 25

Figure 2.8 Experimental data for the contact angle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

Figure 2.9 Empirical relation for the dynamic contact angle. . . . . . . . . . . . . . . . . . . . . . . . . 26

Figure 2.10 ODE solution and literature comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Figure 2.11 Spreading of a water droplet on a polymer brush (80% humidity). . . . . . . . . . . . . . . . 27

Figure 2.12 Empirical function and ODE solution (80% humidity). . . . . . . . . . . . . . . . . . . . . . 27

Figure 2.13 Experimental data for the height compared to the ODE solution (80% humidity). . . . . . . . 28

Figure 2.14 Spreading of a water droplet on a polymer brush (50% humidity). . . . . . . . . . . . . . . . 28

Figure 3.1 Notation, local coordinate system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

Figure 3.2 Spreading droplet with an advancing contact line. . . . . . . . . . . . . . . . . . . . . . . . 36

Figure 4.1 Linear velocity fields satisfying Navier slip with L > 0 in a co-moving reference frame. . . . 53

Figure 5.1 Notation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

Figure 5.2 Roots of the mass flux ṁ1 for α̂1 = 1.6 and α̂2 = 0.51. . . . . . . . . . . . . . . . . . . . . . 73

Figure 5.3 Non-dimensional mass fluxes for α̂1 = 1.6 and α̂2 = 0.51. . . . . . . . . . . . . . . . . . . . 75

Figure 5.4 Dimensionless curvature and pressure at the contact line for α̂1 = 1.6 and α̂2 = 0.51. . . . . . 76

Figure 6.1 Different ways to represent the interface in multiphase flow simulations. . . . . . . . . . . . 79

xv


List of Figures

Figure 6.2 Numerical results for the two-dimensional contact line advection problem. . . . . . . . . . . 83

Figure 6.3 Convergence analysis for the two-dimensional advection problem. . . . . . . . . . . . . . . . 83

Figure 6.4 Numerical results for the three-dimensional contact line advection problem. . . . . . . . . . 84

Figure 7.1 PLIC interface and staggered grid for velocity components. . . . . . . . . . . . . . . . . . . 85

Figure 7.2 Local height function representation of the interface. . . . . . . . . . . . . . . . . . . . . . . 92

Figure 7.3 Ghost-cell based numerical realization of the no-slip and Navier slip boundary conditions
in FS3D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

Figure 7.4 Linear extrapolation of the height function as a numerical driving force for wetting. . . . . . 96

Figure 7.5 Linear extrapolation for θ = π/2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

Figure 8.1 Boundary Youngs reconstruction method for an equidistant mesh. . . . . . . . . . . . . . . 103

Figure 8.2 Boundary ELVIRA method on a 5×3-stencil in 2D. . . . . . . . . . . . . . . . . . . . . . 103

Figure 8.3 Translation of a straight line with a fixed orientation angle of θ = 60◦. . . . . . . . . . . . 104

Figure 8.4 Contact angle and contact line position from PLIC reconstruction. . . . . . . . . . . . . . . 105

Figure 8.5 Computational setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

Figure 8.6 Convergence with respect to the discrete L1-norm for the field (8.20) comparing initial and
final shapes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

Figure 8.7 Numerical motion of the contact line for the field (8.20) using the Boundary Youngs
reconstruction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Figure 8.8 Numerical motion of the contact line for the field (8.20) using the Boundary ELVIRA
reconstruction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Figure 8.9 Numerical contact angle evolution for the field (8.20) using the Boundary Youngs
method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

Figure 8.10 Numerical contact angle evolution for the field (8.20) using the Boundary ELVIRA
method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

Figure 8.11 Jump in the numerical contact angle for the field (8.20). . . . . . . . . . . . . . . . . . . . 109

Figure 8.12 Contact line evolution for the field (8.22). . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

Figure 8.13 Numerical contact angle evolution for the field (8.22). . . . . . . . . . . . . . . . . . . . . 110

Figure 8.14 Convergence behavior for the field (8.22) and different values of CFL. . . . . . . . . . . . 110

Figure 8.15 Contact line motion for the field (8.23). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

Figure 8.16 Contact angle evolution for the field (8.23). . . . . . . . . . . . . . . . . . . . . . . . . . . 111

Figure 8.17 Convergence behavior for the field (8.23) with the Boundary ELVIRA method. . . . . . . . 112

Figure 9.1 Finite differences based curvature for symbolically computed heights. . . . . . . . . . . . . 120

Figure 9.2 Convergence of the curvature for the Gauß-Legendre quadrature. . . . . . . . . . . . . . . 121

Figure 9.3 Convergence of the curvature for the trapezoidal rule. . . . . . . . . . . . . . . . . . . . . 122

Figure 9.4 Standard second-order finite differences vs. 4th-order method by Sussman and Ohta
[SO06]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

Figure 10.1 Setup for the VOF curvature transport. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

Figure 10.2 Stencil used to approximate f , f ′ and f ′′ at the boundary coordinate y = 0. . . . . . . . . . 129

xvi


List of Figures

Figure 10.3 Initial errors for the contact angle and the curvature - VOF and Level Set. . . . . . . . . . . 130

Figure 10.4 Initial curvature error. Higher-order [KB19] vs. standard second-order initialization. . . . . 130

Figure 10.5 Numerical transport of the contact point and contact angle for the field (10.15) with
Boundary Youngs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

Figure 10.6 Contact angle determined from PLIC and from a local height function. . . . . . . . . . . . 132

Figure 10.7 Numerical transport of the curvature for the field (10.15) with Boundary Youngs. . . . . . . 132

Figure 10.8 Numerical transport of the contact point and contact angle for the field (10.15) with
Boundary ELVIRA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

Figure 10.9 Numerical transport of the curvature for the field (10.15) with Boundary ELVIRA. . . . . . 134

Figure 10.10 Maximum error curvature transport T ∈ {0.4,1.2,2.4}. . . . . . . . . . . . . . . . . . . . 135

Figure 10.11 Numerical transport of the contact angle and the curvature for the field (10.19) with
Boundary ELVIRA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

Figure 10.12 Direct comparison Boundary ELVIRA vs. Level Set. . . . . . . . . . . . . . . . . . . . . . 136

Figure 11.1 Sketch of the capillary rise problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

Figure 11.2 Computational domain for the continuum mechanical model. . . . . . . . . . . . . . . . . 142

Figure 11.3 Simulation setup for FS3D. Length in cm . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

Figure 11.4 Mesh study with FS3D using numerical slip and resolved Navier slip. . . . . . . . . . . . . 146

Figure 11.5 Comparison of the methods for Ω = 1 with different slip lengths. . . . . . . . . . . . . . . 146

Figure 11.6 Numerical solutions for different values of Ω. . . . . . . . . . . . . . . . . . . . . . . . . 147

Figure 11.7 Convergence studies with FS3D for Ω=0.1, 0.5, 10, 100. . . . . . . . . . . . . . . . . . . . 148

Figure 11.8 Ghost-cell based numerical realization of the no-slip, Navier slip and staggered Navier slip
boundary conditions in FS3D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149

Figure 11.9 Capillary rise with the staggered slip boundary condition, results for L = R/5. . . . . . . . 151

Figure 11.10 Capillary rise with the staggered slip boundary condition, results for L = R/50. . . . . . . . 151

Figure 11.11 Illustration of the single-phase channel flow problem. . . . . . . . . . . . . . . . . . . . . 152

Figure 11.12 Signed relative error for the mean velocity 〈vx〉. . . . . . . . . . . . . . . . . . . . . . . . 153

Figure 11.13 Convergence of the mean velocity 〈vx〉. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

Figure 12.1 Wetting of a structured surface with a harmonic potential. . . . . . . . . . . . . . . . . . . 158

Figure 12.2 Considered configurations of the system. . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

Figure 12.3 The energy as a function of V0 for the configurations (a), (b) and (c) depicted in
Fig. 12.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161

Figure 12.4 Initial geometry from Surface Evolver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161

Figure 12.5 Time series of the breakup process top and side view. . . . . . . . . . . . . . . . . . . . . 162

Figure 12.6 Breakup dynamics for different choices of the breakup time t0 and the scaling exponent
ν . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165

Figure 12.7 Breakup dynamics of a capillary bridge - qualitative comparison. . . . . . . . . . . . . . . 168

Figure 12.8 Generation of satellite droplets in the VOF simulation. . . . . . . . . . . . . . . . . . . . . 168

Figure 12.9 Minimum bridge width as a function of time before breakup, experiment vs. simulation
(α = 1, L̃ = 500nm). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169

xvii


List of Figures

Figure 12.10 Breakup speed (α = 1, L̃ = 500nm) as a function of the bridge width compared to three
different experiments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170

Figure 12.11 Breakup speed and position of the local minimum relative to the center of the hydrophobic
stripe (numerical simulation data). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170

Figure 12.12 Comparison with numerical data from Li and Sprilles [LS16] for the breakup of a free
capillary bridge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171

Figure 12.13 Mesh and slip-length dependence study for α = 1, θphob = 102◦ and θphil = 28◦. . . . . . . 172

Figure 12.14 Influence of the dynamic viscosity and the staggered slip parameter ω . . . . . . . . . . . . 172

Figure 12.15 Influence of the wetting conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173

Figure 12.16 Breakup dynamics for α = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173

Figure 12.17 Breakup dynamics for α = 1.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174

Figure 12.18 Initial conditions from Surface Evolver compared to experimental images. . . . . . . . . . 174

Figure 12.19 Breakup dynamics for α = 1, simulation with increased surface tension. . . . . . . . . . . 175

Figure 12.20 Local interface generation rate during breakup (α = 1, σ = 72mN/m). . . . . . . . . . . . 176

Figure 12.21 Capillary rise problem for a structured solid surface. . . . . . . . . . . . . . . . . . . . . . 177

Figure 12.22 Approximation of the energy functional (Example 12.10). . . . . . . . . . . . . . . . . . . 178

Figure 12.23 Dynamics of the structured capillary rise (L/R = 1/5). . . . . . . . . . . . . . . . . . . . . 179

Figure 12.24 Dynamics of the structured capillary rise with increased slip length (L/R = 2/5). . . . . . . 180

Figure D.1 Illustration of the setup and the coordinate transformation. . . . . . . . . . . . . . . . . . . 197

Figure D.2 The function ψ for different values of β . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199

Figure D.3 Setup for the Boundary ELVIRA method on a 5×3-stencil in 2D. . . . . . . . . . . . . . . 200

xviii


List of Symbols
R Real numbers
N Natural numbers
Ω Computational domain
∂Ω Boundary of the computational domain
Σ Liquid-gas interface
κ Mean curvature
Γ Contact line
W Wetted area of the solid boundary
nΣ Interface normal vector
n∂Ω Outer unit normal to the domain
nΓ Contact line normal vector
tΓ Tangent vector to the contact line
PΣ Orthogonal projection onto the tangent space of Σ

P∂Ω Orthogonal projection onto the tangent space of ∂Ω

v Velocity vector
vΣ Interface velocity vector
VΣ Speed of normal displacement
VΓ Contact line normal speed
Φ Flow map
E Energy functional
η Dynamic viscosity
ρ Mass density
g Gravitational acceleration
p Pressure
σ Surface tension
T Cauchy stress tensor
S Viscous stress tensor
D Rate-of-deformation tensor
θ Contact angle
L Slip length
a Inverse slip length
Ca Capillary number
Bo Bond number
Oh Ohnesorge number
Re Reynolds number
χ Indicator function of the liquid phase
α Volume fraction
CFL Courant number

xix


Acronyms
CFL Courant-Friedrich-Lewy. 82

CSF Continuum Surface Force. 85, 91

FS3D Free Surface 3D. 85

GNBC Generalized Navier boundary condition. 37

IFM Interface Formation Model. 37

MD Molecular Dynamics. 10

MKT Molecular Kinetic Theory. 10

ODEs Ordinary Differential Equations. 10

PDEs Partial Differential Equations. 10

PLIC Piecewise Linear Interface Calculation. 85

VOF Volume-of-Fluid. 80

xxi


1. Fundamental concepts of wetting
“The force exerted on the solid surface is logarithmically infinite: not even Herakles could sink a solid
if the physical model were entirely valid, which it is not.” [HS71, p. 94]

This famous statement by Huh and Scriven in their seminal paper on the moving contact line problem [HS71]
marks the starting point of a history of research lasting for almost 50 years. The so-called “moving contact line
paradox” refers to the observation that the classical no-slip condition at solid boundaries – a cornerstone of fluid
mechanics – seems to be incompatible with dynamic wetting phenomena. By the term “dynamic wetting”, we
usually mean the motion of a droplet in a surrounding gas or a bubble in an ambient liquid along a solid boundary.

Figure 1.1.: Kinematics of wetting with no-slip.

The moving contact line paradox: The essence of the paradox can be easily described in terms of the pure
kinematics of the flow, i.e. without considering the details of the forces acting on the fluid particles. In the class of
“sharp interface” continuum mechanical models, the fluid-fluid interface is modeled as a mathematical hypersur-
face with zero thickness. In the simplest setting, the hypersurface is a “material surface” for the flow, i.e. it can be
described by a collection of marker particles that stay on the moving interface Σ(t) while being transported with
the flow according to the ordinary differential equation

ẋ(t) = v(t,x(t)). (1.1)

Here, v is the fluid (bulk) velocity, which is assumed to be continuous across the moving interface while its deriva-
tives may jump at the moving interface (see Chapter 3 for the details on the mathematical modeling). Provided
that v is locally Lipschitz continuous and satisfies the Dirichlet (“no-slip”) boundary condition

v|∂Ω = 0, (1.2)

it is clear that (1.1) is locally well-posed1. Physically, the local Lipschitz continuity of the velocity field in the
respective phases is related to the phenomenon of viscous dissipation of energy2. So, the obvious problem is that
– provided that v obeys the regularity described above – a marker particle located at the “contact line” Γ(t), i.e. the
line of intersection of the free surface and the solid boundary, stays at rest and the contact line as a whole cannot
move; see Figure 1.1. Since this property of the solution is clearly incompatible with experimental observations,
there must be some modification to either the model or the solution concept or a combination of both. The quote
above by Huh and Scriven refers to a hypothetical generalized “solution” of the problem which, however, predicts

1Note that the well-posedness of (1.1) has been proven recently in a much more general setting including mass transfer across the interface
[Bot20b].

2In fact, the rate of viscous dissipation in the bulk is proportional to
∫

Ω\Σ D : DdV , where D = 1
2 (∇v+∇vT) is the rate-of-deformation tensor

(see Section 3.2). Hence, the Lipschitz continuity of v is a sufficient (albeit not necessary) condition for a finite viscous dissipation rate.

1


Chapter 1. Fundamental concepts of wetting

an infinite force to move the contact line.

A generalization of the solution concept, i.e. the required regularity of the solution, leads to the study of hy-
drodynamic singularities of different kinds at the moving contact line. The physical implications of singularities
at the moving contact line are discussed controversially in the scientific literature. Eggers and Fontelos understand
singularities as “the fingerprint of a nonlinear PDE” since the structure of the singularity “is universal, i.e. indepen-
dent of the initial or boundary conditions imposed over macroscopic distances” [EF15, p.4]. The latter perspective
motivates to study the singularity itself in great detail to explore the physics of the process. The monograph by
Eggers and Fontelos [EF15] provides an introduction to that line of research applied to various phenomena in hy-
drodynamics.

A contrary position raises the question about the validity of the continuum mechanical model once a singular-
ity in one of the field variables or its derivatives is formed. In particular, Shikhmurzaev formulates the following
criterion for a mathematical theory of fluid flow:

“(B) The solution remains within the limits of applicability of the model specified via the assumptions
made in its formulation. In particular, the flow parameters, that is the components of velocity and
pressure, remain finite everywhere in the flow domain and on its boundaries at all times.”
[Shi06, p. 124]

The latter requirement is justified by the author by the range of validity of the modeling framework itself (see also
the monograph [Shi08] for a comprehensive discussion). More specifically, he states that

”The only element of criterion B we are interested in here is the finiteness of the flow parameters, that
is regularity of the solution. This is more restrictive than just integrability and does not allow solutions
with the integrably singular pressure, which causes no problems mathematically, but physically means
that the pressure in a finite vicinity of the singularity can be arbitrarily high (or low) and hence the
model of an incompressible fluid is no longer valid. Microscopically, this would also imply infinite
energies of molecules.” [Shi06, p. 124]

One of the main contributions of the present work is the development of a mathematical tool that allows to study
qualitative properties of mathematical models of dynamic wetting, provided that the solution is sufficiently reg-
ular. The “kinematic evolution equation for the dynamic contact angle” arises from the systematic study of the
kinematics of the moving contact line as sketched above. Besides its implications for different models of dynamic
wetting, it has also lead to novel numerical methods which are discussed in Part II of the present work. We make
no attempt to decide between the two contrary positions outlined above. This can ultimately only be done based
on careful experimental observations for specific flow configurations. We rather aim at a new perspective on the
field of dynamic wetting, which might lead to new insights in future research.

1.1. Capillarity
“A liquid surface can be thought of as a stretched membrane characterized by a surface tension that

opposes its distortion.” [dGBWQ04, p.1]

The mathematical study of capillarity is an old subject originating from the work of Pierre Simon de Laplace
and Thomas Young [You05] at the beginning of the 19th century, when the famous Young-Laplace equation of
capillarity was derived. The main physical ingredient for capillary phenomena is the surface tension of liquid-
liquid or liquid-gas interfaces. Microscopically, it arises from the attractive (also called “cohesive”) interactions
between liquid molecules. The molecules in the bulk phases are in an energetically favorable state due to the
attractive interaction with all their neighbors. In contrast, a molecule at the surface misses (approximately) half its
attractive interaction partners. Hence, it is energetically unfavorable for a liquid molecule to sit at the surface. As
a result, there is a finite amount of energy per unit area denoted σ > 0, which is necessary to create a unit portion
of the free surface. The latter concept of surface tension as surface energy allows for a variational approach to
capillary phenomena. Note that being an energy per unit area, the surface tension σ has the SI-units

[σ ] =
J

m2 =
N
m
.

2


1.1. Capillarity

Usually, gravity can be neglected if the length scale of the problem is well below the capillary length

lc =
√

σ/(ρg). (1.3)

Here, ρ = ρg− ρl is the density difference between the two phases and g is the gravitational acceleration. In
the absence of gravitational effects, the stationary shape of the capillary surface can be found by minimizing the
functional

E =
∫

Σ

σ dA

subject to appropriate constraints which are typically determined by the geometry of the problem. If the surface
energy σ is constant, the functional can be rewritten as

E = σ |Σ|.

Figure 1.2.: Water strider resting on a water surface (A) and microscopic images of the animals leg (B,C).
Reprinted (adapted) with permission from Feng et al., Langmuir 23(9):4892-4896, 2007.
Copyright 2007 American Chemical Society.

The simplest example is that of a water drop in zero gravity, subject only to the constraint of volume conservation
if incompressibility is assumed. The equilibrium shape of the drop is that of a sphere, minimizing the surface
area for a given volume. A more interesting example is that of a soap film that is geometrically constrained by
some kind of wireframe. In this case, the mathematical problem is to find a hypersurface with minimal area
subject to a prescribed boundary curve. The latter problem, known as the “Plateau-Problem” in mathematics,
fostered a lot of research in abstract differential geometry and related fields. Moreover, the effect of capillarity
is present in various instances in nature and technology. An interesting example from nature is the water strider;
see Figure 1.2 (A). The water strider is a small insect that is able to stand and walk easily on a water surface,
while its weight is supported almost exclusively by surface tension [BH06, p.344]. Remarkably, the legs of the
insect feature a specialized hierarchical structure consisting of oriented needle-shaped microsetae (i.e. tiny hairs)
that themselves possess elaborate nanogrooves [FGW+07, p.4892]; see Fig. 1.2 (B) and (C). As a result, the legs
show superior water repellency and a single leg is able to hold a force of about 15 times the weight of the insect’s
body [FGW+07, p.4893]. A comprehensive review of capillary and wetting phenomena in nature and technology
can be found in the monograph of De Gennes et al. [dGBWQ04].

Mechanical formulation as a surface tension force: It is important to note that there is a second equivalent
formulation of surface tension as a force per unit length acting along the interface. For example, consider a planar
soap film that is placed between a rectangular wireframe and a movable rod of length L, as sketched in Fig. 1.3.
Suppose we are able to measure the force on the rod as it is pulled outwards. The surface tension force acts parallel
to the surface and aims to reduce the surface area of the soap film. The surface tension force acting on the rod is

f =−2σLêx,

where the factor 2 appears because two interfaces are present. Obviously, the variational and the mechanical
formulation are equivalent since the work done to move the rod by an amount ∆x is

∆W = 2σL∆x = 2σ∆A.

3

http://dx.doi.org/10.1021/la063039b


Chapter 1. Fundamental concepts of wetting

Figure 1.3.: Surface tension as a force per unit length.

1.2. Static Wetting
In the context of the present work, we mean by wetting a process where a two-phase flow (typically constituted
by a liquid and a gas phase) interacts with a solid surface; see Fig. 1.4. We assume that the solid is rigid, i.e. we
do not account for deformations of the solid resulting from the wetting process. This is expected to be a good
approximation for surfaces with a sufficiently large elastic modulus. We note, however, that the deformation of the
substrate can be relevant in principle, in particular in the field of “soft wetting” (see, e.g., [AS20]).

Figure 1.4.: Young diagram and equilibrium contact angle.

The Young–Dupré equation: The Young–Dupré equation was established in the early 19th century in the pio-
neering work of Thomas Young [You05]. It describes the equilibrium state of a wetting drop on a flat and homo-
geneous solid surface in terms of the local equilibrium contact angle θeq, i.e. the angle of intersection between the
free surface and the solid boundary as depicted in Fig.1.4 according to

σlg cosθeq+σw = 0. (1.4)

Here, the quantity σw is defined as the difference between the surface tensions σls for the liquid-solid interface and
σsg for the solid-gas interface, i.e.

σw := σls−σsg.

Consequently, it can be understood as the specific energy of the wetted solid surface relative to the “dry surface”.
The latter quantity serves as a driving force for the wetting process. If the “specific wetting energy” σw is lower
than the specific energy σlg of the liquid-gas interface, the system can gain energy by converting liquid-gas into
liquid-solid interface by wetting the solid.

The usual way to read the Young equation (1.4) and the Young diagram in Fig. 1.4 is, however, based on the
mechanical definition of the surface tension as a force per unit length. Then, eq. (1.4) is simply obtained from a
horizontal projection of the forces depicted in Fig. 1.4. It is important to note that the interfacial tensions σls and
σlg can usually not be measured independently. Instead, eq. (1.4) may be used to determine σw from measured
values of σlg and θeq. Remarkably, molecular dynamics simulations indicate that the Young equation (with the
respective local values of the interfacial tensions) holds down to a length scale of a few nanometers [FTBD17]. On
the other hand, there is some debate in the literature whether or not the interpretation of (1.4) as a local balance of

4


1.2. Static Wetting

forces is actually valid [Mak16, Mak18]. In particular, the interpretation of the liquid-solid and solid-gas surface
tension as a force is discussed controversially.

It is, therefore, of interest to give a variational justification of (1.4) which only resorts to the concept of surface
energies. An elementary derivation in two dimensions, which also serves as an introduction to some fundamental
mathematical concepts, is given in Section 1.3. We note that the variational approach also allowed to extend the
theory to the case of rough solid substrates; cf. [Wen36, CB44].

Static wetting model for a homogeneous substrate: In order to determine the equilibrium state, we consider
the energy functional

E = σlg |Σlg|+σw |Σsl |+Eg, (1.5)

where σlg > 0 is the surface tension3 of the liquid-gas interface, σw = σsl−σsg is the specific energy of the wet
solid surface and Eg is the gravitational energy. The equilibrium state of a droplet is a minimizer of the functional
(1.5) under the constraint of volume conservation.

Remark 1.1. From the form of the functional (1.5) one may recognize the following properties:

(i) Clearly, in the absence of gravity and without contact to a solid surface, the equilibrium states of the system
are spheres with a radius determined by the volume. This is a consequence of the fact that the shape of a
sphere minimizes the area enclosing a given volume.

(ii) If the specific wetting energy σw is larger than the liquid-gas surface tension, i.e. for

Wa := σlg−σw < 0, (1.6)

the system does not gain energy by transforming a piece of liquid-gas interface into liquid-solid interface by
wetting the surface. In this case, the fluid is called non-wetting (for the considered solid surface). Otherwise,
the fluid is called wetting.

The quantity Wa defined above is known in the literature as the work of adhesion. It is “defined as the
reversible thermodynamic work that is needed to separate the interface from the equilibrium state of two
phases to a separation distance of infinity” [Ebn11].

(iii) If, one the other hand, the specific wetting energy σw is less than−σlg, i.e. if σlg+σw < 0, a simple argument
shows that no local minimum with finite length L exists and the liquid forms a film: Consider a configuration
of the liquid arranged in a disc of height h and radius R. Clearly, the height is a function of the radius to
satisfy the volume conservation constraint, i.e.

h =
V

πR2 .

If the disk is sufficiently flat, the gravitational energy can be neglected and the energy can be expressed as a
function of the radius according to

E (R) ∝ πR2(σlg+σw)+σlg
2V
R

.

Consequently, the energy is not bounded from below if σlg+σw < 0 (or, equivalently, if Wa > 2σlg) and,
hence, no stable configuration with a finite radius exists in this case. Fluids of this class are called perfectly
wetting for the solid substrate.

(iv) If the work of adhesion is in between the extreme cases describe above, i.e. if the fluid is neither non-wetting
nor perfectly wetting

0≤Wa ≤ 2σlg ⇔ σw

σlg
∈ [−1,1],

the fluid is said to be partially wetting for the solid substrate. Mathematically, there is a unique solution to
the Young–Dupré equation in this case.

5


Chapter 1. Fundamental concepts of wetting

Figure 1.5.: Sketch of the capillary rise problem.

Example 1.2 (Capillary rise). One of the most prominent examples of static wetting is the rise of liquid in a
capillary tube with a diameter comparable to (or smaller than) the capillary length (1.3). The rise height H of
the liquid in the equilibrium state (see Fig. 1.5) can be found by means of energy considerations and is strongly
dependent on the wettability of the solid. Indeed, the gain in energy due to wetting is the driving force of the
process which is opposed by gravity. The height given by

H =
2σlg cosθeq

ρgR
(1.7)

is well-known in the literature as “Jurin’s height”, see, e.g., [dGBWQ04] and Chapter 11 for a more detailed
discussion of the statics and dynamics of capillary rise.

Example 1.3 (Static wetting on structured surfaces). Much more complex wetting patterns arise if the model (1.5)
is generalized to the case of an inhomogeneous solid substrate. In the latter case, the specific wetting energy may
be a function of the position on the solid substrate and the energy functional generalizes to

Ei = σlg |Σlg|+
∫

Σsl

σw dA+Eg. (1.8)

Figure 1.6 shows an example of static and dynamic wetting on a structured surface with a stripe geometry. The
surface has been prepared with photolithographic methods such that stripes with a different surface chemistry, in
particular with a different wettability, are produced.

(a) Surface Evolver calculation. (b) Experiment - early stage. (c) Experiment - close to breakup.

Figure 1.6.: Wetting of a patterned surface with hydrophilic and hydrophobic stripes(Figures by M. Hartmann).

The geometry shown in Fig. 1.6(a) has been calculated with the Open Source tool Surface Evolver [Bra92,
BB12], which is able to numerically minimize the functional (1.8). The droplet (blue color) prefers to wet the hy-
drophilic parts of the substrate (green color), leading to two elongated segments that are connected by a capillary

3Note that the liquid-gas surface tension is always assumed to be positive since otherwise, no stable interfaces exist.

6


1.3. Variational derivation of the Young–Dupré equation in two dimensions

bridge over the hydrophobic region. Note that this shape has a much larger gas-liquid interfacial area compared
to a spherical cap of equal volume. However, the contact area between the liquid and the hydrophobic part of the
substrate is minimized such that the overall energy is lower compared to a spherical cap of equal volume. Depend-
ing on the total volume, it can be energetically favorable for the droplet to break up into two individual droplets,
wetting only the hydrophilic regions. In this case, the capillary bridge becomes unstable and a dynamic breakup
process is initiated [HH19]; see Fig. 1.6(b) and 1.6(c). The dynamics of this breakup process is investigated in
detail with continuum mechanical simulations in Chapter 12 (see also [HFW+20]).

1.3. Variational derivation of the Young–Dupré equation in two
dimensions

Let us start the discussion of the physics of wetting with a proper derivation of (1.4) based on variational principles.
We consider the simplified two-dimensional case and neglect the effect of gravity. Even though this is a rather
limited setting compared to real examples of static wetting, it is still useful to introduce and understand some
basic mechanisms of wetting. For a more detailed mathematical discussion of equilibrium capillary surfaces see,
e.g. [Fin86, Bra92, BB12].

Formulation via height functions: Suppose that the liquid-gas interface in two dimensions can be described as
a graph over the solid boundary, i.e.

Σlg = {(x,h(x)) : 0≤ x≤ L};
see Fig. 1.7. Note that this approach obviously fails for θ ≥ π/2 and another representation of the interface is nec-
essary in this case. However, we leave this case aside for the purpose of this introduction to the topic. Moreover,
we only consider the case of a single droplet and do not consider the process of breakup of a droplet into smaller
droplets (see Chapter 12).

Figure 1.7.: Height function representation of a sessile wetting droplet in two dimensions.

If the interface can be represented by a height function, we may express the energy functional as

E (L,h) = σlg

∫ L

0

√
1+h′(x)2 dx+Lσw+

ρg
2

∫ L

0
h(x)2 dx, (1.9)

where g is the gravitational acceleration and
ρ = ρl−ρg.

We are looking for minima of E in the configuration space

L≥ 0, h ∈ C 1(0,L) with h(x)≥ 0, h(0) = h(L) = 0

subject to the constraint ∫ L

0
h(x)dx =V0,

where V0 > 0 is the prescribed volume of the drop.4

4Note that the energy of the sessile drop is invariant with respect to translations along the boundary since the solid is assumed to flat and
homogeneous. We, therefore, do not need to consider translations.

7


Chapter 1. Fundamental concepts of wetting

Non-dimensional form of the problem: In order to proceed, it is convenient to reformulate the problem in non-
dimensional form, which allows to define the height function h on a fixed interval. The importance of the latter
approach is to allow for an independent variation of the length L and the height function itself. We define

x̃ = x/L, h̃(x̃) = h(x̃L)/L.

Hence, we have the relations
h̃′(x̃) = h′(x̃L), h̃′′(x̃) = Lh′′(x̃L).

Then the transformed minimization problem reads as (dropping the tilde notation):

Given V0,σlg,ρ,g > 0 and σw ∈R, minimize the functional

E = L
(

σlg

∫ 1

0

√
1+h′(x)2 dx+σw

)
+

ρgL3

2

∫ 1

0
h(x)2 dx

= Lσlg

(∫ 1

0

√
1+h′(x)2 dx+

σw

σlg
+

Bo
2

∫ 1

0
h(x)2 dx

) (1.10)

over the configuration space

L≥ 0, h ∈ C 1(0,1) with h(x)≥ 0, h(0) = h(1) = 0

subject to the constraint ∫ 1

0
h(x)dx =

V0

L2 . (1.11)

We note that the non-dimensional Bond number (also known as Eötvös number Eo), defined as5

Bo =
ρgL2

σlg
, (1.12)

measures the importance of gravitational forces compared to surface tension forces. Since in the present work we
are mainly interested in capillary (i.e. surface tension dominated) flows, the Bond number will typically be small.

Euler-Lagrange Equations in the absence of gravity: We now proceed in the simplified case where gravity is
negligible, i.e. we assume Bo = 0. In order to satisfy the volume constraint, we introduce a Lagrange multiplier λ

and consider the functional

Eλ = Lσlg

(∫ 1

0

√
1+h′(x)2 dx+

σw

σlg

)
+λ

(∫ 1

0
h(x)dx− V0

L2

)
, (1.13)

which is now to be minimized over

L > 0, λ ∈R, h ∈ C 1(0,1) with h(x)≥ 0, h(0) = h(1).

By differentiating (1.13) with respect to L, we find the necessary condition

0 = σlg

(∫ 1

0

√
1+h′(x)2 dx+

σw

σlg

)
+

2λV0

L3 . (1.14)

We now take a smooth function ϕ ∈ C ∞
c (0,1) with compact support in (0,1) and consider the Fréchet-derivative

of Eλ , i.e. we consider the function

gϕ(ε) := Eλ (h+ εϕ) = Eλ = Lσlg

(∫ 1

0

√
1+(h′+ εϕ ′)2 dx+

σw

σlg

)
+λ

(∫ 1

0
(h+ εϕ)dx− V0

L2

)
5In this case the foot-length of the droplet is chosen as the typical length scale. In a more general definition, one may define Bo := ρgl2/σlg,

where l is a characteristic length scale of the problem.

8


1.3. Variational derivation of the Young–Dupré equation in two dimensions

and compute its derivative at the origin

g′ϕ(0) = Lσlg

∫ 1

0

h′ϕ ′√
1+h′2

dx +λ

∫ 1

0
ϕ dx.

Using partial integration, we obtain

g′ϕ(0) = Lσlg
h′ϕ√
1+h′2

∣∣∣1
0︸ ︷︷ ︸

=0

+
∫ 1

0

(
−Lσlg

(
h′√

1+h′2

)′
+λ

)
ϕ dx. (1.15)

By the fundamental lemma of the calculus of variations, the necessary conditions g′ϕ(0) = 0 for all ϕ ∈ C ∞
c (0,1)

only hold if the integrand in (1.15) vanishes, i.e. for a solution of

Lσlg

(
h′√

1+h′2

)′
= λ ∀x ∈ (0,1), h(0) = h(1) = 0.

A short calculation shows the relation(
h′√

1+h′2

)′
=

h′′

(1+h′2)3/2 = κ̃(x).

Here κ̃ is the curvature of the curve (x̃, h̃(x̃)), i.e. the curvature of the interface in the transformed coordinates.
Note that we may rewrite κ̂ in terms of the curvature of the original interface thanks to

κ̃(x̃) =
h̃′′(x̃)√

1+(h̃′(x̃))2
=

Lh′′(x̃L)
(1+h′(x̃L)2)3/2 = Lκ(x̃L)

Hence, we have the condition

κ(x) =
λ

L2 σlg
= const. (1.16)

So, in fact, the Lagrange multiplier introduced to satisfy the volume constraint turns out to be the curvature of the
interface. Moreover, we conclude that the minimizer of the functional is a curve with constant curvature, i.e. a
circular cap. Combining the expression (1.16) for λ with the stationarity condition (1.14) yields

0 = σlg

∫ 1

0

√
1+ h̃′(x̃)2 dx̃+σw+

2σlg κV0

L

= σlg

(∫ 1

0

√
1+ h̃′(x̃)2 dx̃+

2κV0

L

)
+σw

= σlg

(
1
L

∫ L

0

√
1+h′(x)2 dx+

2κV0

L

)
+σw

=
σlg

L

(
|Σlg|+2κV0

)
+σw .

(1.17)

Some elementary geometry shows the following relations for a circular cap

κ =− 1
R
, V0 = θR2− LR

2
cosθ , |Σlg|= 2θR,

where R is the radius of the circle and θ is the contact angle. This immediately shows that

|Σlg|+2κV0 = Lcosθ

and (1.17) is in fact equivalent to the Young–Dupré equation (1.4). In summary, we have shown that the minima
of the energy functional are spherical caps with volume V0 and a contact angle that respects (1.4). In particular,
the variations in the foot-length L lead to the Young–Dupré equation (1.4) and the variations in the rescaled height
function h̃ lead to the constant curvature condition (1.16).

9


Chapter 1. Fundamental concepts of wetting

1.4. Dynamic Wetting

1.4.1. Hierarchy of models
Even though the statics of wetting is well-understood for a long time, the mathematical modeling of dynamic
wetting is still a challenge. The research efforts over the last decades lead to a lot of scientific debate and a
great variety of different mathematical models. For a recent survey on the field, see, e.g., [dGBWQ04], [Bla06],
[BEI+09], [Shi08], [SA13] and the references therein. The purpose of this section is to give a brief (and by no
means complete) overview and to put the present work into a broader context.

From the perspective of mathematical modeling, there is a whole hierarchy of models to describe dynamic wet-
ting, ranging from microscopic descriptions via continuum mechanics to simplified/empirical descriptions of the
macroscopic flow. Possible frameworks to describe dynamic wetting are:

(i) Molecular Dynamics,

(ii) Continuum Mechanics and

(iii) Simplified/Empirical descriptions.

Molecular Dynamics: The Molecular Dynamics (MD) point of view accounts for the discrete structure of mat-
ter and models the wetting problem by a large collection of point particles. Typically, the system contains two
types of particles, namely fluid particles and solid particles which are arranged in some lattice structure to model
the solid substrate. The physical interaction between a pair of fluid-fluid or fluid-solid particles is modeled by an
effective interaction potential which is frequently chosen to be the Lennard-Jones potential leading to so-called
“Lennard-Jones liquids”. Once the interaction potentials are modeled, the main task is to solve Newtons equations
of motion which leads to a large system of Ordinary Differential Equations (ODEs) that determines the trajectories
of the particles. The advantage of the molecular dynamics approach is that it allows to study the physical processes
at a moving contact line in great detail. However, the method is limited by the available computational resources.
In order to be computationally feasible, MD simulations are often restricted to very small length and time scales.
Moreover, the computational feasibility may also limit the range in which physical parameters (such as tempera-
ture, density, viscosity or surface tension) can be chosen in practice. Then the question arises how the MD results
translate to different parameter regimes. Despite these limitations, the MD approach gives important insights into
physical processes on the microscale which may lead to closure relations or boundary conditions for higher-level
models such as Navier Stokes (see, e.g., [QWS03, RHE10, CWQS15]). The Molecular Kinetic Theory (MKT) of
dynamic wetting by Blake and Hayes [BH69], being one of the first quantitative theories of dynamic wetting on
the microscale, still leads to new insights into the physics of wetting. For example, recent work shows that key
parameters of dynamic wetting can be obtained from thermal fluctuations at equilibrium [FTBDC19].

Continuum Mechanics: The continuum mechanical approach employs Partial Differential Equations (PDEs)
that govern the temporal evolution of macroscopic quantities, such as the mass density and the average particle
velocity. Typically, the latter is done in the framework of the Navier Stokes equations in a free-surface or sharp-
interface two-phase formulation. The relevant physical effects on the microscale are modeled through constitutive
equations and boundary conditions. The latter is done on a (semi-) empirical basis and/or motivated by more
fundamental physical models. While the continuum mechanical models have shown their ability to describe a
large number of wetting flows quite accurately, they may still be demanding in terms of computational costs. This
motivates to look for simplified continuum mechanical models. A frequent choice to obtain a simplified continuum
mechanical model is the lubrication approximation. The lubrication approximation applies in a geometry where
one dimension is significantly smaller than the others. In this case, the Navier Stokes system can be simplified
and the two-phase flow problem can be reduced to a PDE that only involves the free surface shape (see, e.g,
[Sno06, Lea07, EF15]).

Simplified and Empirical descriptions: A substantial amount of research is dedicated to finding empirical cor-
relations describing some macroscopic features of the flow with only a few parameters. The latter type of descrip-
tion might be useful to design systems for practical applications. A typical application example is curtain coating.

10


1.4. Dynamic Wetting

It is known that undesired air inclusions are generated if the coating speed is too high. The latter phenomenon
is called “dynamic wetting failure” in the literature. The maximum speed at which the transition to the wetting
failure occurs is of great practical interest. Experimental data [BBS99] can be a basis to find some empirical math-
ematical formula that describes the maximum coating speed. According to Weinstein and Ruschak, “genuinely
predictive modeling of complex coating processes is not yet possible and coating practice remains largely empiri-
cal. Nonetheless, coating science is sufficiently advanced that physical insights and mathematical models greatly
benefit design and practice.” [WR04] However, recently the wetting failure in certain coating became accessible
also to detailed numerical simulations [LVCK16]. Another example of a simplified model of a dynamic wetting
process is the ODE model for the spreading of thin perfectly wetting drops described in [dG85]. A generalization
of the latter model to partial wetting is discussed in Chapter 2.

Hybrid models: To conclude this short tour of mathematical models, we note that there are also hybrid ap-
proaches that combine some of the general frameworks listed above. For example, Diffuse Interface models employ
the continuum mechanical description but model the interface as a transition region of a small but finite thickness
(see [AMW98] for a review). The technical advantage of this approach is that the governing equations can be for-
mulated on a fixed domain and discontinuities of the field variables are avoided. The physical motivation for diffuse
interface models is the fact that real fluid interfaces have a small but finite thickness, determined by molecular in-
teractions. According to Nold, diffuse interface models can be categorized as mesoscopic models “because they
employ ideas from the nanoscale, and are used to compute multiphase scenarios at the macroscale.” [Nol16, p.57].

1.4.2. Motion of viscous fluids: The Navier Stokes equations
The basis of the continuum mechanical description is formed by the classical Navier Stokes equations. The Navier
Stokes equations are a system of PDEs that describe the conservation of mass and linear momentum for a viscous
fluid. Starting from the continuity equation

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

one immediately arrives at the condition

∇ · v = 0 (1.18)

if the flow is considered incompressible, i.e. if the density of fluid particles is required to be constant. The momen-
tum balance equations, which are obtained from Newtons laws of motion and an appropriate closure relation for
the stresses in the fluid, read as

ρ(∂tv+ v ·∇v)−η∆v+∇p = ρb. (1.19)

The vector field b describes external forces that act on the fluid such as gravity, the scalar field p is called the
pressure and η is called the dynamic viscosity of the fluid. The dynamic viscosity is a measure for the internal
friction of the fluid which leads to a resistance of the medium to shear stresses. We note that (1.19) is obtained
from a constitutive law that relates the shear linearly with the stress in the fluid. The latter assumption may
be generalized to study so-called non-Newtonian fluids. We do not consider Non-Newtonian fluids in the present
work, but note that the breakdown of the linear closure model is one of the suggested ways to regularize the moving
contact line singularity (see, e.g., [AG02]). The system of equations (1.19) and (1.18) must be accompanied by
appropriate boundary and initial conditions. The most common choice for the velocity of a viscous fluid at a solid
boundary is the no-slip condition which states that the fluid particles right next to the solid boundary do not move
relative to the boundary. Mathematically, the no-slip condition is a homogeneous Dirichlet boundary condition,
i.e.

v = 0 at ∂Ω. (1.20)

1.4.3. Continuum mechanical models of dynamic wetting
In 1971, Huh and Scriven [HS71] showed that the usual no-slip condition for the velocity at solid boundaries is
not appropriate for the modeling of moving contact lines (see also [PS82], [Sol97]). Depending on the model

11


Chapter 1. Fundamental concepts of wetting

and the solution concept, this means that the no-slip condition either leads to the non-existence of solutions with
moving contact lines or to an infinite viscous dissipation rate. Since then, many attempts have been made to
solve this problem, for instance by means of numerical discretizations6 (for an overview over numerical methods
see [SDS14]) and/or by replacing the boundary conditions in the model. Essentially, some mechanism is introduced
which allows for a tangential slip of the interfacial velocity at the solid wall.

Slip at the fluid-solid interface: The most common choice for the velocity boundary condition is the Navier
slip condition, already proposed by Navier in the 19th century, which relates the tangential slip to the tangential
component of normal stress at the boundary according to

v‖+L(Dn∂Ω)‖ = 0, (1.21)

where D = 1
2 (∇v+∇vT) is the rate-of-deformation tensor and L > 0 is called slip length.

Figure 1.8.: Sketch of the channel flow problem.

Example 1.4 (Single-phase channel flow with Navier Slip). Consider the simple example of the steady flow of a
Newtonian fluid described by the incompressible Navier Stokes equations through a two-dimensional channel of
width H; see Fig. 1.8. The flow is driven by a prescribed, spatially constant pressure gradient ∂x p = −G. The
Navier slip condition is prescribed at the solid boundaries at y = 0 and y = H. So, we are looking for a solution of
the Navier Stokes equations (1.19) - (1.18) with b = 0 of the form

p(x,y) = p0−Gx, v(x,y) = (vx(y),0)

subject to the boundary conditions

vx−L∂yvx = 0, y = 0,H.

Substituting the ansatz for the velocity field into the Navier Stokes equations yield the ordinary differential equation

v′′x (y) =−
G
η
.

The two constants of integration are determined by the boundary conditions leading to the solution

vx(x,y) =
GH2

2η

(
y
H
−
( y

H

)2
+

L
H

)
, (1.22)

which is a modification of the well-known Hagen–Poiseuille flow. We note that the slip condition introduces a
simple additive modification to the no-slip solution, i.e.

vx(x,y) = vno-slip
x +

GHL
2η

.

6A typical approach to circumvent the problem numerically is to use the so-called numerical slip. The main observation is (see [RRL01]) that
in many cases an artificial slip is introduced by the discretization itself; see Section 7.5 for more details. Due to this numerical effect, the
contact line is able to move even though the no-slip condition is used. However, the numerical slip is typically strongly dependent on the
grid size. Moreover, the model which is supposed to describe the physics of dynamic wetting is in this case purely numerical and has no
physical foundation.

12


1.4. Dynamic Wetting

In particular, the limit L→ 0 is regular and Navier slip with a small slip length L on a molecular length scale is
practically indistinguishable from no-slip. For example, the average velocity throughout the cross-section of the
channel (which is proportional to the mass transport rate) is given by

〈vx〉=
GH2

12η

(
1+

6L
H

)
.

So, a slip length L = 10nm gives rise to a relative increase of the mass transport in a 1mm channel of 6 ·10−5 which
is irrelevant for any practical purposes. Remarkably, when we generalize this example to the case of a two-phase
flow in a narrow channel involving a moving contact line, the influence of the slip length even on a scale as small
as 10nm on the macroscopic mass transport rate is tremendous (at least in the framework of the “standard model”
of wetting). From a mathematical point of view, this is reflected in the fact that the moving contact line problem
with no-slip is ill-posed in the class of physically reasonable solutions with a finite rate of viscous dissipation. A
discussion of the influence of slip on the rise of liquid in a capillary is provided in Chapter 11.

Modeling the wettability: Besides the regularization of the no-slip paradox, the modeling of the wettability is
the second main ingredient for mathematical models of dynamic wetting. In the case of sharp interface models,
this is usually done (with some exceptions) by explicitly prescribing the contact angle at the solid boundary. The
experimentally measured contact angle, which is always subject to a finite measurement resolution, typically shows
a strong correlation to the capillary number [Hof75],

Ca =
ηVΓ

σ
,

where η and VΓ denote the dynamic viscosity and the contact line velocity, respectively. Motivated by this ob-
servation, many models prescribe the contact angle as some function of the capillary number and the equilibrium
contact angle, i.e.

θ = f (θeq,Ca), (1.23)

or vice versa

Ca = ψ(θ ;θeq). (1.24)

Since models of the above type are usually motivated by empirical correlations, we refer to (1.23) and (1.24)
as “empirical contact angle models”. Even though (1.23) and (1.24) are consistent closure relations for the en-
tropy production if the functions f and ψ satisfy an additional thermodynamic condition (see Theorem 3.6 in
Section 3.2), their validity is discussed controversially in the literature (see, e.g., [Shi20, Bot20c]). In particular,
it is argued that the dynamic contact angle cannot be prescribed a priori. Instead, experimental data suggest that
the dynamic contact angle depends, among others, on the details of the flow field in the vicinity of the contact
line [BBS99, Shi08].

Remark 1.5 (Dynamic contact angles and regularity). It can be shown by means of an asymptotic analysis for the
stationary Stokes equations that, for a fixed contact angle, the Navier slip condition with a finite slip length makes
the viscous dissipation rate finite, while the pressure is still logarithmically singular at the moving contact line (see,
e.g., [HM77], [Shi06]). This integrable type of singularity is commonly referred to as a weak singularity. It is an
interesting question, under which circumstances even the weak singularity is removed from the description. In the
publications [RE07, p.15] and [E11], Ren and E formulate the expectation that the weak singularity is removed
if, instead of a fixed contact angle, a certain model for the dynamic contact angle is applied. The present work
shows that this is not the case. Instead, it is shown that the dynamic behavior for sufficiently regular solutions to
the simplest model is unphysical (if the slip length is finite); see Chapter 4.

13


Part I.

Mathematical Modeling
&

Qualitative Analysis

15


2. A geometrical model for spreading drops in the
partial wetting regime

2.1. Introduction

(a) Glycerol-water droplet on a silicon wafer. (b) Water droplet on a PNIPAm polymer brush.

Figure 2.1.: Experimental data for the spreading of a droplet (experiments by B. Fickel and M. Hartmann
[FFH+20]).

The spreading of droplets on solid surfaces is a prototypical situation in the field of dynamic wetting and has
been studied extensively in the literature, see [dG85, BEI+09, dGBWQ04] and the references therein for a review.
As discussed in Chapter 1, the equilibrium state of a droplet on a surface depends strongly on the surface energy
of the material which is closely related to the surface properties, chemistry and microstructure. By modifying
the surface properties, one can therefore drastically change the wettability of the surface with a potentially large
impact on technical processes.

A particularly simple model of the spreading dynamics of a perfectly wetting drop is discussed in [dG85]. An
ODE model is proposed based on the observation that “in the regimes where gravitation is negligible, the macro-
scopic shape of the droplet is found to be rather close to a spherical cap” [dG85, p.851]. A geometrical relation is
used that links the base radius of a thin droplet of a given volume to the contact angle. Then one immediately ob-
tains a closed model for the evolution of the base radius if the contact line velocity and hence the capillary number
is a function of the contact angle as the only independent variable (or vice versa), i.e. if a functional relation of the
type (see page 13)

θ = f (Ca) or Ca = ψ(θ) (2.1)

holds for a given set of physical parameters. Since the geometrical relation assumes the contact angle to be small,
the model is usually applied only in the case of perfect wetting. The goal of the present chapter is to generalize the
geometrical model by de Gennes to the case of arbitrary contact angles.

Note that the condition (2.1) is discussed controversially in the literature and it is argued that the function f
(or ψ) should at least contain further independent variables. In particular, one can expect a much more complex
dynamics in the case of adaptive surfaces:

17


Chapter 2. A geometrical model for spreading drops in the partial wetting regime

“The challenge becomes even bigger for adaptive or responsive surfaces, that is, surfaces which
change their properties with external conditions such as temperature, humidity, electric and magnetic
fields or the presence of the liquid itself.” [BBS+18, p.11292]

In the present chapter1, we assess the ability of (2.1) to describe the spreading dynamics of a droplet in contact
with swellable polymer brushes. In particular, different amounts of ”pre-swelling” in humid air are considered.
The pre-swelling influences the imbibition of fluid into the polymer brush and hence the spreading dynamics
significantly. Note that a continuum mechanical model for the wetting on polymer brushes has been introduced
recently in [TH20]. The present chapter aims at an empirical mathematical model to describe the process of
spreading in terms of the base radius of the droplet as a function of time. Figure 2.1 shows experimental data for
the spreading of a glycerol-water droplet on a silicon wafer and the spreading of a pure water droplet on a PNIPAm
brush (see Section 2.3.2). Note that the timescale for the spreading dynamics on the polymer brush in Figure 2.1(b)
differs by orders of magnitude from the timescale for the water-glycerol droplet in Figure 2.1(a).

Assumptions: The reasoning which has been formulated for the special case of flat drops in [dG85] is based on
the following assumptions:

(I) The shape of the droplet throughout the spreading process is a spherical cap.

(II) The speed of the contact line can be expressed as a function of the contact angle2 θ (or vice versa).

(III) The volume of the drop is known (as a function of time if it is not conserved).

The ordinary differential equation (2.13) derived in Section 2.2 is a direct consequence of (I)-(III) without further
approximations or assumptions.

Remarks on assumption (I): The first assumption (I) greatly simplifies the description since in this case, the
shape is completely determined by only two parameters. Physically, the spherical cap is expected to be a good ap-
proximation in the late stage of spreading, i.e. for θ close enough to equilibrium if surface tension forces dominate
over both gravity and viscous forces, i.e. if the Bond and capillary numbers defined as

Bo =
ρgl2

σ
=

ρg
σ

(
3V
2π

)2/3

, Ca =
ηVΓ

σ
(2.2)

are sufficiently small3. However, in experiments, the shape of the drop is strongly influenced by the needle which
is used to deposit the drop on the substrate (see Figure 2.2 for the deposition of a pure water droplet). While
the contact line advances fast after the droplet gets in contact with the substrate (see Fig. 2.2 (b)-(c)), there is
a complex process of detachment from the needle which is accompanied by the propagation of capillary waves
and subsequent oscillations of the contact line (see Fig. 2.2 (d)-(f)). For partially wetting liquids (θeq > 0) with
low viscosity, the oscillations typically continue until the stationary radius is reached. Therefore, the assumption
(I) is typically not applicable in this case. However, it has been reported in the literature that for highly viscous
liquids such as water-glycerol mixtures, there is a viscous regime of spreading where viscous friction becomes the
main source opposing capillarity [CB14]. During this final stage of wetting after detachment from the needle, the
oscillations are mainly damped out by viscosity and the drop shape is quite close to a spherical cap. Moreover, the
spherical cap shape is also observed for the spreading on a polymer brush since in this case, the time scale for the
spreading is quite large (see Section 2.3.2).

1Please note: The results presented in the present chapter have been published as a preprint in [FFH+20].
2Here, by “contact angle” we mean the “macroscopic” contact angle that corresponds to the spherical cap shape.
3Note that the length scale for the Bond number is chosen to be the radius of a volume-equivalent drop with a contact angle of 90◦, i.e.

l := 3
√

3V/2π .

18


2.1. Introduction

(a) (b) (c)

(d) (e) (f)

Figure 2.2.: Deposition of a pure water droplet on a silicon wafer (Figure by M. Hartmann).

Remarks on assumption (II): Starting from the work by Hoffman [Hof75] and Jiang et al. [JSGS79], where it
is suggested that there is a universal relationship of the form

θapp = f (Ca), (2.3)

assumption (II) is a frequent choice to model dynamic wetting, even though its validity is discussed controversially
(see, e.g., [Shi08,Shi20,Bot20c]). Here θapp denotes the “apparent” (i.e. optically observable) contact angle and Ca
is the capillary number proportional to the contact line speed VΓ. Mathematically, it is more convenient to express
the capillary number as a function of the contact angle, i.e.

Ca = ψ(θapp), (2.4)

since there might be a whole interval of contact angles [θr,θa] where the contact line does not move. This well-
known effect is called contact angle hysteresis.

Example 2.1 (Cox-Voinov law). A prominent example of a relation of the type (2.4) is the Cox-Voinov law
(see [BEI+09] chapter III for a discussion) given by

G (θapp) = G (θeq)+Ca ln
( x

L

)
. (2.5)

Here x/L is the ratio of macroscopic to microscopic length scales, an unknown parameter that has to be found by
fitting experimental data (Bonn et al. [BEI+09] report x/L = 104). The function G defined as

G (θ) =
∫

θ

0

x− sinxcosx
2sinx

dx

can be well approximated by θ 3/9 for θ < 135◦ (see [BEI+09]), so that in the following we will refer to the
simplified equation

θ
3
app = θ

3
eq +9Ca ln

( x
L

)
(2.6)

as the “Cox-Voinov law”. Equation (2.6) results from an asymptotic solution of a hydrodynamic model by
Cox [Cox86], which is valid in the limit Ca,Re→ 0 for an ideal surface with no heterogeneities, if the microscopic
contact angle equals the equilibrium contact angle θeq. According to [BEI+09], the approximation θm ≈ θeq is
justified if viscous dissipation is large compared to local dissipation at the contact line.

Example 2.2 (Kistler’s empirical function). Another commonly used example is the empirical function by Kistler
[Kis93] given by the expression

θapp = fHoff(Ca)

= cos−1

(
1−2tanh

[
5.16

(
Ca

1+1.31Ca0.99

)0.706
])

.
(2.7)

19


Chapter 2. A geometrical model for spreading drops in the partial wetting regime

The latter function has been obtained from a fit to the experimental data by Hoffman [Hof75]. In order to describe
partial wetting, the function is shifted according to

θapp = fHoff[Ca+ f−1
Hoff(θ0)]. (2.8)

A direct comparison of the empirical models (2.6) for x/L = 104 and (2.8) in the case of complete wetting shows
that the two models are very similar; see Figure 2.3. Both empirical relations have been shown to be compatible
with experimental data over a range of capillary numbers for different fluid/substrate combinations.

Figure 2.3.: Kistler’s empirical function and the Cox-Voinov law for θeq = 0.

2.2. Spreading dynamics of spherical drops
Spreading dynamics in phase space: It can be advantageous to consider the spreading dynamics in phase
space4. The idea is to express the time-derivative of the base radius L̇ as a function of the base radius itself.
Formally, we define

L̇(L0) :=
dL
dt

(
L−1(L0)

)
. (2.9)

Note that this is possible since L(t) is a monotonically increasing function with inverse L−1. The advantage of
the latter approach is that it eliminates the dependence of the data on the arbitrary choice of t0. This is helpful
since the instant of contact t0 with L(t0) = 0 might, in practice, be hard to determine experimentally with sufficient
precision. Mathematically, the instant of contact marks a topological change of the flow which is followed by a
very fast initial dynamics, see [Bot20c] for a discussion.
In fact, the function L̇(L) eliminates the dependence on t0 since it is invariant with respect to a shift in the time
coordinate. It is easy to show that if L is given by

L(t) =C(t− t0)β

for some C, β > 0 and t0 ∈R, then
L̇(L) = (βC1/β )L1−1/β .

For the famous Tanner law [Tan79] (β = 0.1) this means L̇(L) ∝ L−9.

A Geometry-based model for spreading drops: For a spherical droplet, there is a purely geometric relationship
between the base-radius L, the volume V and the contact angle according to

L
V 1/3 = g(θ) := sinθ

(
π(1− cosθ)2(2+ cosθ)

3

)−1/3

; (2.10)

20


2.2. Spreading dynamics of spherical drops

(a) Notation. (b) The geometric relation g.

Figure 2.4.: Some basic geometry.

see Fig. 2.4(a) for notation. The function g is monotonically decreasing with θ (see Figure 2.4(b)) and hence
invertible on [0,∞). For small values of θ , i.e. for flat drops, the function g can be approximated by (see [dG85])

g(θ)≈ ĝ(θ) =
(

4
πθ

)1/3

.

Inverting relation (2.10) allows to express the contact angle as a function of the base radius and the volume accord-
ing to

θ = g−1
(

L
V 1/3

)
. (2.11)

Using the approximation for flat drops allows to simplify the above equation as

θ =
4
π

V
L3 . (2.12)

Assuming that the contact angle is related to the contact line speed VΓ by an empirical relation of the type (2.4)
allows to derive an ordinary differential equation for the spreading dynamics of spherical drops. Note that this ob-
servation has already been made by de Genne