Active control of liquid film flows: beyond reduced-order models
AActive control of liquid film flows: beyond reduced-order models
Radu Cimpeanu , , Susana N. Gomes Demetrios T. Papageorgiou Mathematics Institute, University of Warwick, Coventry, CV4 7AL, UK Department of Mathematics, Imperial College London, London, SW7 2AZ, UK Mathematical Institute, Andrew Wiles Building, Radcliffe Observatory Quarter,Woodstock Road, Oxford, OX2 6GG, UKAugust 31, 2020
Abstract
The ability to robustly and efficiently control the dynamics of nonlinear systems lies at theheart of many current technological challenges, ranging from drug delivery systems to ensuringflight safety. Most such scenarios are too complex to tackle directly and reduced-order modellingis used in order to create viable representations of the target systems. The simplified settingallows for the development of rigorous control theoretical approaches, but the propagation oftheir effects back up the hierarchy and into real-world systems remains a significant challenge.Using the canonical setup of a liquid film falling down an inclined plane under the action ofactive feedback controls in the form of blowing and suction, we develop a multi-level modellingframework containing both analytical models and direct numerical simulations acting as an insilico experimental platform. Constructing strategies at the inexpensive lower levels in the hier-archy, we find that offline control transfer is not viable, however analytically-informed feedbackstrategies show excellent potential, even far beyond the anticipated range of applicability of themodels. The detailed effects of the controls in terms of stability and treatment of nonlinearityare examined in detail in order to gain understanding of the information transfer inside the flows,which can aid transition towards other control-rich frameworks and applications.
The ability to manipulate physical systems has applications in all areas of engineering, medicine, andmany other fields. Examples range from path tracking and planning for self-driving cars [36, 50], in-flight measurement and manoeuvre adjustments in unmanned aerial vehicles (e.g. drones) [12, 29],management of crowds during mass gatherings or evacuations [9, 19], to controlled-release deliverysystems for targeted drug delivery [44], to name but a few. However, before being able to controlreal-life systems, one must be able to model them accurately. There have been a number of recentinvestigative efforts on data-driven modelling and control of real-life systems, such as the derivationof coarse-grained PDEs from macroscopic observations [28], the use of machine learning to predictbehaviour and learn models in fluid mechanics [11], or the design of controls based on modelswith uncertainty [25], to name but a few. However, these settings do not allow for the rigoroustheoretical development of control strategies and for these reasons we chose the reduced-order modelapproach that follows. Real-life phenomena usually lead to extremely intricate models which aredifficult to tackle analytically and are challenging computationally; this motivates the developmentof physically relevant reduced-order modelling approaches resulting in systems that are amenable tomathematical and computational analysis. Typical examples include fluid dynamics [27] and image1 a r X i v : . [ phy s i c s . f l u - dyn ] A ug estoration [3], see [5] for a recent review. The resulting models are usually accurate in certainregions of parameter space and are used as surrogates of the original system, thus enabling detailedexploration of underlying phenomena without recourse to physical experiments and direct numericalsimulations ( in silico or virtual experiments). Most importantly, such reduced-order models can beused as the starting point in controlling the dynamics of the more general physical model, and thismulti-way synergy between hierarchical models constitutes our present methodology.This study has two main goals: (i) Establish the advantage of using control methodologiesdesigned for reduced-order models in an in silico experiment; (ii) Use this framework as a test-bedfor informing the control of large scale systems. Our approach of using controls based on modelslower in a hierarchy and moving up to the full physical problem, circumvents challenges present inthe latter class including lack of rigorous analytical results and reliance on expensive trial and errornumerical experiments. Proceeding directly with the full model could predict qualitative featuresin certain situations, but to our knowledge this is the first attempt at propagating quantitativeinformation between vastly different complexity levels.The physical model we select to study is that of a thin film of water flowing down an inclinedplane. This has been studied extensively both theoretically and experimentally and is ideal for ourgoals since accurate and efficient direct numerical simulations (DNS) of the Navier–Stokes equationscan be utilized as in silico experiments. Furthermore, falling film flows have benefited from a strongbody of analytical efforts for reduced-order modelling over the last 50 years (see, for example, [6, 39];we also point out that the simplest model, the Kuramoto-Sivashinsky equation is also widely usedas a model for flame front propagation, reaction-diffusion systems and other relevant scenarios[26, 38, 37]), as well as very recent comparisons between reduced-order models, full computationsand experiments in relevant regimes [14].Invoking the multiscale nature of the problem (the film is thin compared to its wavelength)leads to reduced-order long-wave or weakly nonlinear models [20]. These have been generalised toincorporate various effects that can act as active control mechanisms, such as electric fields [2, 45, 48],wall heating [7, 41] and same-fluid blowing and suction [43]. Alternative and/or supplementarypassive control can be achieved, for example, by including substrate topography [15, 47, 35] or theinclusion of surfactants [8, 10, 16, 22, 21]. In this paper, we focus on same-fluid blowing and suction,which has recently been used to develop efficient feedback controls based on interface observationsfor a hierarchy of models of increasing complexity [4, 18, 46, 42].Our results show that controls learned from reduced-order models based on observations of theinterface, are not directly applicable to the virtual experiment if they are applied pre-experiment inthe real-world system. Nevertheless we can still obtain crucial information from the reduced-ordermodels to be used in the full problem. The actuation rules (which are based on properties of thereduced order models, often their linear stability properties) remain the same, only now the controlsare predicated on observations of the interfacial shape in the virtual experiment. Furthermore, thepredictions on the strength and number of controls needed to stabilise any desired solution, as wellas their location and shift-displacement from observations is correct as demonstrated by extensivenumerical simulations. The virtual experiments also allow us to explore the control methodologybeyond the parameter regimes where the reduced order models are valid, hence opening a new avenuefor control development learned from faster and parametrically rich computations. We emphasisethat our results and methodology are not restricted to the particular physical model consideredhere, and are expected to be applicable to other similar systems, such as self-driving cars, smartrobots, or active matter. 2 The control methodology
In this section, we outline the theory behind feedback control, which is the basis for the design ofthe controls developed in [18, 42]. We then present a hierarchy of models for the dynamics of theinterface of a thin film of water flowing down a plane in two dimensions (2D), which is the physicalsystem of interest – this hierarchy is obtained using asymptotic expansions that lead to reducedorder models.
The theory of linear feedback control was first established for (systems of) ordinary differentialequations (ODEs) in [49]. To introduce this methodology we consider the simplest example ofcontrolling a scalar ODE ˙ y = λy, y (0) = y , (1)where the dot represents derivative with respect to time. It is well known that the solution of (1)is y ( t ) = y e λt and that y ( t ) → λ < y ( t ) → ∞ if λ is positive.The main goal of linear control theory is to introduce a control to equation (1), i.e.˙ y = λy + f, y (0) = y , (2)and choose f in such a way that the solution is stabilised , i.e. , driven to y = 0 even for λ > proportional feedback control would be to choose f ( t ) = − αy ( t ) for somepositive constant α . It is clear that if α is such that λ − α <
0, then y ( t ) → t → ∞ and wesay that the control f ( t ) stabilises the solution to the ODE. The term feedback is used because thecontrol uses information on the current state of the system; the control is called proportional sinceit is proportional to the current solution. A similar idea can be used in systems of ODEs˙ y = A y + f , y (0) = y , (3)where now y , y , f ∈ R d and A is a d × d matrix. Similarly to the one dimensional case, when f = ,if the eigenvalues of A all have negative real part, the solution is asymptotically stable , i.e., y ( t ) → as t → ∞ . The analogue of the previous control here is to use f ( t ) = − α y ( t ) = − αI y ( t ), where I is the d × d identity matrix. A simple calculation can be used to find the smallest α necessary tostabilise the system, namely, choose α such that the eigenvalues of A − αI all have negative realparts.It is naturally more efficient to choose the controls using more information about the model.This can be done by designing the control as follows˙ y = A y + B f , y (0) = y . (4)Here, B is a d × M matrix that encodes some information about how one applies the controls (e.g.,if the controls are localised in some part of the domain), and the controls are now f ∈ R M . Notethat we can have M = d and B = I , which is the case outlined above. It can be shown that undercertain assumptions on the matrices A and B (namely, the Kalman rank condition [49]), one canfind a matrix K such that the eigenvalues of A + BK all have negative real part, and therefore thecontrols f = K y stabilise the system. The matrix K can be computed using, e.g., a pole placementalgorithm [23] or by solving a linear-quadratic regulator problem [49].In more realistic physical systems, the dynamics of the solution are modelled by a nonlinearsystem of ODEs, ˙ y = F ( y ) + B f , y (0) = y , (5)3igure 1: Schematic of a falling liquid film down an inclined plane allowing for blowing and suctioncontrols. The dynamics of the interface at y = h ( x, t ) is controlled not only by fluid parameters,but also by the geometry inclination angle θ and the imposed control values v = F ( x, t ) at y = 0 ina suitably aligned coordinate system.where F is some nonlinear function of y . In addition, physical continuum models provide (linear ornonlinear) partial differential equations (PDEs), e.g., a reaction-diffusion equation for the evolutionof a population, tumour growth or other biological and chemical applications. Such PDEs take thegeneral form u t = L u + N ( u ) + f, (6)along with appropriate initial and boundary conditions. The subscript t denotes time derivative,and L , N are linear and nonlinear differential operators, respectively. By projecting this equationto an appropriate basis (e.g., taking Fourier transforms), one can write the PDE as an infinite-dimensional system of ODEs such as (5). Passing to this limit is not straightforward, even for linearPDEs [49], however, in certain cases this is possible, even for nonlinear PDEs, as was shown for theso-called Kuramoto–Sivashinsky equation in [4, 18]. In this case, in order to prove that the controlsstabilise the full nonlinear system, one needs to use nonlinear stability analysis techniques, such asfinding a Lyapunov function [24]. This section is devoted to a brief description of the physical problem, its mathematical modellingand the asymptotic work that leads to the hierarchy of reduced-order models we are focussing on.Consider a thin film of water flowing down a 2D plane inclined at an angle θ to the horizontal, seeFigure 1. Throughout this discussion, we will be using the ( · ) ∗ notation to distinguish dimensionalquantities from their undecorated dimensionless counterparts. The film thickness is denoted by h ∗ ( x, t ) and our main goal is to drive the system to its undisturbed flat interface, h ∗ ( x, t ) = h ∗ .We will apply proportional feedback controls which are actuated by means of blowing or suctionof water via slots in the wall, and the controls will be designed based on readings of the interfacialheight.The full mathematical model for the fluid motion is given by the Navier–Stokes equations forboth liquid film and gas above it with a suitable nonlinear coupling at the interface. Due to thepassive nature of the gas in our applications (in view of large density/viscosity ratios), one canrestrict attention to the liquid film alone. The blowing/suction conditions at the wall y = 0 are u ∗ = 0 , v ∗ = F ∗ ( x ∗ , t ∗ ) , (7)4here u ∗ , v ∗ are the streamwise (parallel to the inclined plane) and transverse (perpendicular to theplane) velocities, respectively, as depicted in Figure 1. The uncontrolled system admits a uniformflat film solution known as the Nusselt solution [20], given by h ∗ ( x, t ) = h ∗ and a semi-parabolicin y ∗ horizontal fluid velocity with surface value U ∗ s = ρ ∗ g ∗ sin θ ( h ∗ ) µ ∗ , where g ∗ , ρ ∗ and µ ∗ are theacceleration of gravity and the fluid’s density and viscosity, respectively. An appropriate non-dimensionalisation of this problem allows us to define two important dimensionless parameters thatcharacterise the system. The Reynolds number Re = ρ ∗ U ∗ s h ∗ µ ∗ , and the capillary number Ca = µ ∗ U ∗ s γ ∗ ,which measure the relative importance between inertia and viscosity, and between gravity andsurface tension (represented by γ ∗ ), respectively.As explained previously, full models such as the Navier–Stokes equations are computationallyexpensive to simulate. However, in the case of thin liquid films, the mean interface height h ∗ ismuch smaller than the length of the domain, L ∗ = nh ∗ , so that we can define a long wave parameter (cid:15) = h ∗ L ∗ = n − (cid:28)
1. This disparity of scales facilitates a multiscale approach to derive from firstprinciples hierarchies of reduced-order models. For the remainder of the modelling discussion weuse h ∗ and U ∗ s as reference length and velocity scales alongside the defined groupings to transfer thesystem to its dimensionless counterpart (and drop the decorations accordingly). In this context,the requisite assumptions are: (A1) (long-wave assumption) the geometrical aspect ratio (cid:15) is small; (A2) The Reynolds number Re is O (1); (A3) Surface tension is sufficiently strong to appear at leading order, i.e., the capillary number issmall, and Ca = O ( (cid:15) ) is the appropriate distinguished limit; (A4) The controls F are small F = O ( (cid:15) ), implying weak blowing/suction.Using assumptions (A1) - (A4) and asymptotic analysis techniques, Thompson et al. [42] derivedtwo different long-wave models. Both models satisfy a mass conservation equation h t + q x = F ( x, t ) , (8)and couple with an equation for the flux q ( x, t ). In the first model, the Benney equation, theyobtain an explicit expression for q ( x, t ) and the model is a single PDE for the interfacial height h ( x, t ): q ( x, t ) = h (cid:18) − h x cot θ + h xxx Ca (cid:19) + Re (cid:18) h h x − h F (cid:19) . (9)The second model is the weighted residuals model, which describes the evolution of the interfacialheight h ( x, t ) and the flux q ( x, t ):2 Re h q t + q = h (cid:18) − h x cot θ + h xxx Ca (cid:19) + Re (cid:18) q h x − hqq x
35 + hqF (cid:19) . (10)We note that the controls appear as an inhomogeneous term F ( x, t ) in the mass conservationequation (8), and this structure plays a crucial role in the efficiency of these controls.Due to the asymptotic reduction, the models do not directly provide the evolution of bulkquantities such as the streamwise and transverse velocities u and v , but these are known in termsof the interfacial height h ( x, t ) and the flux q ( x, t ) from the analysis. At leading order, and usingthe results from the weighted residuals model, for example, the fluid velocities are u ( x, y, t ) = 3 qh (cid:18) yh ( x, t ) − y h ( x, t ) (cid:19) , v ( x, y, t ) = F ( x, t ) − (cid:90) y u ( x, y (cid:48) , t ) dy (cid:48) , (11)5igure 2: Validation summary for a film of undisturbed thickness of h ∗ = 175 µ m on an inclinedplane described by θ = π/
3. The top panel, shows, from left to right, the interfacial film height asobtained from calculations based on the Kuramoto–Sivanshinsky equation, the weighted residualslong-wave model, as well as direct numerical simulations. Dashed lines indicate a virtual measure-ment station for the interfacial height showing the evolution of the film in the bottom left panel.The bottom right panel compares saturated interfacial shapes for the long-wave model and theDNS.thus providing a description of the flow field in the whole domain that can be compared withdirect numerical simulations, for instance. Details on the numerical methodology behind solvingthese reduced-order models are provided in Appendix B. We underline however that we have useda unimodal perturbation of sufficiently small amplitude, typically of O (10 − ), as initial interfaceshapes in all our numerical solutions, including higher up in the hierarchy at the DNS level to ensureconsistency in comparison.The above long-wave models are significantly more accessible computationally than the fullNavier–Stokes equations, but they are still highly nonlinear and to-date have not been tackledanalytically to obtain rigorous results. Hence, further simplifications are necessary in order to makeanalytical progress. One can, for example, perform weakly nonlinear analysis to obtain a Kuramoto–Sivashinsky (KS) equation for small but nonlinear perturbations from a flat interface [20, 42]. TheKS equation is a fourth order nonlinear PDE having the same form as (6) with L u = u xxxx + u xx and N ( u ) = uu x . The KS equation appears in a plethora of applications and is widely studiedsince it is one of the simplest model PDEs which exhibit spatiotemporal chaotic behaviour. Overthe last few decades existence and uniqueness of solutions has been explored [40], different typesof attractors have been characterized [13], and the route to chaos for solutions of the KS equationhave been reported [30], thus exemplifying the range of interesting analytical and computationalresults that can be achieved even at this lowest member of the model hierarchy.Some uncontrolled computations are presented next to showcase the underlying nonlinear dy-6amics. Figure 2 compares the results for the evolution of the interfacial height using the KSequation (top-left panel), our most comprehensive long-wave model (the weighted residuals model,top-middle panel) and the virtual experiment solving the Navier–Stokes equations (top-right panel)for a film thickness of h ∗ = 175 µ m and an inclination angle of θ = π/
3. These x − t plots are colourcoded according to the film thickness with darker colours representing thinner regions. In all casesafter a short transient a nonlinear travelling wave emerges moving from left to right as seen fromthe straight lines in the x − t plane followed by wave troughs and crests. The bottom left panelpresents the evolution of the film thickness measured at a fixed station positioned at the centre ofthe L ∗ = 64 h ∗ -sized periodic domain used for the computations, highlighted by a black dashed linein the figure. The time-periodic signal (after an initial transient) once again verifies the presenceof a nonlinear travelling wave of permanent form. Finally, the bottom right panel presents a com-parison between the saturated interfacial profiles of the resulting nonlinear wave, as obtained fromboth long-wave model and DNS predictions. Notably, from both the transient and the final statecomparisons, there is excellent agreement between the long-wave model and the DNS - this is thecase as long as the assumptions made in the derivation of the model are used for the DNS. This is infact a stringent scenario for all the reduced-order models, as inertial effects and nonlinear featureslead to specific forms of breakdown (which we will soon describe) in this region of the parameterspace. However a weighted-residuals approach is still reliable at this stage. By contrast, a settingweighed down by restrictive assumptions as for the KS equation, leads to quantitatively differentsolutions unless the resulting dynamical behaviour is simple; this is illustrated by the convergenceof the solution to the KS equation to a bimodal travelling wave, instead of the unimodal wave thatboth the weighted residuals and the DNS converge to, see bottom-right panel of figure 2.Given these results, one could question the appropriateness of some of the reduced-order modelsin direct comparisons with DNS and experiments (of course the dynamics supported are rich andthe equations are of fundamental mathematical importance). Notably, both the Benney and theKS equations are valid in very limited parameter regimes making comparisons with experimentsdifficult. Even the weighted-residuals model solutions included in the comparisons with DNS infigure 2 are close to the boundary of their applicability, and numerical solution is already hinderedby stiffness.Simpler mathematical models, however, are key players in mathematical studies and help us topush conceptual boundaries to the point where the developed methodologies can be applied higherup in the model hierarchy. This is the approach taken here and in particular we subsequentlyuse such methodologies in our virtual experiments with the aim of utilizing them in real-worldscenarios, e.g. physical experiments and applications. We should point out that small discrepanciesstill exist between any model and the full DNS. As illustrated, the long-wave model and DNS are inquantitative agreement (even during transient dynamics towards equilibrium coherent structures),hence we are confident that comparisons and hybrid use of the to frameworks across a wide rangeof scenarios are appropriate.More recently, there has been interest in the study of feedback control for the Kuramoto–Sivashinsky equation. Armaou and Christofides [4] explored the control of the zero solution (flatstate) in small domains, while Gomes et al. [17, 18] generalised their results to long domains (wherechaotic behaviour is observed) and to stabilising solutions with a chosen non-uniform interfacialshape. The results in [17, 18] show that any possible solution to the KS equation can be stabilisedusing a finite number of point actuated controls whose strength only depends on the differencebetween the observed and desired interfacial shapes. The number of control actuators dependsonly on the domain length, and the control rule can be computed using a standard pole placementalgorithm [23]. Furthermore, the controls are robust to uncertainty in the problem parameters, aswell as to small changes in the number of controls used. Motivated by the similar linear stability7roperties between the KS equation and the Benney equation (the simplest long-wave model),Thompson et al. [42] studied the control problem for two long-wave models: the Benney equationand the first order weighted residual model, which acted as a test for the robustness of the controlsacross the full hierarchy of models. The authors start by showing that in the unrealistic scenariowhere one can observe the whole interface and actuate everywhere , the simplest proportional controlsof the form f ( x, t ) = − α ( h ( x, t ) − , (12)for some constant α > h ( x, t ) = 1 (or indeed any desired solution H ( x, t ), by replacing 1 by H ( x, t )). The critical valueof α can be computed from linear stability analysis of the Benney or KS equations and it dependsonly on the Reynolds and capillary numbers. It is also shown that the critical α for the Benneyequation is sufficient to obtain linear stability of the weighted residuals model and indeed the fullNavier–Stokes equations, by solving an Orr–Sommerfeld system. Nonlinear stability is confirmedby numerical simulations of the initial value problem. In the more realistic case of a finite numberof observations of the interface and a finite number of control actuators, Thompson et al. [42] useproportional feedback controls of the form f ( x, t ) = − α M (cid:88) j =1 δ ( x − x j )( h ( x j − φ, t ) − , (13)where δ ( · ) is the Dirac delta function, the control actuators are located at the positions x j , j =1 , . . . , M, and observations of the interface are made at x = x j − φ for some displacement φ . Such acontrol protocol was shown to be efficient in stabilising the flat solution h ( x, t ) = 1 for M sufficientlylarge (in practice, M = 5 is usually sufficient for Reynolds and capillary numbers found in relevantflows – see Appendix A). The authors explore other forms of controls, such as the case when thewhole interface is observed (in which case one can use pole placement algorithms, similarly to theones used for the KS equation), or when the number of observations and controls are different, inwhich case these can also be combined using dynamic observers [42, 49].The latter control strategies are the most efficient in stabilising the flat solution for the Benneyequation, but since their design is model-dependent, their applicability across the hierarchy ofmodels is unlikely to be accurate. For this reason, in this paper we chose to use DNS to studythe applicability of full proportional controls (12) and point-actuated proportional controls (13)developed for the long-wave models. We will see that we cannot simply “translate” the controlsdesigned for the long wave models (even the weighted residuals model) into the virtual experimentframework, since there are physical effects that appear at the DNS level which are not mitigated inthe weighted residuals model. However, our study enables us to design a simple adaptation of themodel-based control rules to attain desired control in the DNS, thus acting as a valuable guidingtool within the multi-dimensional parameter space and reducing computational time requirementsby several orders of magnitude. We have constructed a state-of-the-art computational framework in which we can conduct highlyaccurate in silico experiments of a real-world scenario, controlling a falling film down an inclinedplane. This framework does not require any restrictive assumptions and is capable of resolvingall the relevant scales and nonlinearities, thus enabling direct comparisons with real-world physicalexperiments, and indeed theoretical predictions that are devoid of experimental errors and challengesin imaging, measurement and data acquisition. As such, it constitutes a powerful environment to8valuate mathematical model prior to refining the control methodology for specific applications.In addition, we can construct databases containing the entirety of the flow information withoutexperimental restrictions and errors. Importantly, everything above holds in a general setting;however, we selected a classical fluid problem for which a range of well-known hierarchy of reducedmodels, numerical and experimental results are available. This allows us to focus on the mostdelicate aspect of our work, namely the efficiency and accuracy of controls constructed on reducedmodels as they are used in the control of progressively more complex systems, e.g. Navier-StokesDNS. Furthermore, the fluid flow problem is of intrinsic importance and offers a rich landscape ofsolutions and interplay of physical effects and pertains to a wide range of industrial applicationsfrom coating technologies to cooling systems in microchips and multi-physics solutions for heat andmass transfer.Details of the open-source computational platform we used, the
Gerris
Flow Solver [31, 32], areprovided in Appendix C. There we discuss the general direct numerical simulation methodology,details about the discretisation scheme, the volume-of-fluid method used to represent the fluidinterface, as well as technical aspects related to the large scale solution effort, data gathering andpost-processing. A rigorous validation procedure has also been implemented utilizing both convergedand transient model solutions, and also information regarding stability and regime boundaries inthe target parameter space (discussed in detail in the following section).The implementation is based upon monitoring a periodic domain of sufficient length comparedto the film height (typically 64 times the film thickness) while being able to investigate all flowquantities of relevance in an unsimplified setting. This still enables us to import the full controltoolkit presented at the end of subsection 2.2, including either distributed (full-surface) or, morerealistically, point-actuated controls based on discrete interfacial height measurements.Distributed controls are approximated as piece-wise constant strips – numbering between 4 and64 across the length of the periodic domain – attempting to mirror a realistic setup with changeableparts and modular elements. We have noticed that beyond 16 elements the results no longer varywithin the tested conditions, thus amounting to a sufficiently accurate representation of the fullcontinous setup.As in Thompson et al. [42], for the point-actuated scenario we require a special treatment ofthe localised control region. The Dirac δ ( x )-functions introduced in equation (13) are convertedto smooth finite counterparts via s ( x ) = exp[(cos(2 πx ) − /w ], where w denotes the smoothingwindow. As w →
0, we find s ( x ) → δ ( x ), but in practice we choose w to be of O (10 − − − ) topreserve the nature of the intended effect while allowing for an efficient numerical solution of theresulting system given resolution and multi-scale constraints. Here we present our results of controlling a falling liquid film in a DNS framework, as informedby control-theoretical approaches for reduced-order models. Figure 3 encapsulates our multi-stageapproach to the control problem, where our fundamental aim is to be able to accurately and effi-ciently control the dynamics of a complex physical setup. Trial and error in such a delicate andexpensive environment is inefficient, hence modelling and simulation-based approaches are to beused in conjunction to inform a cogent strategy. Full model representations of the target systemstill introduce complexity and efficiency barriers, often offering no advantages apart from shiftingthe difficulty and burden to expensive computational clusters. Our approach is to use a variety ofreduced-order models to construct controls and use the obtained knowledge to pinpoint controllableregimes and guide the heavier simulation machinery to pre-determined regions in parameter space.9igure 3: Multi-layer control methodology outline, from reduced-order models to in silico experi-ments and real-world setup.The computational platform can then be expanded to overcome modelling assumptions and em-ploy the same control strategies in a targeted way, emulating challenging regimes in the real-worldphysical system itself.For our particular model choice of multi-phase system, water and air at room temperature arechosen as reference. Other fluids could be used (see Appendix A for a brief discussion), however wehave opted for a more challenging yet realistic scenario that supports a balance of competing forcesin the flow including inertia that is typically absent in an oil-based flow. The controls are builtbased on linear stability analysis and the numerical solution of a first-order weighted residuals modelwhose suitability has been discussed in subsection 2.2. Similar results can be obtained using differentreduced order models, but with notable restrictions in terms of applicability in the parameter space.Given the surprisingly good agreement between simulations of the reduced-order models andDNS (see Appendix C and Figure 2), we would perhaps expect to be able to use the controlsdeveloped for the model directly in the simulation of the full Navier–Stokes equations. If this wereto work, the methodology would require no complicated data measurements (e.g., interfacial height)which would be needed in real time, and the control could be “pre-coded”, fast, and efficient. Thisis what we call “offline” controls. Alternatively, if unsuccessful, one would need to use feedbackcontrols, which are built using linear stability information from the reduced-order models, but whichrequire readings from the full system. We explore each of these approaches next.
It is tempting, given the potential for complexity reduction, to conduct a full control feasibilitystudy at the level of the reduced-order model and then simply export the resulting controls intothe DNS. Should this work, all values of interest could be stored and relatively easily adapted toflow conditions at the level of the full setup. The initial motivation for this line of investigationfollowed an interesting observation regarding model behaviour; when the control actuation is ofthe form (12), after a very short-lived initial adjustment, the dynamics is steered towards a linearregime and for the rest of the evolution towards a flat state, the interface decays in amplitude andtranslates with a constant speed, keeping the same shape to within a translation and lateral rescaling.Furthermore, the decay in time of its amplitude is exponential since the controlled dynamics arestable and linear. In fact, for the Kuramoto-Sivashinsky equation which is the simplest model for10igure 4: Control results for a film of undisturbed thickness of h ∗ = 125 µ m on an inclined planedescribed by θ = π/
3. (Left) Evolution of the interfacial shape maximum h ∗ − h ∗ for model andtwo DNS control approaches (offline and feedback). (Right) Measured instantaneous velocity basedon the interface maximum location for both model (black line, constant) and DNS (symbols) forthe feedback control case. In both panels the start time of the controls ( t ∗ ≈ .
62 s) is shown via adashed vertical line.falling liquid films presented in section 2.2, it is possible to prove that this decay is exponential andto estimate its decay rate, see [18]. Thus it would seem reasonable that, using normalisation andspatial readjustment (left/right shifts), the precise interfacial profiles would inform a very simplesetup for the proportional control functional f ( x, t ) as defined in eq. (12). This can be summarisedas a separation of variables via f ( x, t ) = g ( x − ct ) · exp( − ωt ), where the shape g would be givenby the nonlinear saturated wave shape moving with a characteristic constant velocity c , while thetime-dependent amplitude (informed by decay rate ω ) could be readily extracted from the model.Making sure of the alignment, i.e. starting with a suitable shift in both space and time, would bethe only sensitive aspect of the procedure, since it does require a reasonable level of synchronisation.A typical result of this procedure is illustrated in red as part of Figure 4 for a relatively mild casegiven by a liquid film of height h ∗ = 125 µ m falling down an plane inclined at an angle of θ = π/ t ≈ .
62 s (this time is indicated bya vertical dashed line), after convergence to a nonlinear saturated shape. The initial evolutionlooks promising, with a rapid decay of the perturbation being observed. This occurs, however,only during a very brief window and is followed by an equally rapid return to the nonlinear waveshape we were aiming to control to begin with (though this is may not necessarily be the caseand a convergence to a different solution is also possible in principle). Albeit a specific example,this behaviour is characteristic of every case study attempted with this strategy. The action ofthe imposed control eventually falls out of sync with the interface due to nonlinearities and thericher physics in the full model. During the respective time, the imposed control continues to decayexponentially and its strength eventually becomes insufficient to affect the flow, which reverts toan uncontrolled configuration. Further tailoring of the imposed control parameters (in terms ofvelocity, amplitude etc.) results in variations around this trend, and without a perfect co-evolution,the precise self-similar interfacial structure is not preserved, resulting in a strategy which is notrobust and ultimately fails.Further details are best explained in contrast to full proportional control, in which we followand use the interface shape at every timestep to inform the construction of f ( x, t ) rather than its11odelled offline counterpart. The results in question are shown in blue in Figure 4. The panel on theleft provides a useful indication that in this case, the decay of the interfacial perturbation is indeedstill exponential in time, however at a slower rate than in the reduced order-model. But arguablyeven more revealing is a study of the horizontal velocity of the interface as a function of time (rightpanel). Even in a successful control scenario, the variation once the control is applied is significantand becomes progressively more sensitive as the perturbation decays in time. There should thus beno surprise that a fixed shape, along with its velocity and a decay rate would be near impossible todetermine reliably in a practical manner, even with perfect information. Any such fixed parameterswould rapidly lead to an out-of-phase behaviour, with the exponentially decaying control ultimatelybecoming null while the system state is still considerably perturbed. Whilst unsuccessful, this studyhas given us valuable insight in terms of the difference in behaviour between the reduced and fullmodels. It has also provided a strong early indication that, while offline controls can only achievea transient positive response, if the capabilities were enhanced towards the integration of dynamicobservations via feedback controls, the desired convergence towards a target state can be recovered.We explore this further in the following subsection. The findings of the previous section indicate that time-dependent flow information needs to beincorporated into the control methodology. Broadly speaking, there are two types of feedbackcontrol we can use: 1. distributed controls, which can act in principle on the entire bottom surfaceof the system and are given by (12), and 2. point-actuated controls given by (13), which are onlydefined at discrete points along the control surface. Both of these can rely on either complete orlimited information on the interfacial shape. There is, however, a clear difference in applicabilitybetween the two. The former is more idealized in view of the need of full information for theinterfacial shape and continuous actuation on the boundary. On the other hand, a practical settingwill only offer the possibility of limited flow information (interfacial shape, velocities, pressure) anda specific range of actuation points are feasible in order to preserve the structural integrity of theboundary. Both scenarios are explored in what follows.For the DNS implementation each of the above methods requires a suitable discretization treat-ment. For distributed controls, we have opted to construct a full-surface actuation method basedon equally-sized segments that would represent actuation strips on the bottom boundary. Forour target parameter regimes, we have noticed results no longer vary beyond 16 such segments,however in most cases 64 have been used in order to ensure smoothness even in more challengingnonlinear scenarios. In the case of point-actuated controls, the theoretical δ − function profiles weretreated following Thompson et al. [42], using smoothed actuator shape functions (with weights w = O (10 − − − )) extending over several computational grid cells in order to suitably incorpo-rate the required variation, whilst maintaining the total flux of fluid in the domain constant via anintegral normalisation.With the above computational framework in place, we need to fix an observation strategy thatis both powerful and practically achievable, with the interfacial height location proving sufficientto satisfy the required criteria. This is sufficient to fully inform the reduced-order models, whileobtaining and processing additional information regarding the flow field on the fly would be tooexpensive (with currently available technology) to allow for a timely control adaptation. Assumingfor simplicity that a full interfacial shape approximation h ( x, t ) is available (this need not be thecase, and limited observers are also possible, see [4, 42]), the control procedure can be summarisedusing the following parameters:1. control strength α , typically a multiplier of the difference in absolute value between the current12nd the target states;2. an upstream or downstream displacement φ , allowing for the possibility of incorporating spa-tial delay information in the actuation procedure;3. for the point-actuated setup, the number of application points M – in this case, we alsoconsider a limited number of observations, equal to the number of actuators.We note that, for the point actuated case, a more complicated (and possibly more efficient)control strategy can be obtained with observations of the full interface. This would involve usingthe full interface information to design an optimal feedback control strategy, which can be achievedusing pole placement or a linear quadratic regulator – this was done in [42], where a feedback controlstrategy was designed for the Benney equation and tested in the weighted residuals model. Further-more, if the number of observations is different than the number of actuators (or if a distributedcontrol were to be designed using a finite number of observations of the interface), similar strategiescan be used to design dynamic observers (see [4, 42]). Both of these strategies are the subjectof current work, as they are harder to translate across the hierarchy of models and are thereforebeyond the scope of this paper. In the distributed control case, for given parameters of the system (such as the Reynolds number Re , capillary number Ca and the inclination angle θ ), we can perform linear stability analysis of thesimplest reduced order model – the Benney equation – and obtain an analytical expression for theminimum value of α needed to obtain stability; it was shown in [42] that these values are sufficientto guarantee linear stability of the full system. In contrast, if the controls are point actuated, theeigenvalues of the linearised system need to be computed numerically. However, the reduced ordermodels still provide valuable information in this case, such as an estimate on the number of controlactuators sufficient to obtain stability. For the Kuramoto–Sivashinsky equation this is shown to belarger than or equal to the number of unstable Fourier modes, see [17, 18], and this was observedto be sufficient for both long wave models in [42].We perform linear analysis for both control cases in Figures 5 and 6. Figure 5 presents theresults for distributed controls; in the left panel we plot the linear stability regions predicted by thefirst order weighted-residuals model for films with undisturbed thickness spanning h ∗ = 125 µ m –200 µ m and inclination angle θ = π/
3. These regions were obtained by computing the eigenvalues ofthe linearised system by including the effects of the spatial delay φ in the analysis in [42, SectionsIII.B,E], and plotting neutral curves determined by the maximum real part of the resulting eigen-values being zero. We note that due to the longwave assumptions, we only considered eigenvaluescorresponding to long-wave modes for the analytical computation. However, since linear stabilitydoes not guarantee that the nonlinear evolution of the controlled system is stabilised, we also per-formed a number of tests whereby we fix the undisturbed thickness h ∗ = 175 µ m and the controlstrength α = 0 .
15, and vary the shifted observer values by − ≤ φ ≤
8, corresponding to the dashedline in the left panel. The evolution of the interfacial maxima as a function of time for this familyof tests is plotted in the right panel, and the plot also includes the subsequent dynamics after thecontrols are switched on at t = 1 . h ∗ = 125 µ, µ m,and for varying α, φ and M the number of actuators (and observers). For this case, we tested thenonlinear evolution for h ∗ = 150 µ m, M = 5, and various values of α , for fixed φ = 4, and of φ for13igure 5: Distributed control study summary for a test case described by films of undisturbedthickness h ∗ = 125 µ m-200 µ m falling down a plane with inclination angle θ = π/ α = 0 .
15 and undisturbedinterface height h ∗ = 175 µ m. The regions above each curve denote areas of expected linear stability.Coloured symbols indicate a successful control scenario, while black circles denote convergence toa non-flat state in the DNS. (Right) Evolution of interfacial maxima as a function of time for thefamily of tests with fixed h ∗ = 175 µm , α = 0 .
15 and different shifted observer values − ≤ φ ≤ α = 0 .
7. The right panel plots the evolution of the interfacial maxima as a function of timefor fixed α and varying φ , where again we include the dynamics after the controls are switched onindicating that the value of φ is crucial in attaining a controlled flat state.As pointed out earlier, prediction of linear stability (or instability) does not guarantee that thetime evolution of the nonlinear system results in convergence to a flat state. Indeed we observethat in the distributed case presented in figure 5, the test described by h ∗ = 175 µ m, α = 0 .
15 and φ = 0 is predicted to be linearly unstable as it sits outside the corresponding neutral curve depictedin orange, i.e., it has at least one eigenvalue with a positive real part. However, the time evolutionshows that the controlled interface evolves towards a flat state. On the other hand, for the samefilm thickness h ∗ = 175 µ m, control strength α = 0 .
15, but when the observation shift is φ = 8,linear stability predicts that controls stabilise the flat solution, however the controlled dynamicsevolves to a non-uniform equilibrium solution instead of the flat target state. In the latter case, wefound (not shown) that the dynamics evolves to a trimodal wave, which corresponds to an unstableshort-wave mode in the reduced-order model that was not considered in the stability plots due toits short-wave nature. This short wave destabilisation is most likely triggered by the fact that φ = 8is a large shift, and hence information propagates from the wrong place. Similarly, in figure 6 weobserve that in the tests for the point actuated case with fixed h ∗ = 150 µ m and M = 5 controls,there are two values for α , for fixed φ = 4, for which linear stability analysis predicts the flat stateto be unstable to perturbations, but for which the nonlinear simulations show that the flat solutionis stabilised. Furthermore, there is one value of φ for fixed α = 0 . h ∗ = 125 µ m-150 µ m falling down a plane with inclination angle θ = π/ M = 3 (dotted), M = 5 (dash-dotted), M = 7 (dashed) and M = 9 (solid), and associated direct numerical simulation results with M = 5controls for fixed h ∗ = 150 µ m and either control strength α = 0 . φ , or fixed shiftedobserver value φ = 4 and varying α . The regions bounded by each curve denote areas of expectedlinear stability. Coloured symbols indicate a successful control scenario, while black circles denoteconvergence to a non-flat state in the DNS. (Right) Evolution of interfacial maxima as a function oftime for this family of tests with undisturbed thickness h ∗ = 150 µm , M = 5 controls, fixed controlstrength α = 0 . − ≤ φ ≤ h ∗ = 150 µ m, and reasonably strong control strength α = 1 . M = 5 application points and a shift of φ = 4 and calculated the discrepancies arising from blow-ing/suction controls as we move up the hierarchy of nonlinearities. Figure 7 presents the horizontalvelocity fields (top panels) and the corresponding vertical velocity fields (bottom panels) for themodel and DNS in the left and middle panels as labelled, along with the difference between the twovelocity fields (right panels) at 0 .
01 s after the control is initialised. The horizontal velocity fieldsare reasonably close to each other, both qualitatively and quantitatively. For the vertical velocityfields, looking beyond the general qualitative and quantitative agreement, the most poignant featureis the stronger vertical banding of the observed structures in the model. The transfer of informationfrom the blowing/suction imposed at the bottom boundary reaches the liquid-gas interfaces morerapidly, while in the case of the DNS a dampening of this effect is observed as we move laterally.These discrepancies originate in the simplifying assumptions that were necessary in the derivationof the reduced-order model but not in the DNS. These include, (a) depth-averaging, (b) neglect ofhigher-order effects leading to what is known as viscous dispersion in falling film flows [33], and (c)inertial effects that are fully present in the DNS, which also make a difference given the consideredmoderate Reynolds number regimes. Interestingly, the most discernible differences occur at thecenter of the fluid region rather than near the interface, and generally take a much less localisedform. Detailed comparisons between model and DNS for particular cases provides insight into theaction of the control mechanism and helps identify predictive bottlenecks, e.g. by overplaying thecontrol effects in the case of thicker films in which DNS indicates that diffusion and stronger inertialfeatures are likely to damp the intended control strength, resulting either in a delayed convergenceto a steady state or insufficient energy input altogether.15igure 7: Comparison between horizontal (top row) and vertical (bottom row) velocity fields forboth reduced-order model (left panels) and direct numerical simulation (central panels) results,with their difference in absolute value shown on the right panels. The snapshots are taken 0 .
01 safter the application of a point actuated control of strength α = 1 . φ = 4 to thedeveloped solution of a liquid film of undisturbed thickness of h ∗ = 150 µ m falling down a planewith inclination angle θ = π/ Up to this point, our results have have focused on the control towards a flat uniform state. However,in many applications such as heat and mass transfer or directed assembly in micro-manufacturing,increasing interfacial area in a precise manner is a desirable feature. We have carried out extensivenumerical computations to model a wide spectrum of scenarios, and focused on the control ofboth shapes and travelling wave velocities of the resulting interfacial profile. The control strategiesimplemented here were found to be successful as long as the target states do not deviate significantlyfrom what should be supported by the underlying physics. For example, surface tension will typicallyoppose interfacial shapes with small wavelengths and large amplitudes, and we found that suchcases result in smoothed out shapes far from the desired target state. It is interesting to notethat short-wave features that overcome such surface tension limitations can be robustly controlledusing externally applied electric fields, for example, with applications in micro- and nano-scale softlithography [34]. However scenarios with primarily unimodal uncontrolled solutions can easily beswitched to higher modes, as described for example by using a target sinusoidal interfacial shapewith wavenumber k = N π/L . A DNS example using distributed controls with strength α = 0 . h ∗ target ( x ∗ ) = sin (cid:16) πx ∗ h ∗ (cid:17) . Provided theactuation process is sufficiently strong, we found that prescribed solutions ranging from interfacialprofiles effectively frozen in space (the shape no longer changes, but the bulk velocities are non-zeroand the flow is moving down the plane subject to gravity), to an accelerated flow are all possibleprovided the energy input is deemed acceptable. In all cases this was achieved with α = O (1) andwould in principle be a viable mechanism to specialise towards more specific applications.The scenario summarised in Figure 8 is highly nonlinear due to the relatively large interfacialperturbation amplitude. Nevertheless, following the initial transient, an exponential decay of the16igure 8: Direct numerical simulation results for feedback control towards a non-uniform sinusoidaltarget solution h ∗ target ( x ∗ ) = sin (cid:16) πx ∗ h ∗ (cid:17) from the developed solution of a liquid film of undisturbedthickness of h ∗ = 225 µ m falling down a plane with inclination angle θ = π/ α = 0 .
5. (Left) L − norm of the error between obtained and target interfacial shape. (Right)Horizontal velocity u ∗ (top) and vertical velocity v ∗ (bottom) at time t ∗ f ≈ α (depending on physical parameter values), local extremain the interfacial profile align with the position of maximum blowing/suction at the centre of theboundary strips. This type of patterning can be achieved at scale and is not restricted to a singletype of surface profile. Microfluidic devices which require different features at different locations (orstages) could in principle take advantage of this flexibility. Finally, we also note that while in theinterface flattening scenarios the control will eventually be reduced close to (but never identically)zero, in non-uniform scenarios it has to be maintained throughout in order to sustain non-trivialshapes.Given the results above, it is natural to question whether we can control the system towardsany designed steady state or saturated solution and what controls would be needed for that. Thediscussion in [42, Section V] shows that we can reach any unstable steady state of our choice, withthe caveat that these are model dependent and consequently the controls obtained cannot be readily“translated” across the hierarchy of models. It is therefore suggested that it is more viable to controltowards non-solutions H ( x ). This can be achieved by adding a forcing term S ( x ) that guaranteesthat H is a (likely unstable) solution of the new forced equation, with the (full) control now being17f the form F = − α ( h ( x, t ) − H ( x )) + S ( x ) . The constant in time forcing S ( x ) can be worked out by solving a linear PDE, in the case of aBenney equation, or a nonlinear PDE for the weighted residuals model – these PDEs are obtainedby replacing H ( x ) in the respective model and obtaining a differential equation for S . Once S isfound, it is possible to show that the steady state H is reached for large α – which in itself can alreadybe an impediment to the controllability of the system. We note that the forcing S ( x ) needs to be constantly applied for the controls to work, even after the steady state H is reached (otherwise H isno longer a solution of the equation and the interface becomes unstable again). Therefore knowing S can inform us of the feasibility of controls for a chosen H : if S ( x ) is prohibitively large whencompared to, e.g., the characteristic speed of the system, the controls are no longer “imposed” butare rather “enforced” and one could argue that they are not a realistic strategy. This study highlights a route towards the application of robust control strategies derived rigorouslyat the level of reduced-order mathematical models towards realistic “in silico” experiments (pre-sented in the form of direct numerical simulations). This allows, for the first time, the study offeasibility of both offline and feedback approaches, alongside a systematic scrutiny of what systemfeatures are retained as the hierarchy in nonlinearity is traversed.Using the classical setup of a falling liquid film with imposed blowing/suction as active controlmechanism, we find that offline (or statically tailored) strategies provide a poor candidate for prac-tical purposes due to the difficulty in synchronising exponentially decaying solutions with highlynonlinear systems. By contrast, feedback strategies (both distributed and point-actuated) provide aviable setup for convergence to desired states which can be simple/trivial (flat interfaces desirable incoating technologies) or spatially non-uniform (with increased interfacial area, of use in heat/masstransfer applications, for example). The analytically-informed high performance computing frame-work has been built on the basis of interrogating both stability and control theoretical results atthe reduced-order model level. These included linear stability predictions on the minimum controlstrength α needed to stabilise a flat system, and on the best value for the spatial delay φ in thedistributed case, complemented with the number of control actuators and their location for pointactuated controls. The theoretical results provided a priori quantification on the behaviour of thecontrols which proved to be substantiated at the full model level. After confirming the validity ofthe model predictions, we could build on this knowledge and propagate similar strategies towardsregimes outside the reach of the analytical (and numerical) methods.In summary, the proposed methodology (1) is unique in its multi-faceted approach towardsefficient modelling, with both analytical and computational components acting in tandem, (2) safe-guards against the overly simplistic use of offline control methods, (3) offers a paradigm shift interms of control information transition between modelling approaches, (4) informs viable validatedcontrol strategies within the range of applicability of analytical methods, and (5) provides access towide regions in parameter space, and is a natural departure point for feedback control applicationseven in cases beyond the validity of the asymptotic models, with capabilities of informing a power-ful integrated computational platform and ensuring a step change in transitioning towards real-lifesystems.The particular interfacial flow setup used to demonstrate this interplay between techniques is awell-known and canonical one, however, more importantly, the approach has tremendous potentialfor generalisation towards complex systems in fluid mechanics and other fields. The toolkit has18lear extensions in terms of classical control theoretical aspects, for example full feedback control,combining different observations to obtain a more efficient control rather than associating one ob-servation to its own actuator, and using optimal control techniques to design controls which arerobust to the hierarchy of models. Extensions in the application area include multi-physics scenar-ios such as electrohydrodynamics and soft lithography, heat and mass transfer, photo-excitation,and acoustics, to mention a few. In such applications the control design procedure is most oftenperformed either very early (and thus with inaccurate propagation of information and insight) orvery late (after expensive trial and error in early stages) in the process, with the feedback controltechnology becoming of limited value, especially in view of further generalisation and applicability.The proposed hierarchical control design approach is anticipated to considerably reduce expensiveexperimental and manufacturing pipelines and is well suited towards specialisation within a widevariety of physical, engineering and cross-disciplinary contexts. R.C. gratefully acknowledges the resources and continued support of the Imperial College Lon-don Research Computing Service. He also notes the important role of the University of OxfordMathematical Institute funding in the form of the Hooke Research Fellowship, as well as resourcesfrom the UK Fluids Network (EPSRC Grant EP/N032861/1), which have enabled significant cross-institutional interactions. S.G. is supported by the Leverhulme Trust through the Early CareerFellowship ECF-2018-536. D.T.P. was partly supported by EPSRC Grant EP/L020564/1.
A Parameter values
In all our numerical calculations we considered parameters pertaining to a water-air system whichis easily realisable experimentally (see e.g. [20, 14]). We stress, however, that the approachesdiscussed here are applicable to liquid-gas systems in general, as long as the assumptions behindthe modelling framework hold. The inclination angle is fixed at θ = π and the fluid thicknessis varied from h ∗ = 125 µ m to h ∗ = 225 µ m, recalling that the ( · ) ∗ notation is used to denotedimensional quantities. The periodic domain is constructed with an aspect ratio of 64 : 1, suchthat L ∗ = 64 h ∗ , which is in accordance with the long-wave assumption. These settings translate to moderate Reynolds number flows where the uncontrolled dynamics exhibits a range of behaviours,including nonlinear travelling waves, trains of solitary wave pulses, three-dimensional dynamicsand spatiotemporal chaos. The thinnest film h ∗ = 125 µ m is just beyond the threshold for linearstability and hence resulting in the formation of mildly nonlinear travelling waves, while the upperlimit h ∗ = 225 µ m pertains to sufficiently thick films in more inertially-dominated flow regimes sothat the reduced-order models start losing their validity and computations become restrictive dueto numerical stiffness.Relevant parameters for the liquid are its density ρ ∗ l = 998 kg/m and dynamic viscosity µ ∗ l =8 . × − kg/ms , while the equivalent gas parameters are ρ ∗ g = 1 . kg/m and µ ∗ g = 1 . × − kg/ms . The water-air surface tension coefficient is taken to be γ ∗ = 0 . kg/s , while theacceleration of gravity is given by g ∗ = 9 . m/s . For the ranges of interface heights we consider inour numerical experiments, we obtain the reference velocities and dimensionless groupings presentedin Table 1 below, with the liquid values considered as reference (see details in subsection 2.2).19 ∗ ( µm ) U s ( m/s ) Re Ca
125 0 . . . . . . . . . . . . . . . Re and Capillary number Ca . B Numerical solution of the reduced order models
The solutions to the reduced order models presented in Figure 2 were obtained as follows. TheKuramoto-Sivashinsky equation was solved in an appropriately rescaled periodic domain x ∈ [0 , π ]using a pseudo-spectral method for the spatial discretisation based on a Fourier series decompositionof the solution. To time-step, we used second order backward differentiation formulae (see [1])with a sufficiently small time step. The weighted residuals model was solved using a pseudo-spectral spatial discretisation and a variable-step, variable-order time-stepping scheme based on thenumerical differentiation formulas (NDFs) of orders 1 to 5, provided by Matlab ’s inbuilt function ode15s . C Direct numerical simulations and validation
We have developed the code for the direct numerical simulations conducted in this investigationas part of the open-source solver
Gerris [31, 32]. The main implementation, which can be foundat http://gfs.sourceforge.net , is ideally suited for interfacial flows in both fundamental andapplied contexts and has seen adoption and growth over the past two decades.Here we briefly sketch the approach for our particular problem and refer the reader to [32] forfurther details. The equations of motion are ρ ( ∂ t u f + u f · ∇ u f ) = −∇ p + ∇ · (2 µ D ) + σκδ s n + F e ,∂ t ρ + ∇ · ( ρ u f ) = 0 , (14) ∇ · u f = 0 , where D is the rate of strain tensor D ij = (1 / ∂u i /∂x j + ∂u j /∂x i ). We stress that the fullNavier-Stokes equations are solved in both the liquid and the gas, with the decoration ( · ) f being aplaced holder for both gas-related - ( · ) g - or liquid-related - ( · ) l - quantities. In the interfacial regionbetween the two, discrete counterparts of the normal and tangential stress balances, the kinematiccondition, as well as continuity of velocities are imposed. All interfacial forces are transferred tothe momentum equations in what is commonly called the “one-fluid” formulation. The physicalproperties describing each fluid (density, viscosity, permittivity etc.) are included by singular distri-butions and the same set of equations (14) accounts for the entire domain. The Dirac distribution δ s isolates the surface tension effects to the interface alone, and any external forces such as gravityare included via the F e term.The relevant physical properties of the fluids (such as density and viscosity) are represented interms of a volume fraction c ( x , t ), where c is a generic colour function which takes the value 0 inone fluid (e.g. the gas) and 1 in the other (e.g. the liquid, our reference fluid). More specifically we20rite ρ ( c ) ≡ cρ l + (1 − c ) ρ g , µ ( c ) ≡ cµ l + (1 − c ) µ g , (15)Under this treatment, a density equation (the other properties are treated in a similar manner) ofthe general form ρ t + ∇ · ( ρ u f ) = 0 (16)becomes c t + ∇ · ( c u f ) = 0 , (17)which is solved for c and the results substituted into (15). The value of c is interpolated across theinterface by introducing a small transition layer in its vicinity to smooth the variation of quantitiesfrom one region to the other, thus relaxing an otherwise singular transition.The discretisation schemes are 2 nd -order adaptive in both time and space, with a series of criteriaranging from refinement settings to capillarity- and inertia-based constraints dictating the requiredtimestep to advance the solution within a tolerated error bound. The critical numerical step inthe process is a multi-level Poisson solver for pressure (following projection steps). Efficient solu-tions thereof are facilitated by excellent spatial adaptivity which can span 3 orders of magnitude incell sizes, thus reducing the number of degrees of freedom considerably and allocating them whereneeded. In our particular case, we employed changes in the velocity fields and vorticity, as wellas interfacial location as refinement criteria, leading to systems with O (10 ) cells. Furthermore,the structured mesh based on an underlying quadtree (octree in 3D) setup results in strong par-allelisation properties. Most numerical experiments have been run on 4 −
16 CPUs on local highperformance computing facilities at Imperial College London and required roughly 100 − O (10 − ) in the constructed norms pertaining to O (1) changes inthe flow properties), which allowed initial resolution and refinement criteria tailoring. Secondly,the reason for choosing this particular application lies in the fact that both model solutions (re-visit subsection 2.2) and previous numerical and experimental solutions, e.g. [14], are available forcomparison. The former were employed as stringent target cases for relatively small thicknesseswere the models are known to be valid, while the latter resources allowed additional testing in morechallenging and often inertially-dominated regimes outside the reach of the models.21he summary in Figure 2 pertains to a scenario which constitutes the most difficult test casefor which modelling efforts still converge to a solution. Similar studies have been performed forsmaller film thicknesses. We noticed that the saturated nonlinear wave shapes compared very wellin all such cases and that the evolution itself (based on tracking the interfacial height at specificmeasurement stations) matches surprisingly well even at moderate Reynolds numbers (of up to Re ≈ . References [1] Georgios Akrivis, Demetrios T Papageorgiou, and Yiorgos-Sokratis Smyrlis. Linearly implicitmethods for a semilinear parabolic system arising in two-phase flows.
IMA journal of numericalanalysis , 31(1):299–321, 2011.[2] Thomas G Anderson, Radu Cimpeanu, Demetrios T Papageorgiou, and Peter G Petropoulos.Electric field stabilization of viscous liquid layers coating the underside of a surface.
Phys. Rev.Fluids , 2:054001, 2017.[3] Denise L Angwin and Howard Kaufman. Image restoration using reduced order models.
SignalProcessing , 16(1):21–28, 1989.[4] Antonios Armaou and Panagiotis D Christofides. Feedback control of the kuramoto–sivashinskyequation.
Physica D: Nonlinear Phenomena , 137(1-2):49–61, 2000.[5] Peter Benner, Serkan Gugercin, and Karen Willcox. A survey of projection-based model re-duction methods for parametric dynamical systems.
SIAM review , 57(4):483–531, 2015.[6] DJ Benney. Long waves on liquid films.
Journal of mathematics and physics , 45(1-4):150–155,1966.[7] MG Blyth and Andrew P Bassom. Flow of a liquid layer over heated topography.
Proceedings ofthe Royal Society A: Mathematical, Physical and Engineering Sciences , 468(2148):4067–4087,2012.[8] MG Blyth and C Pozrikidis. Effect of surfactant on the stability of film flow down an inclinedplane.
Journal of Fluid Mechanics , 521:241–250, 2004.[9] Nikolai WF Bode, Armel U Kemloh Wagoum, and Edward A Codling. Human responses tomultiple sources of directional information in virtual crowd evacuations.
Journal of The RoyalSociety Interface , 11(91):20130904, 2014.[10] V Bontozoglou. The effect of adsorption modeling on the stability of surfactant-laden liquidfilm flow.
Acta Mechanica , 229:535–547, 2018.[11] Steven L Brunton, Bernd R Noack, and Petros Koumoutsakos. Machine learning for fluidmechanics.
Annual Review of Fluid Mechanics , 52:477–508, 2020.2212] Young-Pil Choi, Dante Kalise, Jan Peszek, and Andr´es A Peters. A collisionless singularcucker-smale model with decentralized formation control. arXiv preprint arXiv:1807.05177 ,2018.[13] Pierre Collet, Jean-Pierre Eckmann, Henri Epstein, and Joachim Stubbe. A global attract-ing set for the kuramoto-sivashinsky equation.
Communications in mathematical physics ,152(1):203–214, 1993.[14] Fabian Denner, Alexandros Charogiannis, Marc Pradas, Christos N Markides, Berend GM vanWachem, and Serafim Kalliadasis. Solitary waves on falling liquid films in the inertia-dominatedregime.
Journal of Fluid Mechanics , 837:491–519, 2018.[15] PH Gaskell, PK Jimack, M Sellier, HM Thompson, and MCT Wilson. Gravity-driven flowof continuous thin liquid films on non-porous substrates with topography.
Journal of FluidMechanics , 509:253–280, 2004.[16] A Georgantaki, M Vlachogiannis, and V Bontozoglou. Measurements of the stabilisation ofliquid film flow by the soluble surfactant sodium dodecyl sulfate (sds).
International Journalof Multiphase Flow , 86:28–34, 2016.[17] SN Gomes, Marc Pradas, Serafim Kalliadasis, DT Papageorgiou, and GA Pavliotis. Controllingspatiotemporal chaos in active dissipative-dispersive nonlinear systems.
Physical Review E ,92(2):022912, 2015.[18] Susana N Gomes, Demetrios T Papageorgiou, and Grigorios A Pavliotis. Stabilizing non-trivial solutions of the generalized kuramoto–sivashinsky equation using feedback and optimalcontrollighthill–thwaites prize.
IMA Journal of Applied Mathematics , 82(1):158–194, 2017.[19] Anders Johansson, Michael Batty, Konrad Hayashi, Osama Al Bar, David Marcozzi, and Ziad AMemish. Crowd and environmental management during mass gatherings.
The Lancet infectiousdiseases , 12(2):150–156, 2012.[20] Serafim Kalliadasis, Christian Ruyer-Quil, Benoit Scheid, and Manuel Garc´ıa Velarde.
Fallingliquid films , volume 176. Springer Science & Business Media, 2011.[21] G Karapetsas and V Bontozoglou. The primary instability of falling films in the presence ofsoluble surfactants.
Journal of Fluid Mechanics , 729:123–150, 2013.[22] G Karapetsas and V Bontozoglou. The role of surfactants on the mechanism of the long-waveinstability in liquid film flows.
Journal of Fluid Mechanics , 741:139–151, 2014.[23] Jaroslav Kautsky, Nancy K Nichols, and Paul Van Dooren. Robust pole assignment in linearstate feedback.
International Journal of control , 41(5):1129–1155, 1985.[24] Petar Kokotovic, Hassan K Khali, and John O’reilly.
Singular perturbation methods in control:analysis and design , volume 25. Siam, 1999.[25] Boris Kramer, Benjamin Peherstorfer, and Karen Willcox. Feedback control for systems withuncertain parameters using online-adaptive reduced models.
SIAM Journal on Applied Dy-namical Systems , 16(3):1563–1586, 2017.[26] Yoshiki Kuramoto and Toshio Tsuzuki. On the formation of dissipative structures inreaction-diffusion systems: Reductive perturbation approach.
Progress of Theoretical Physics ,54(3):687–699, 1975. 2327] Toni Lassila, Andrea Manzoni, Alfio Quarteroni, and Gianluigi Rozza. Model order reductionin fluid dynamics: challenges and perspectives. In
Reduced Order Methods for modeling andcomputational reduction , pages 235–273. Springer, 2014.[28] Seungjoon Lee, Mahdi Kooshkbaghi, Konstantinos Spiliotis, Constantinos I Siettos, and Ioan-nis G Kevrekidis. Coarse-scale pdes from fine-scale observations via machine learning.
Chaos:An Interdisciplinary Journal of Nonlinear Science , 30(1):013141, 2020.[29] Kwang-Kyo Oh, Myoung-Chul Park, and Hyo-Sung Ahn. A survey of multi-agent formationcontrol.
Automatica , 53:424–440, 2015.[30] Demetrios T Papageorgiou and Yiorgos S Smyrlis. The route to chaos for the kuramoto-sivashinsky equation.
Theoretical and Computational Fluid Dynamics , 3(1):15–42, 1991.[31] S. Popinet. Gerris: A tree-based adaptive solver for the incompressible Euler equations incomplex geometries.
Journal of Computational Physics , 190:572–600, 2003.[32] S. Popinet. An accurate adaptive solver for surface-tension-driven interfacial flows.
Journal ofComputational Physics , 228:5838–5866, 2009.[33] Christian Ruyer-Quil and Paul Manneville. Improved modeling of flows down inclined planes.
The European Physical Journal B-Condensed Matter and Complex Systems , 15(2):357–369,2000.[34] Eric Sch¨affer, Thomas Thurn-Albrecht, Thomas P Russell, and Ullrich Steiner. Electrohydro-dynamic instabilities in polymer films.
EPL (Europhysics Letters) , 53(4):518, 2001.[35] M Sellier. Inverse problems in free surface flows: a review.
Acta Mechanica , 227:913–935, 2016.[36] Alex Shum, Kirsten Morris, and Amir Khajepour. Direction-dependent optimal path planningfor autonomous vehicles.
Robotics and Autonomous Systems , 70:202–214, 2015.[37] GI Sivashinsky. Nonlinear analysis of hydrodynamic instability in laminar flames—i. derivationof basic equations.
Acta astronautica , 4:1177–1206, 1977.[38] Gregory I Sivashinsky. Instabilities, pattern formation, and turbulence in flames.
AnnualReview of Fluid Mechanics , 15(1):179–199, 1983.[39] Gregory I Sivashinsky and DM Michelson. On irregular wavy flow of a liquid film down avertical plane.
Progress of theoretical physics , 63:2112–2114, 1980.[40] Eitan Tadmor. The well-posedness of the kuramoto–sivashinsky equation.
SIAM Journal onMathematical analysis , 17(4):884–893, 1986.[41] Alice B Thompson, Susana N Gomes, Fabian Denner, Michael C Dallaston, and Serafim Kalli-adasis. Robust low-dimensional modelling of falling liquid films subject to variable wall heating.
Journal of Fluid Mechanics , 877:844–881, 2019.[42] Alice B Thompson, Susana N Gomes, Grigorios A Pavliotis, and Demetrios T Papageorgiou.Stabilising falling liquid film flows using feedback control.
Physics of Fluids , 28(1):012107,2016.[43] Alice B Thompson, Dmitri Tseluiko, and Demetrios T Papageorgiou. Falling liquid films withblowing and suction.
Journal of Fluid Mechanics , 787:292–330, 2016.2444] Gaurav Tiwari, Ruchi Tiwari, Birendra Sriwastawa, L Bhati, S Pandey, P Pandey, andSaurabh K Bannerjee. Drug delivery systems: An updated review.
International journalof pharmaceutical investigation , 2(1):2, 2012.[45] Ruben J Tomlin, Radu Cimpeanu, and Demetrios T Papageorgiou. Instability and drippingof electrified liquid films flowing down inverted substrates.
Physical Review Fluids , 5:013703,2020.[46] Ruben J Tomlin and Susana N Gomes. Point-actuated feedback control of multidimensionalinterfaces. arXiv preprint arXiv:1901.09223 , 2019.[47] Dmitri Tseluiko, Mark G Blyth, and Demetrios T Papageorgiou. Stability of film flow over in-clined topography based on a long-wave nonlinear model.
Journal of Fluid Mechanics , 729:638–671, 2013.[48] Dmitri Tseluiko and Demetrios T Papageorgiou. Wave evolution on electrified falling films.
Journal of Fluid Mechanics , 556:361–386, 2006.[49] J Zabczyk.
Mathematical Control Theory: An Introduction Birkh¨auser Boston . 1992.[50] Pan Zhao, Jiajia Chen, Yan Song, Xiang Tao, Tiejuan Xu, and Tao Mei. Design of a controlsystem for an autonomous vehicle based on adaptive-pid.