Contents lists available atScienceDirect
SoftwareX
journal homepage:www.elsevier.com/locate/softx
Original software publication
WaveMaker: The three-dimensional wave simulation tool for falling
liquid films
Wilko Rohlfs
a,*
, Manuel Rietz
a, Benoit Scheid
baInstitute of Heat and Mass Transfer, RWTH Aachen University, Augustinerbach 6, 52056 Aachen, Germany bTIPs Laboratory, Université Libre de Bruxelles, C.P. 165/67, Avenue Franklin Roosevelt 50, 1050 Bruxelles, Belgium
a r t i c l e i n f o
Article history:
Received 15 May 2018
Received in revised form 13 July 2018 Accepted 13 July 2018
Keywords:
Falling films
Integral boundary layer model Spectral code
a b s t r a c t
WaveMaker is a MATLAB-based software for the simulation of spatially periodic hydrodynamic waves in
the two and three-dimensional domain. It utilizes the full or simplified second-order weighted residual integral boundary layer (WRIBL) model based on a pseudo-spectral scheme. A graphical user interface GUI is included for easy access to the simulation tool. Further, GPU’s for fast simulation are supported.
© 2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
Code metadata
Current code version 1.0
Permanent link to code/repository used of this code version https://github.com/ElsevierSoftwareX/SOFTX_2018_53
Legal Code License Apache License, Version 2.0
Code versioning system used none
Software code languages, tools, and services used Matlab 2015b - Matlab 2018a Compilation requirements, operating environments & dependencies none
If available Link to developer documentation/manual none
Support email for questions rohlfs@wsa.rwth-aachen.de
1. Motivation and significance
Falling liquid films are encountered in many industrial pro-cesses. Owing to their large contact areas, they offer small transfer resistances. Applications include phase separation, like for instance in sugar refinery or sea water desalination, and heat removal, e.g. in cooling towers of power plants or in cooling applications for electronic devices.
On the fundamental front, a falling liquid film can serve as a canonical system to study spatio-temporal chaos and weak/ dissipative turbulence [1]. Indeed, hydrodynamic waves arising at the surface of a falling film exhibit several spatio-temporal transitions from 2D traveling waves to 3D wave patterns, leading to solitary pulses and eventually interfacial turbulence [2]. This system has been studied for decades, since the seminal work of Kapitza [3], who identified the benefit of using low-dimensional
*
Corresponding author.E-mail address:rohlfs@wsa.rwth-aachen.de(W. Rohlfs).
models by averaging the balance equation over the thickness of the falling film, which is much smaller than the two other dimen-sions. There have been efforts to develop low-dimensional models. There are two main categories: (i) long-wave equations obtained by asymptotic expansion of the Navier–Stokes equations [4,5]; (ii) integral boundary layer equations obtained by averaging the boundary layer equations [6]. Nevertheless, these models suffer from bad quantitative comparisons with direct numerical simula-tions in the intermediate range of Reynolds numbers, say between 1 and 100. Improved models have then been proposed almost two decades ago, which combined an asymptotic expansion and an averaging procedure relying on a weighted residual method [7,8]. The model, referred to as the weighted residual integral boundary layer (WRIBL) method, exists in its full or simplified version, the latter being obtained by assuming the velocity profile to remain parabolic across the film thickness. First introduced for two-dimensional flow, it has been later expanded to three dimen-sions [9]. The excellent accuracy in predicting the dynamics of non-linear waves has been verified against direct numerical simulations
https://doi.org/10.1016/j.softx.2018.07.003
Fig. 1. Graphical user interface of WaveMaker.
Table 1
Different versions of the second-order WRIBL model that can be solved with
WaveMaker. The independent variables are indicated in parentheses. The
depen-dent variables are given for each case, together with the reference to the original equations.
WRIBL 2D 3D
(x,t) (x,z,t)
Simplified h,q h,q,p
Eqs. (11) and (41) in [7] Eqs. (6.2) in [9] withGx=1 Full h,q,s1,s2 h,q,p,s1,s2,r1,r2
Eqs. (11) and (38–40) in [7] Eqs (6.2a) and (C1) in [9]
(DNS) in various situations, as for instance for solitary waves in 2D [10,11] or 3D wave patterns [12]. As a result, the WRIBL model has become today a standard in simulating hydrodynamic waves at low computational cost.
This paper presents a software that solves the WRIBL model in two different versions as outlined inTable 1, along with the inde-pendent and deinde-pendent variables. The indeinde-pendent variables are the time t and the streamwise direction x for 2D, complemented by the spanwise direction z for 3D (seeFig. 3for an illustration of the different axes). Reminiscent to averaged equations, the cross-stream dimension is represented by a dependent variable, namely the film thickness h. The other dependent variables are the streamwise flow rate q, the spanwise flow rate p, the flow rate corrections s1and s2due to the departure from a parabolic velocity
profile in the streamwise direction and the corresponding flow rate corrections r1and r2in the spanwise direction.
2. Software description
2.1. Software principle
WaveMaker performs spatio-temporal simulations of the two and three-dimensional dynamics of falling liquid films using the full and simplified second-order models (seeTable 1). Periodic
boundary conditions in both x and z directions are imposed. This allows for the use of a pseudo-spectral scheme and thus benefits from the good convergence properties of spectral methods. The derivatives are evaluated in Fourier space and the nonlinearities in physical space. The time dependence is accounted for by a fifth-order Runge–Kutta method, which allows controlling the error and thus the adaptive time step from the difference with an embedded fourth-order scheme [13]. Details of the scheme, such as represen-tation of the variables in Fourier space, as well as the treatment of the associated ‘‘aliasing phenomenon’’ (see below), can be found in Appendix F of the book by Kalliadasis et al. [14].
2.2. Software architecture
The software architecture consists of a graphical user interface (GUI), and routines for the numerical simulation, and data post-processing. All routines are stored in the subfolder ./scr and orga-nized according to their functionalities. A brief description of each file is given in the file’s header. An image of the GUI is presented in
Fig. 1. The GUI of WaveMaker (Wavemaker.m) is designed to
•
define all relevant physical and numerical parameters as well as the output options, for 2D or 3D geometries;•
start the simulation using either the simplified or full second-order WRIBL model;•
monitor the wave profile and relevant physical parameters while computing.Each sub-panel of the GUI is described below. 2.2.1. Physical parameters
The relevant physical properties of the falling film flow can be defined in three different ways:
(i) In dimensional form using the kinematic viscosity
ν [
m2s−1]
,density
ρ [
kg m−3]
, surface tensionσ [
N m−1]
, magnitude|
g| [
m s−2]
and orientationθ [
◦]
and in the range of 0
<θ <
180) of gravitational acceleration, and Nusselt film thickness hN[
m]
. In addition, the domaindimensions in streamwise Lx
[
m]
and spanwise Lz[
m]
direc-tions need to be specified.
(ii) In dimensionless form using the Nusselt scaling with the Reynolds number Re, Kapitza number Ka, and inclination number Ct together with the streamwise kxand spanwise
kzwavenumbers, in units of h −1
N .
(iii) In dimensionless form using the Shkadov scaling with the reduced Reynolds number
δ
, reduced inclination numberζ
, and extensional viscous numberη
, which is associated with extensional viscous stresses. Note, that kxand kzare givenin units of (
κ
hN)−1, whereκ = η
−1/2.The Shkadov scaling accounts for the balance of gravity and viscous drag with inertia and surface tension, which is relevant for large-amplitude waves and gathers all second-order viscous effects in the boundary layer equations in front of the parameter
η
. The GUI also allows for a conversion between dimensional and dimensionless quantities. The definition of the dimensionless parameters can be displayed in the GUI by pressing the question-mark button next to the scaling.A non-trivial parameter is the dimensionless film thickness correction h0(in units of hN) which relates the Reynolds number
Re based on the Nusselt film thickness hNto the Reynolds number
based on the actual flow rate, denoted Req. The flow rate in a wavy
film is usually higher compared to a flat film. If simulations were conducted in an open domain with a defined volume flow at the inlet, the flow rate based Reynolds number had to be given, but in the present case of a periodic domain only the volume of liquid enclosed in this domain is prescribed. Thus, the correction factor can be used to obtain the desired flow rate based Reynolds number with periodic boundary conditions. For a flat film, this correction factor is unity and decreases with increasing waviness. A value of around h0
=
0.
9 has been shown [12] to give a good relationbetween the two Reynolds numbers in many cases (Re
=
40–60), particularly in the case of traveling waves. Traveling waves are stationary solutions in a frame of reference moving at the phase speed of the wave (see Section2.2.7). In such a case, Reqreaches aconstant value which can then be modified by adjusting the value of h0iteratively until Reqmatches the targeted Reynolds number
Re. For this purpose, the software tool allows to monitor the flow rate based Reynolds number Req(see again Section2.2.7).
To provide the user with some guidance on how to choose a rea-sonable domain size, the wavenumber kx,max,Nu(in Nusselt scaling)
corresponding to the most amplified wavelength in streamwise direction is displayed in the GUI. This wavenumber is derived from the linear dispersion relation of the full second-order model [15]. 2.2.2. Numerical parameters
The numerical parameters are the anti-aliasing filter, the streamwise and spanwise resolution, and the tolerance on the variables.
Nonlinearities are well known to generate aliasing errors. These are distortions of high frequency Fourier modes due to the trun-cation in Fourier space (sampling) when pseudo-spectral meth-ods are used. The anti-aliasing filter removes the high frequency modes in the solution process. A standard value for the aliasing filter is
α =
0.
66 according to the two-thirds rule that applies for quadratically nonlinear equations (as for the Navier–Stokes equations) [16]. Despite the fact that the order of nonlinearity of the WRIBL equation system is five (much larger than two) it has been shown that keeping the first two-thirds of the Fourier modes in each direction before each iteration is in most cases sufficient to reach convergence.The spatial resolution, M streamwise and N spanwise grid points, defines the number of modes in the frequency spectrum
of the spectral code, which are M
/
2×
N/
2, orα
M/
2×
α
N/
2 with the aliasing treatment. Therefore, the computational domain of size Lx×
Lz is discretized with regularly spaced grid pointswith coordinates xi
=
iLx/
M and zj=
jLz/
N. For the simulationof two-dimensional waves, the dropdown menu for the spanwise resolution can be set to 2D (seeFig. 1) corresponding to only two spanwise grid points being used.
The relative tolerance on the variables is used as an input for the fifth-order Runge–Kutta adaptive time-step routine ode45 in Matlab. The relative tolerance for each variable has been fixed to 10−4, and is thus not accessible in the GUI.
2.2.3. Initial condition
The initial condition can be defined in two different ways. (1) By clicking on ‘Define IC’ and initiating a film profile using a superposition of sine and cosine functions with defined amplitudes in streamwise and spanwise directions. By default, the wavelength equals the domain size. However, it is also possible to generate multiple waves in the domain in order to investigate wave-to-wave interactions. In addition, either two-dimensional noise or spanwise noise can be added. Amplitudes of the trigonometric functions and noise are defined in units of the mean film thickness. (2) By loading an existing solution.
2.2.4. Output parameters
Similar to the different options for the dimensional or dimen-sionless input parameters, the user can define the scaling for the output. Additionally, the final time (tend) of the simulation and the
print interval (tprint) can be specified in seconds (‘Dimensional’)
or in dimensionless form in units of h2
N
/
(3ν
Re) (‘Nusselt scaling’)or in units of hN
/
(3ν
Reη
1/2) (‘Shkadov scaling’). Output optionsinclude Matlab-based figures of the simulated wave profile (*.fig) and associated raw data (*.mat). Additionally, the velocity and the wave shape can be exported in the .vtk format to be post-processed using freeware tools like ParaView. By default, the solution for every print interval is stored and can be loaded as initial condition. 2.2.5. Load/save parameters
The entire set of parameters can be saved and imported. Param-eter files are stored in the subfolder ./paramParam-eterfiles.
2.2.6. Simulation
In the simulation panel, the user can select the model type which is either the simplified model or the full second-order model. The simulations can be started, stopped and continued with the respective buttons. If a graphical processor unit (GPU) exists in the machine, this option can be activated by the check box. 2.2.7. Results
The results screen on the right side allows to monitor the wave profile according to the parameter Tprint. In the 2D case, the
wave profile is shown together with a plot of the streamlines. If WaveMaker detects a traveling wave the plot is given in the moving frame of reference together with values of the wave velocity and the wave frequency. If the wave is non-stationary, the streamlines are given in the laboratory frame of reference. In the 3D case, the surface topology of the film is shown. For both, the 2D and the 3D cases, the time evolution of the maximum/minimum film thickness and the flow rate based Reynolds number are displayed in the bottom plot.
3. Application examples
Multiple test cases taken from literature are provided for ex-amination in the subfolder ./parameterfiles. The different cases are listed inTable 2, including respective dimensionless parameters in Nusselt scaling and the film-thickness correction h0.
Table 2
Example cases, which can be loaded from the WaveMaker GUI. For consistency, flow rate, dimensionless fluid properties, inclination, and the domain dimensions are given in Nusselt scaling. The parameter files in WaveMaker also provide dimensional values. Name Re Ka Ct kx kz h0 Reference Case I 40.8 3923 0 0.0729 0.0729 0.9147 [17] Case II 59.3 3923 0 0.0661 0.0826 0.9033 [17] Case III 15 509.5 0 0.101 – 0.8922 [18] Case IV 6.8 509.5 0 0.078 – 0.8922 [19] Case V 6.2 17.8 0 0.263 0.263 0.9012 [12] Case VI 20.1 515 0 0.25 – 0.91 [20] Case VII 69 515 0 0.125 – 0.91 [20] Case VIII 1.05 5.1 −0.4 0.0193 0.0966 0.91 [21] Case IX 1.05 5.1 −1 0.0212 0.106 0.91 [21]
Fig. 2. Comparison of experimental [19] and numerical [11] results for cases III and IV outlined inTable 2, i.e. for Re=6.8 (left) and for Re=15 (right). The upper plots illustrate the film thickness profile. The lower plots show the streamwise velocity given at a distance of 0.1 mm from the wall.
3.1. Simulation of two-dimensional surface waves
Fig. 2provides a direct comparison between experimental data from Dietze et al. [18], fully resolved numerical simulations from a previous study [11], and two-dimensional simulations using the full second-order WRIBL model implemented in WaveMaker. The data is taken from a recent study by Rohlfs et al. [11]. The first example case (IV) is characterized by a nearly-sinusoidal wave with no recirculation in the main hump (Re
=
6.
8 and f=
24 Hz, seeFig. 2, left). The second exemplary two-dimensional case (III) is characterized by a solitary-like wave with the presence of a recirculation zone, evident through a negative velocity in front of the main wave hump (Re=
15 and f=
16 Hz, seeFig. 2, right). Very good agreement is obtained in the wave shape and the velocity profile between the fully resolved simulations, and the full second-order WRIBL model.3.2. Simulation of three-dimensional surface waves
The benefits of using low-dimensional models with WaveMaker for the simulation of falling films become especially evident when three-dimensional films are investigated. In this context, simula-tions using the WRIBL model can be conducted with significantly reduced computational power compared to fully resolved 3D nu-merical simulations. Moreover, large domain simulations, e.g. in-teracting wave patterns without pronounced symmetries, become manageable without an extensive computational architecture.
Information about the difference in computational effort be-tween fully resolved DNS and simulations using the WRIBL model is available for Case I of the supplied example cases. In Dietze et al. [12], respective DNS have been performed on the JUROPA supercomputer of Forschungszentrum Jülich. The computational
effort for the simulation of one second of film evolution amounted to about 80000 processor hours. Simulations had been performed on 1024 processors.
Respective simulations of Case I using the WRIBL model with 512 x 512 grid points (convergence was shown for this resolution in [12]) require a computational effort of about 1500 processor hours on a personal computer (Intel Core i5-4570 CPU @ 3.20 GHz), which is less than 2% of the time required for DNS.
Fig. 3shows WaveMaker simulation results for Case I ofTable 2
in a relatively small domain, comprising one wavelength in the streamwise direction and two wavelengths in the spanwise direc-tion. The film topology is displayed after a time period of 0.175 s, namely shortly after the destabilization of quasi 2D wave fronts into regular 3D surface waves with capillary ripples preceding the main wave hump. A 2D wave profile along a symmetry plane of the film topology is presented. The surface topology compares very well with experimental data of Park et al. [17]. A full discussion and further comparisons between experimental data and WRIBL model results is presented in the work of Dietze et al. [12].
InFig. 4, large domain WaveMaker simulations (0
.
4×
0.
08 m2)of a film on the outside of a vertical rotating cylinder of large radius are presented (cases VIII and IX) together with their respec-tive experimental data [21]. In this context, negligible curvature of the cylinder justifies comparability to WRIBL simulations. The different body forces acting on the fluid have been accounted for in WaveMaker by assuming the flow to be subject to gravitational acceleration with its components in streamwise and cross-stream directions equaling standard gravitational acceleration and cen-trifugal acceleration, respectively. A related problem would be a falling film hanging beneath an inclined wall, where the wall normal component of gravitational acceleration acts as the desta-bilizing body force acting on the system following the Rayleigh– Taylor instability mechanism [15]. In both cases, the long-term
Fig. 3. Surface topology of a falling film calculated with WaveMaker (Case I, t=0.175 s). The simulation has been conducted on a personal computer (Intel Core i5-4570 CPU @ 3.20 GHz) within a few minutes.
Fig. 4. Experimental results (a) and WaveMaker simulations (b) of a falling film on the outside of a rotating vertical cylinder of large radius showing rivulet formation for
different values of the parameter G, which is the quotient of centrifugal acceleration and gravitational acceleration, t being the dimensionless timescale in Shkadov scaling andλRiv,iis the spanwise wavelength of rivulets.
nonlinear film evolution displays rivulet formation with distinct wavelength selection, visible both in experiments (Fig. 4left) and WaveMaker simulations (Fig. 4 right). In this context, the large simulation domain on the one hand minimizes influences of do-main size on wavelength selection and on the other hand enables investigation of larger scale rivulet and wave interactions. Note,
that experimental and WaveMaker results inFig. 4are shown for different rotation speeds of the cylinder, i.e. different values of the destabilizing body force. With WaveMaker, rivulet formation could be investigated at low rotational speeds, meaning rivulet inception lengths larger than the length of the experimental setup (
≈
0.
8 m).4. Software limitations and expandability
4.1. Software limitations
The nature of the low dimensional WRIBL model implemented in WaveMaker involves some inherent limitations of applicability. WaveMaker allows for the simulation of isothermal falling films on a plane wall. Furthermore, as the gaseous phase over the film is not modeled, no respective interactions are considered. Being a single value model, dripping from falling films, expected for hanging or centrifuged films as presented in Fig. 4, cannot be resolved. Hence, inaccuracies are expected close to the dripping limit. However, validity of simulation results in the linear and early nonlinear regime of film evolution has been demonstrated for cen-trifuged films [21] and improvement of the WRIBL model including the full-curvature has been proposed near the dripping limit [22].
The approach used in the WRIBL model assumes a parabolic ve-locity profile across the film. Although correction terms accounting for a deviation from that assumption are part of the full second-order model, inaccuracies are expected for higher values of the Reynolds number. A good agreement of model results with the experimentally observed evolution of waves in the linear and non-linear regime has been demonstrated up to a Reynolds number of approximately 100 [8,12] even though this limit strongly depends on all other parameters.
4.2. Expandability
Despite the inherent limitations, the WRIBL model imple-mented in WaveMaker has proven to be a powerful tool for the investigation of film flows and was widely used in recent years. Adding to this, many efforts have been made to increase the range of applicability by extending the modeling approach to more com-plex situations. For instance, non-isothermal falling films [23,24], falling films subject to a counter-current gas flow [25,26], falling films undergoing phase change [27] or dielectric films in an electric field [28]. An upgrade of WaveMaker by implementing further variations and extensions of the integral boundary layer model is intended.
5. Impact
WaveMaker is a software able to compute three-dimensional wave patterns in falling films with unprecedented rapidity as compared to direct numerical simulations. This is allowed thanks to the combination of three factors: (i) the use of validated low-dimensional WRIBL models; (ii) the use of a spectral method with the FFTW (the Fastest Fourier Transform in the West) algo-rithm [29] reducing the computing time; (iii) the implementation of GPU-accelerated computing. It is thus believed that the software can be used by researchers as well as engineers to quickly study the behavior of a falling film under various conditions, including films hanging on the underside of an inclined wall. The software also enables simulating large areas, barely accessible by DNS because of prohibitive computing time. This advantage allows exploring complex wave-interaction behaviors such as wave coalescence, wave locking or interfacial weak turbulence. Additionally, the ex-plicit scheme makes the implementation of other long-wave type equations easy, opening new possibilities for further investigations of thin film dynamics.
Acknowledgments
The authors are grateful to Pierre Colinet who provided the backbone routines based on Numerical Recipes for building the
present code and Pierre-Alexandre Petitjean for debugging a first version of WaveMaker. B.S. and W.R. acknowledge the financial support of EU-FP7 ITN Multiflow (No. 214919) as well as the BELSPO agency under the Grant No. IAP-7/38 MicroMAST. B.S. also thanks the F.R.S.-FNRS (No. 197890) for financial support. This research has been performed under the umbrella of the COST Action No. MP1106.
References
[1] Denner F, Charogiannis A, Pradas M, Markides CN, van Wachem BGM, Kalliada-sis S. Solitary waves on falling liquid films in the inertia-dominated regime. J Fluid Mech 2018;837:491–519.
[2] Pradas M, Kalliadasis S, Nguyen PK, Bontozoglou V. Bound-state formation in interfacial turbulence: direct numerical simulations and theory. J Fluid Mech Rapids 2013;716:R2.
[3] Kapitza PL. Wave flow of thin layers of a viscous fluid: I. Free flow - II. Fluid flow in the presence of continuous gas flow and heat transfer. In: Haar DT, editor. Collected Papers of P. L. Kapitza (1965). Pergamon (Oxford); 1948. p. 662–89. (Original paper in Russian: Zh. Eksp. Teor. Fiz. 18, I. 3–18, II. 19–28). [4] Benney DJ. Long waves on liquid films. J Math Phys 1966;45:150–5. [5] Ooshida T. Surface equation of falling film flows with moderate reynolds
number and large but finite weber number. Phys Fluids 1999;11:3247–69. [6] Shkadov VYa. Wave flow regimes of a thin layer of viscous fluid subject to
grav-ity. Izv Akad Nauk SSSR, Mekh Zhidk Gaza 1967;1:43–51. English translation in Fluid Dynamics 2, 29–34, 1970 (Faraday Press, N.Y.).
[7] Ruyer-Quil C, Manneville P. Improved modeling of flows down inclined planes. Eur Phys J B 2000;15(2):357–69.
[8] Ruyer-Quil C, Manneville P. Further accuracy and convergence results on the modeling of flows down inclined planes by weighted-residual approxima-tions. Phys Fluids 2002;14(1):170–83.
[9] Scheid B, Ruyer-Quil C, Manneville P. Wave patterns in film flows: modelling and three-dimensional waves. J Fluid Mech 2006;562:183–222.
[10]Chakraborty S, Nguyen PK, Ruyer-Quil C, Bontozoglou V. Extreme solitary waves on falling liquid films. J Fluid Mech 2014;745:564–91.
[11]Rohlfs W, Pischke P, Scheid B. Hydrodynamic waves in films flowing under an inclined plane. Phys Rev Fluids 2017;2:044003.
[12]Dietze GF, Rohlfs W, Nährich K, Kneer R, Scheid B. Three-dimensional flow structures in laminar falling films. J Fluid Mech 2014;743:75–123.
[13]Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in C - The art of scientific computing. 2nd ed. New-York: Cambridge University Press; 1992.
[14]Kalliadasis S, Ruyer-Quil C, Scheid B, Velarde M. Falling liquid films. London: Springer-Verlag; 2012. p. 440.
[15]Scheid B, Kofman N, Rohlfs W. Critical inclination for absolute/convective instability transition in inverted falling films. Phys Fluids 2016;28:044107. [16]Boyd JP. Chebyshev and fourier spectral methods. Dover publication; 2001. [17]Park CD, Nosoko T. Three-dimensional wave dynamics on a falling film and
associated mass transfer. AIChE J 2003;49:2715–27.
[18]Dietze GF, Al-Sibai F, Kneer R. Experimental study of flow separation in laminar falling liquid films. J Fluid Mech 2009;637:73–104.
[19]Dietze GF. Flow separation in falling liquid films [Ph.D. thesis], RWTH Aachen University; 2010.
[20]Gao D, Morley NB, Dhir V. Numerical simulation of wavy falling film flow using VOF method. J Comput Phys 2003;192:624–42.
[21]Rietz M, Scheid B, Gallaire F, Kofman N, Kneer R, Rohlfs W. Dynamics of falling films on the outside of a vertical rotating cylinder: waves, rivulets and dripping transitions. J Fluid Mech 2017;832:189–211.
[22]Kofman N, Rohlfs W, Gallaire F, Scheid B, Ruyer-Quil C. Prediction of two-dimensional dripping onset of a liquid film under an inclined plane. Int J Multiph Flow 2018;104:286–93.
[23]Ruyer-Quil C, Scheid B, Kalliadasis S, Velarde M. Thermocapillary long waves in a liquid film flow. Part 1. Low-dimensional formulation. J Fluid 2005;538:199– 222.
[24]Trevelyan P, Scheid B, Ruyer-Quil C, Kalliadasis S. Heated falling films. J Fluid Mech 2007;592:295–334.
[25]Tseluiko D, Kalliadasis S. Nonlinear waves in counter-current gasliquid film flow. J Fluid Mech 2011;673:19–59.
[26]Dietze GF, Ruyer-Quil C. Wavy liquid films in interaction with a confined laminar gas flow. J Fluid Mech 2013;722:348–93.
[27]Pillai D, Narayanan R. Rayleigh-Taylor stability in an evaporating binary mix-ture. J Fluid Mech 2018;848.
[28]Wray AW, Matar OK, Papageorgiou DT. Accurate low-order modeling of electrified falling films at moderate Reynolds number. Phys Rev Fluids 2017;2:063701.
[29] Frigo M, Johnson SG. (1998) FFTW: an adaptive software architecture for the FFT. In: Proceedings of the 1998 IEEE international conference on acoustics, speech and signal processing, vol. 3, Seattle, WA. p. 1381–4.
