Pancake making and surface coating: optimal control of a gravity-driven liquid film

TL;DR

This study explores optimal surface kinematics to achieve uniform liquid film spreading during solidification, using advanced control methods to improve film uniformity in processes like pancake making and surface coating.

Contribution

It introduces an optimal control framework with adjoint equations for better film uniformity, surpassing previous harmonic parameterization approaches.

Findings

Optimal control significantly improves film uniformity.

Adjoint method reduces the number of evaluations needed.

Harmonic parameterization yields moderate improvements.

Abstract

This paper investigates the flow of a solidifying liquid film on a solid surface subject to a complex kinematics, a process relevant to pancake making and surface coating. The flow is modeled using the lubrication approximation, with a gravity force whose magnitude and direction depend on the time-dependent orientation of the surface. Solidification is modeled with a temperature-dependent viscosity. Because the flow eventually ceases as the liquid film becomes very viscous, the key question this study aims to address is: what is the optimal surface kinematics for spreading the liquid layer uniformly? Two methods are proposed to tackle this problem. In the first one, the surface kinematics is assumed a priori to be harmonic and parameterized. The optimal parameters are inferred using the Monte-Carlo method. This ``brute-force'' approach leads to a moderate improvement of the film uniformity compared to the reference case when no motion is imposed to the surface. The second method is formulated as an optimal control problem, constrained by the governing partial differential equation, and solved with an adjoint equation. Key benefits of this method are that no assumption is made on the form of the control, and that significant improvement in thickness uniformity are achieved with a comparatively smaller number of evaluations of the objective function.

Tables

Table 1. Table 1: Parameter table.

Parameter	$\rho$	$k$	${\mu }_0$	$\alpha$	$R$
	(kg m-3)	(W m-1 K-1)	(Pa s)	(oC-1)	(m)
Value	1000	1	$1.668$	$2.56\times {10}^{-2}$	$1.5\times {10}^{-1}$
Parameter	$V_D$	$h_c$	$T_{\infty }$	$T_s$	$h_{\text{opt}}=V_D/\left(\pi R^2\right)$
	(m3)	(W m-2 K-1)	(oC)	(oC)	(m)
Value	$1.414\times {10}^{-4}$	variable	200	20	$2\times {10}^{-3}$

Table 2. Table 2: Final uniformity and smoothness. Improvements in uniformity and smoothness follow the same trend.

	Final uniformity $𝒰\left(t_f\right)$	Final smoothness $𝒮\left(t_f\right)$
	${\iint }_D{\left(h(𝐱,t_f)-h_{\text{opt}}\right)}^2\text{d}𝐱$	${\iint }_D{\left|\nabla h(𝐱,t_f)\right|}^2\text{d}𝐱$
	(m4)	(m2)
Uncontrolled	$6.8\times {10}^{-8}$	$3.6\times {10}^{-4}$
Best harmonic (MC)	$5.5\times {10}^{-8}$	$1.4\times {10}^{-4}$
2nd best harmonic (MC)	$4.3\times {10}^{-8}$	$1.6\times {10}^{-4}$
Adjoint $\gamma ={10}^{-6}$	$2.4\times {10}^{-8}$	$0.31\times {10}^{-4}$
Adjoint $\gamma ={10}^{-8}$	$1.2\times {10}^{-8}$	$0.20\times {10}^{-4}$

Equations

$$\frac{\partial }{\partial z}\left(\mu \frac{\partial u}{\partial z}\right)$$

$$\frac{\partial }{\partial z}\left(\mu \frac{\partial v}{\partial z}\right)$$

$$0$$

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}$$

$$\frac{d^{2}T}{d{z}^2}$$

$${T(z=0)=T_s,k\frac{\partial T}{\partial z}|}_{z=h}=h_c\left(T_{\infty }-T(h)\right),$$

$$\mu (T)={\mu }_0{e}^{-\alpha T}.$$

$$T(x,y,z,t)=A^{\prime }\frac{z}{h}+T_s,\text{where }{A}^{\prime }=\frac{h_c\left(T_{\infty }-T_s\right)}{\frac{k}{h}+h_c}.$$

$$\mathbf{v}(x,y,z,t)$$

$$\mathbf{K}(x,y,t)$$

$$\frac{\partial h}{\partial t}+\nabla \cdot \mathbf{q}=0,$$

$$\mathbf{q}(x,y,t)=-\frac{\rho g{e}^{\alpha T_s}}{{\mu }_0}\frac{h^3}{{\left(\alpha A^{\prime }\right)}^3}\left(2e^{\alpha A^{\prime }}-{\left(\alpha A^{\prime }\right)}^2-2\alpha A^{\prime }-2\right)\mathbf{K}.$$

$$e_a\frac{{\partial }^{2}h}{\partial t^2}+d_a\frac{\partial h}{\partial t}+\nabla \cdot \left(-c\nabla h-\alpha h+\gamma \right)+\beta \cdot \nabla h+ah=f,$$

$$\mathcal{U}(t)={\iint }_D{\left(h(x,y,t)-h_{\text{opt}}\right)}^2\text{d}x\text{d}y.$$

$$\theta (t)=A_1\sin \left(\frac{2\pi }{T_1}t\right),\beta (t)=A_2\sin \left(\frac{2\pi }{T_2}t\right),$$

$$\theta (t)=A_1\sin \left(\frac{2\pi }{T_1}t\right),\beta (t)=A_2\cos \left(\frac{2\pi }{T_2}t\right),$$

$${\mathcal{J}}_t(h,\mathbf{c})$$

$$=\mathcal{U}(t_f)+\frac{\gamma }{t_f}\mathcal{C},$$

$$\underset{\mathbf{c}(t)}{\min }{\mathcal{J}}_t(h,\mathbf{c})\text{mbox}subjectto(\text{ref}eq:governing),$$

$${\mathcal{L}}_t(h,\mathbf{c},h^{†})=$$

$$=$$

$$-{\int }_I{\iint }_D{h}^{†}\left(\frac{\partial h}{\partial t}+\nabla \cdot \mathbf{q}\right)\mathrm{d}\mathbf{x}\mathrm{d}t,$$

$$\frac{\partial f(s)}{\partial s}\delta s=\underset{ϵ\to 0}{\lim }\frac{f(s+ϵ\delta s)-f(s)}{ϵ}\forall \delta s.$$

$$\frac{\mathrm{d}{\mathcal{L}}_t}{\mathrm{d}\mathbf{c}}=\frac{\partial {\mathcal{L}}_t}{\partial \mathbf{c}}+\frac{\partial {\mathcal{L}}_t}{\partial h}\frac{\mathrm{d}h}{\mathrm{d}\mathbf{c}}+\frac{\partial {\mathcal{L}}_t}{\partial h^{†}}\frac{\mathrm{d}{h}^{†}}{\mathrm{d}\mathbf{c}},$$

$$\frac{\mathrm{d}{\mathcal{L}}_t}{\mathrm{d}\mathbf{c}}=\frac{\partial {\mathcal{L}}_t}{\partial \mathbf{c}}$$

$$\frac{\mathrm{d}{\mathcal{J}}_t}{\mathrm{d}\mathbf{c}}=\frac{\partial {\mathcal{L}}_t}{\partial \mathbf{c}}.$$

$$\frac{\mathrm{d}{\mathcal{J}}_t}{\mathrm{d}\mathbf{c}}$$

$$\displaystyle\mathbf{q}\simeq-\dfrac{e^{\alpha T_{s}}\rho g}{3\mu_{0}}h^{3}\left({\cos\theta}{\cos\beta}\bm{\nabla}h-{\left(\begin{array}[]{c}{\sin\theta}\\ {\sin\beta}\end{array}\right)}\right).$$

$$\displaystyle\dfrac{\partial\left(\bm{\nabla}\bm{\cdot}\mathbf{q}\right)}{\partial\theta}=W\bm{\nabla}\bm{\cdot}\left[h^{3}\left({\sin\theta}{\cos\beta}\bm{\nabla}h+\left(\begin{array}[]{c}{\cos\theta}\\ 0\end{array}\right)\right)\right],$$

$$\displaystyle\dfrac{\partial\left(\bm{\nabla}\bm{\cdot}\mathbf{q}\right)}{\partial\beta}=W\bm{\nabla}\bm{\cdot}\left[h^{3}\left({\cos\theta}{\sin\beta}\bm{\nabla}h+\left(\begin{array}[]{c}0\\ {\cos\beta}\end{array}\right)\right)\right].$$

$$\displaystyle\dfrac{\partial{h^{\dagger}}}{\partial t}+W\bm{\nabla}\bm{\cdot}\left(h^{3}{\cos\theta}{\cos\beta}\bm{\nabla}{h^{\dagger}}\right)-3Wh^{2}\left({\cos\theta}{\cos\beta}\bm{\nabla}h-{\left(\begin{array}[]{c}{\sin\theta}\\ {\sin\beta}\end{array}\right)}\right)\bm{\cdot}\bm{\nabla}{h^{\dagger}}=0,$$

$$h^{†}(\mathbf{x},t_f)=2(h(\mathbf{x},t_f)-h_{\text{opt}}),$$

$$\nabla h^{†}\cdot {\mathbf{n}}_D=0\text{mbox}on\partial D.$$

$${\mathcal{J}}_r(h,\mathbf{c})$$

$$={\int }_I\mathcal{U}(t)\mathrm{d}t+\gamma \mathcal{C}.$$

$$\underset{\mathbf{c}(t)}{\min }{\mathcal{J}}_r(h,\mathbf{c})\text{mbox}subjectto(\text{ref}eq:governing),$$

Full text

I Introduction

One of the motivations for the work presented here is the process of crêpe making. Crêpe making involves pouring a fixed amount of batter on a hot pan, letting or forcing the batter to spread on the hot pan to obtain optimal coverage, and letting the batter cook. For crêpes, optimal coverage means uniformly thin, hole-free, and perfectly circular. Achieving this goal can however be quite challenging since as the batter spreads, it cooks at the same time and if the pan is left horizontal, the batter tends to solidify before reaching uniformly the rim of the pan. There are two main strategies to circumvent this issue. The first strategy involves using a blade to force spread a uniform layer of batter on the pan in a process reminiscent of blade coating. The other strategy involves tilting the pan in a swirling motion and forcing the batter to spread preferentially in the downslope direction of the span. As soon as this tilting is initiated, the axial-symmetry of the problem is broken and one cannot help but wonder what the optimal swirling motion is to achieve an optimally thin and round crêpe. This process has other important implications in chocolate manufacturing technologies, the coating of surfaces, or the production of thin elastic shells Lee16 .

The key physical phenomena underpinning crêpe making involves the interaction of the liquid layer with the substrate kinematics and the solidification of the liquid layer. These are relevant in other related contexts of science and engineering. An archetypical illustration of the interaction between substrate kinematics and a liquid layer is spin coating for which a liquid is first deposited on a substrate which is then spun at high angular velocity to produce a thin liquid film which is subsequently cured, see for example Emslie58 ; Meyerhofer78 ; Sahu10 . In that context, the fluid is driven by the centrifugal forces and the substrate kinematics is a simple rotation about the vertical axis, possibly with variable angular velocity. In contrast, in the context of crêpe making, the batter is driven by gravity and the pan kinematics is much more complicated since it involves transient rotation around multiple axes. To the best of the authors’ knowledge, modeling such flows has not been attempted before. Liquid layer solidification is relevant in a number of processes including coating Weinstein04 , ice accretion on surfaces Myers04 ; Moore17 , paint drying Schwartz10 ; Gaskell06 . The present work builds mostly on the literature related to geophysical flow solidification such as lava flows for which the viscosity is highly temperature dependent Bercovici94 ; Sansom04 and cooling is associated with eventual flow arrest.

A key question the current contribution is aiming to address is what is the optimal pan kinematics to achieve a desired film profile? Such problems fall in the realm of optimal control in fluid mechanics for which substantial literature already exist, see Gunzburger02 , for example. Applications of optimal control in thin layer flows are reviewed in Sellier16 . Relevant contributions include that of Grigoriev where the author finds optimal heating strategies to actively suppress evaporative instabilities in thin liquid films Grigoriev02 , that of Sellier and Panda where the optimal substrate shape to control the free surface of a liquid film is inferred Sellier10 , or that of Papageorgiou et al. where the authors find the optimal source/sink distribution to control and stabilize falling liquid films Thompson16 ; Tomlin18 . Because thin liquid films are commonly described by the long-wave approximation which leads to second-order or fourth-order parabolic partial differential equations (PDEs), depending on whether the effect of surface tension are prevalent or not, it is natural to apply the standard framework of optimal control of PDEs as described in Troltzsch10 ; Neitzel09 . This is what the present contribution proposes to do.

The paper is organized as follows. Section II describes the mathematical model for the spreading of a non-isothermal gravity current and the numerical solution strategy. Results illustrating the gravity current dynamics (solution of the “direct” problem) are provided in Section III. The formulation of the optimal control problem (“adjoint” problem and optimization) is given in Section IV with the corresponding results in Section V. The paper concludes with discussion and conclusions in Section VI.

II Description of the mathematical model

II.1 Mathematical model

The liquid layer is treated as a non-isothermal gravity current spreading over a surface whose inclination varies over time, see figure 1. In order to model a solidification process, the effect of the phase transition on the dynamics is replaced by a temperature-dependent viscosity, and the layer is assumed to remain fluid at all times.

Gravity currents have received considerable attention in the past and an excellent account of the current knowledge can be found in Leal10 . For an incompressible liquid film whose characteristic thickness is much smaller than its characteristic lateral dimensions $H/L=\epsilon\ll 1$ (lubrication theory) and whose motion at characteristic velocity $U$ is of small Reynolds number $\mbox{{Re}}=\rho UL/\mu\ll 1$ (inertial effects much smaller than viscous effects), the conservation equations for momentum, mass and energy reduce to:

[TABLE]

where $(u,v,w)$ are the components of the velocity vector, $p$ the pressure, $T$ the temperature, $\mu(T)$ the temperature-dependent dynamic viscosity, $\rho$ the fluid density, and $(g_{x},g_{y},g_{z})$ the components of the gravity acceleration vector. The coordinate system $(x,y,z)$ is attached to the substrate. The liquid layer is bounded by the $xy$ plane from below and the free surface $h(x,y,t)$ from above.

The substrate rotation is assumed to be sufficiently slow that the centrifugal force and Coriolis force may be neglected (this requires the characteristic rotation rate to be small compared to $\sqrt{\mu U/(\rho LH^{2})}$ and to $\mu/(\rho H^{2})$, respectively; we have verified that this is satisfied at almost all times for the substrate motions considered in this study; see Sec. III).

Eq. (5) is valid for thin films when unsteady thermal effects are negligible (such that the film is at thermal equilibrium in the vertical direction at all times) if the product of the square of the dimensionless film thickness $\epsilon^{2}$ with the Péclet number $\mbox{{Pe}}=UL/\kappa$ (with $\kappa$ the thermal diffusivity) is small Sansom04 ; we assume that this condition is satisfied.

The system depends on time via two effects. First, as the substrate rotates, the orientation of the gravity vector changes with time. We represent the substrate kinematics by the two angles $\theta(t)$ and $\beta(t)$ such that $g_{x}=g\sin\theta$, $g_{y}=g\sin\beta$, and $g_{z}=g\cos\theta\cos\beta$. Second, (5) is complemented with time-dependent boundary conditions: while the substrate is assumed to be at a fixed temperature $T_{s}$, the upper surface of the liquid layer is subject to convective heat transfer with the surrounding atmosphere at temperature $T_{\infty}$ (valid for small Biot numbers Oron1997 , which is satisfied for almost all the conditions considered in this study),

[TABLE]

where $h_{c}$ is the convective heat transfer coefficient and $k$ the liquid thermal conductivity.

As a result, the temperature-dependent viscosity of the liquid is also time-dependent. We build on the analysis of Bercovici94 ; Sansom04 to account for the temperature-dependence of the viscosity in the lubrication approximation framework, and assume an exponential dependence of viscosity on temperature:

[TABLE]

The heat equation (5) subject to the boundary conditions (6) can readily be integrated to yield an explicit relation between the temperature profile across the liquid layer $T(x,y,z,t)$ and its thickness $h(x,y,t)$. Accordingly

[TABLE]

This equation shows that the temperature distribution is “enslaved” to the film thickness distribution and the thermo-physical properties of the liquid. Integrating equations (1)-(2) with respect to $z$, subject to the no-slip boundary condition ($\mathbf{v}=\mathbf{0}$) at the substrate and to the stress-free boundary condition at the free surface ($-p\mathbf{n}_{s}+\mu\nabla\mathbf{v}\cdot\mathbf{n}_{s}=\mathbf{0}$, with $\mathbf{n}_{s}$ the unit vector normal to the surface, and neglecting capillary effects) yields the following expressions for the velocity:

[TABLE]

The conservation of mass requires that

[TABLE]

where the discharge $\mathbf{q}=\int\limits_{0}^{h}\mathbf{v}\,\text{d}z$ is given by

[TABLE]

This second-order, non-linear, parabolic partial differential equation forms the basis of the forthcoming analysis. The computational domain $D$ is a disk of radius $R$ and the corresponding boundary is denoted by $\partial D$ with unit outward-pointing normal vector $\mathbf{n}_{D}$. On the boundary, a no-flux condition is imposed such that $\mathbf{q}\cdot\mathbf{n}_{D}=0$ on $\partial D$.

II.2 Solution procedure

Computing the evolution of the film thickness $h(x,y,t)$ consists in solving Eq. (11) for a prescribed substrate kinematics (prescribed $\beta(t)$ and $\theta(t)$) and given initial/boundary conditions.) We call this problem the “direct” problem, as opposed to the “adjoint” problem involved in the optimal control (section IV). In order to solve this non-linear partial differential equations, we use the COMSOL Multiphysics software which offers an integrated Finite Element environment to solve a wide range of PDEs. We use the Coefficient Form PDE module which first requires casting the PDE in the following standard form

[TABLE]

where the coefficients $e_{a}$, $d_{a}$, $c$, $\mathbf{\alpha}$, $\mathbf{\beta}$, $a$, $\mathbf{\gamma}$, $f$ to be specified are identified from (11). The computational domain $D$ is meshed with approximately 9,200 Lagrange quadratic elements. The variable-order, variable-step-size backward differentiation formula (BDF) is used for time integration and the time step is limited to $10^{-2}$ s. These settings ensure a robust, mesh-independent solution.

III Direct problem results

The goal of the present work is to identify the optimal substrate kinematics to achieve a uniform coverage by the liquid layer of the substrate. Because the substrate is a disk of radius $R$, the optimal layer thickness is $h_{\text{opt}}=V_{D}/(\pi R^{2})$ for a given volume of liquid $V_{D}$ spread uniformly over the whole disk. Therefore, a good measure of thickness uniformity is

[TABLE]

A uniform film yields $\mathcal{U}=0$, while large deviations from $h_{\text{opt}}$ yield large values of $\mathcal{U}$. The question therefore arises as to how the various physical parameters involved in the problem affect the film thickness uniformity. The next two subsections illustrate:

• •

first, the effect of the heat transfer and the resulting viscosity variation on the levelling dynamics when the substrate is held fixed and horizontal;

• •

second, the effect of the substrate kinematics on the film thickness uniformity for prescribed heat transfer conditions.

The parameter space is vast, but to fix ideas and for illustration purposes, we set the following parameters inspired by the process of crêpe making (see also Table 1).

The viscosity parameters $\mu_{0}$ and $\alpha$ are chosen such that $\mu(T=200^{o}$C$)=5\times 10^{-3}$ Pa s and $\mu(T=20^{o}$C$)=10$ Pa s. The thermo-physical properties are such that the liquid layer cools down from a temperature close to the surrounding temperature $T_{\infty}$ at $t=0$, down to the substrate temperature $T_{s}$ as $t\rightarrow t_{f}$. The volume is chosen such that $h_{\text{opt}}=2\times 10^{-3}$ m. The initial condition is a liquid column of radius $R_{i}=R/3$ and thickness $h_{i}$ required to match the volume constraint. This liquid column rests on a precursor film of thickness $h_{i}/100$ to regularize the degenerate PDE, and the initial thickness $h(x,y,t=0)$ is smoothed over a narrow band such that the profile is twice differentiable to ease convergence without the need for an exceedingly fine mesh.

As mentioned in Sec. II, the centrifugal and Coriolis forces can be neglected if the characteristic rotation rate $\Omega$ is small compared to $\Omega_{1}=\sqrt{\mu U/(\rho LH^{2})}$ and $\Omega_{2}=\mu/(\rho H^{2})$, respectively. We note that, because the film thickness, fluid velocity and viscosity vary in space and time, so do $\Omega_{1}$ and $\Omega_{2}$. Given the above choice of parameters, typical values may be taken in the range $H\sim[2,5]\times 10^{-3}$ m, $U\sim[2,5]\times 10^{-2}$ m s*-1*, and $\mu\sim[0.1,1]$ Pa s (corresponding to $T\sim[20,100]^{o}$C); the typical lateral extension is taken as $L\sim 0.15$ m. This yields $\Omega_{1}\sim[0.8,9]$ rad s*-1* and $\Omega_{2}\sim[5,250]$ rad s*-1*. In the following, the control will impose rotation rates of the order of $\Omega\sim 0.5$ rad s*-1* (see e.g. Figs. 5 and 9), generally smaller than both $\Omega_{1}$ and $\Omega_{2}$. Therefore, at the exception of short times when $\mu$ and $U$ are small and $H$ is large, the centrifugal and Coriolis forces can be reasonably neglected.

III.1 Effect of heat transfer

First, the substrate is held fixed and horizontal. The fluid layer, initially at high-temperature, has a low viscosity and therefore spreads well on the substrate. As the fluid spreads and cools down, its viscosity increases, resulting in a slowdown of the spreading. From this reasoning, it is clear that the cooling rate plays a major role on the spreading dynamics. Two processes are at play: (i) heat convection with the surrounding atmosphere at high temperature $T_{\infty}$, and (ii) cooling from the substrate at low temperature $T_{s}$. The weaker the heat convection, the faster the cooling. Therefore, as $h_{c}$ is reduced, cooling occurs faster, slowing the spreading and reducing the substrate coverage.

Fig. 2 shows the evolution of the free surface for a high value of $h_{c}$ for which cooling is slow (Fig. 2) and a low value of $h_{c}$ for which cooling is fast (Fig. 2). The advancing front is clearly shown to progress faster for the higher value of $h_{c}$ reaching the domain boundary at the latest times. Consequently, the leveling is better for that case, as one would expect.

Fig. 3 confirms that the free surface temperature calculated according to Eq. (8) decreases faster for $h_{c}=3$ W m*-2* K*-1* than for $h_{c}=300$ W m*-2* K*-1*, and therefore the film temperature gets closer to the surface temperature $T_{s}$. This faster cooling corresponds to a faster increase of the viscosity.

The effect of $h_{c}$ on the layer uniformity (14) is shown in Fig. 4. At $t=30$ s, $\mathcal{U}(t)$ has decreased down to $7.1\times 10^{-8}$ m4 for $h_{c}=3$ W m*-2* K*-1*, and to $2.3\times 10^{-8}$ m4 for $h_{c}=300$ W m*-2* K*-1*. This confirms the slower leveling for smaller values of the convective heat transfer coefficient.

In the remainder of this study, the heat transfer coefficient is set to $h_{c}=3$ W m*-2* K*-1*. This comparatively low value, associated with a fast cooling and poor leveling by gravity when the substrate is fixed, calls for a suitable substrate motion in order to spread the film more uniformly.

III.2 Effect of substrate kinematics

We now investigate the effect of the substrate motion. Without any a priori knowledge of what kind of kinematics can efficiently spread the film, we choose two families of time-harmonic angle laws,

[TABLE]

and

[TABLE]

and perform a broad exploration of the space parameter $(A_{1},A_{2},T_{1},T_{2})$. A Monte-Carlo algorithm picks random combinations of amplitudes $A_{1}$, $A_{2}$ defined in $[5,45]^{o}$ and periods $T_{1}$, $T_{2}$ defined in $[5,50]$ s; for each combination, it then solves the governing equations over $t\in[0,30]$ s and evaluates the final uniformity $\mathcal{U}(t_{f})$. The best results obtained after 1000 evaluations with (15) and 1000 evaluations with (16) are shown in Fig. 5.

From the set of Monte-Carlo realizations, the best harmonic kinematics minimizing $\mathcal{U}(t_{f})$ corresponds to Eq. (15) with $A_{1}$=17.9o, $A_{2}$=22.3o, $T_{1}$=8.26 s, $T_{2}$=16.6 s, see Fig. 5. For comparison, the second best kinematics obtained with $A_{1}$=7.70o, $A_{2}$=8.65o, $T_{1}$=6.91 s, $T_{2}$=33.3 s in Eq. (16) is shown in Fig. 5. The value of $\mathcal{U}(t_{f})$ is $4.19\times 10^{-8}$ m4 and $5.69\times 10^{-8}$ m4 for the best and second best control, respectively. This moderate improvement in film uniformity is confirmed on Fig. 4 showing the time evolution of $\mathcal{U}(t)$. For both cases, we note that the amplitude of the variations in $\theta$ and $\beta$ are of the same order. Moreover, both signals are out of phase and the longest period is an approximate multiple of the shortest one (the period of $\beta$ is approximately twice that of $\theta$ for the best harmonic control, and four times for the second best). In the context of crêpe making, this kinematic corresponds to a slow rocking motion in one direction composed with a transverse rocking motion at a higher frequency.

The corresponding evolution of the film thickness distribution for the best and second best harmonic controls is illustrated in Figs. 6 and 7, respectively. In order to better grasp the substrate kinematics on these figures, the orientation of the substrate is indicated by a vector in the direction of the projection of the gravity vector in the plane of the substrate. Videos are also available as supplementary material SM . In both figures, a positive $\beta$ corresponds to an upwards pointing arrow and a positive $\theta$ to a rightwards pointing arrow. These figures and videos suggest that a good strategy involves draining the fluid to one end of the surface (low-frequency rocking around the first axis) and redistributing it by a sideways tilt of the surface (higher-frequency rocking around a second, perpendicular, axis).

IV Optimal control formulation

The previous section illustrated the spreading behavior of the fluid on a substrate at rest (sec. III.1), and on a substrate moving with time-harmonic kinematics (sec. III.2). We now turn our attention to the question of finding better kinematics in order to improve the uniformity of the film. In this section we present the adjoint framework used to modify the kinematics so as to improve uniformity iteratively, in an efficient way.

IV.1 Optimization problem

Let us define our objective: we wish to find a control law $\mathbf{c}(t)=(\theta(t),\beta(t))$ over the time interval $I=[0,t_{f}]$ that minimizes a measure of the difference between the final film thickness and the ideal uniform thickness $h_{\text{opt}}=V_{D}/(\pi R^{2})$. We may wish to achieve this at minimal control cost, i.e. by spending as little energy as possible to move the substrate. Therefore, we introduce a scalar objective function that includes both the final uniformity and the total cost of the control

[TABLE]

where the weight $\gamma\geq 0$ can be seen as the unit cost of the control and allows one to choose how much the control should be penalized. The subscript t in ${\mathcal{J}}_{t}$ stands for “terminal control” because the uniformity is evaluated at the final time $t_{f}$ only (see Section IV.4 for “regulation control”, where the uniformity is evaluated at all times). The optimization problem reads

[TABLE]

together with the initial and boundary conditions specified in sections II-III.

This constrained optimization problem can be solved with gradient-based methods, which requires computing the gradient $\mathrm{d}{\mathcal{J}}_{t}/\mathrm{d}\mathbf{c}$ of the objective function with respect to the control. Given the infinite number of degrees of freedom in the control $\mathbf{c}(t)$ (which becomes a finite but large number after discretization), evaluating the gradient numerically with a finite-difference approach (where the value of ${\mathcal{J}}_{t}$ is evaluated repeatedly for each degree of freedom perturbed once at a time) has a prohibitive computational cost. Instead, as classically done in optimal control problems, one can use an adjoint-based approach and compute the gradient very efficiently.

IV.2 Adjoint equation

Computing the gradient $\mathrm{d}{\mathcal{J}}_{t}/\mathrm{d}\mathbf{c}$ efficiently relies on solving an adjoint equation. Let us detail now the derivation of this equation. We transform the constrained problem (18) into an unconstrained problem by introducing the Lagrangian

[TABLE]

where the adjoint variable ${h^{\dagger}}(\mathbf{x},t)$ is a Lagrange multiplier that enforces the governing equation (11). In the following, derivatives on functional spaces are understood as Fréchet derivatives, which we denote loosely for any scalar or vectorial quantity $s$:

[TABLE]

The total derivative of the Lagrangian with respect to the control reads

[TABLE]

and it reduces to

[TABLE]

when both $\partial{\mathcal{L}}_{t}/\partial h=0$ and $\partial{\mathcal{L}}_{t}/\partial{h^{\dagger}}=0$. If, in addition, the governing equation (11) is satisfied, then ${\mathcal{L}}_{t}={\mathcal{J}}_{t}$ by construction, and thus the gradient of interest is obtained as

[TABLE]

Specifically, the gradient reads

[TABLE]

It has two time-dependent components, that will be used to update the control $(\theta(t),\beta(t))$.

The yet unknown adjoint variable ${h^{\dagger}}(\mathbf{x},t)$ remains to be specified via the condition $\partial{\mathcal{L}}_{t}/\partial h=0$, which is discussed below. Note that the other condition to be satisfied, $\partial{\mathcal{L}}_{t}/\partial{h^{\dagger}}=0$, is by construction the governing equation (11).

At this point, we note that one can simplify the problem to a large extent when $\alpha A^{\prime}\ll 1$, which is well verified in practice for our choice of parameters $\alpha$, $k$, $h_{c}$ and $T_{\infty}-T_{s}$. In this case, the discharge (12) is well approximated by

[TABLE]

This expression will be used to derive an approximate adjoint equation, and therefore an approximate gradient. Note however that the direct equation for $h(\mathbf{x},t)$ is solved in its full form (11). For simplicity, in the following we denote $W=e^{\alpha T_{s}}\rho g/(3\mu_{0})$, such that the vector in the second term of (27) reads

[TABLE]

The condition $\partial{\mathcal{L}}_{t}/\partial h=0$ yields, after integration by parts, the adjoint equation associated with this terminal control problem,

[TABLE]

together with terminal and boundary conditions for ${h^{\dagger}}(\mathbf{x},t)$:

[TABLE]

Given the terminal condition (40), the adjoint equation (39) must be solved backward in time from $t=t_{f}$ to $t=0$. Note that the adjoint equation is linear in ${h^{\dagger}}(\mathbf{x},t)$, and its coefficients depend on the direct solution $h(\mathbf{x},t)$.

IV.3 Solution procedure

The iterative optimization procedure is the following:

Given a tentative control $\mathbf{c}_{k}(t)=(\theta_{k}(t),\beta_{k}(t))$, solve the governing equation (11) $\rightarrow$ obtain the solution $h(\mathbf{x},t)$; 2. 2.

Given the control $\mathbf{c}_{k}(t)$, and the solution $h(\mathbf{x},t)$ from step (1), solve the adjoint equation (39) $\rightarrow$ obtain the adjoint solution ${h^{\dagger}}(\mathbf{x},t)$; 3. 3.

Given the control $\mathbf{c}(t)$, the solution $h(\mathbf{x},t)$ from step (1) and the adjoint solution ${h^{\dagger}}(\mathbf{x},t)$ from step (2), evaluate the time-dependent gradient $\mathrm{d}{\mathcal{J}}_{t}/\mathrm{d}\mathbf{c}$ with (27); 4. 4.

Given the gradient $\mathrm{d}{\mathcal{J}}_{t}/\mathrm{d}\mathbf{c}$ from step (3), update the control $\mathbf{c}_{k}(t)\to\mathbf{c}_{k+1}(t)$ with a gradient-based method. Go back to step (1) and iterate until convergence.

A few remarks are in order. At the first iteration $k=1$, we initialize the control $(\theta(t),\beta(t))$ either with an arbitrary guess, or with one of the promising harmonic control laws obtained with the Monte-Carlo method. The choice of the initial control can have a substantial effect on the outcome of the optimization because ${\mathcal{J}}_{t}(\mathbf{c})$ is not convex and may exhibit local minima, where gradients methods can stay trapped. To increase the chances of achieving a satisfactory reduction of ${\mathcal{J}}_{t}$, we repeat the optimization from different initial guesses.

At each iteration, the full solution $h(\mathbf{x},t)$ is saved in step (1) and later used to solve the adjoint equation in step (2). Given the reasonable size of the solution (number of mesh points and number of time steps), we do not need to use intermediate checkpoints.

We update the control in step (4) with the Polak-Ribiere variant of the conjugate gradient method (for the choice of the descent direction) and Brent’s method (for the choice of the descent distance in that direction).

Finally, the above steps are repeated until the relative variations of the two terms in ${\mathcal{J}}_{t}$ are both smaller than $10^{-6}$. We observe that when this convergence criterion is satisfied, the norm of the gradient $\mathrm{d}{\mathcal{J}}_{t}/\mathrm{d}\mathbf{c}$ is generally small too.

IV.4 Regulation control

The objective function (17) measures uniformity at the final time $t_{f}$ only (“terminal control”). It is possible to use an alternative measure of uniformity that is distributed over the whole time interval $[0,t_{f}]$ (“regulation control”):

[TABLE]

Which formulation eventually leads to the best uniformity at $t=t_{f}$ is not obvious a priori. It may be expected that terminal control, which allows poorer short-term performance provided it yields a long-term advantage, is less restrictive and therefore more effective Bewley2001 . In this study we focus on terminal control, but for the sake of completeness we nonetheless consider the optimization problem

[TABLE]

with the same initial and boundary conditions as for (18). The Lagrangian is

[TABLE]

and the adjoint equation associated with this regulation control problem reads:

[TABLE]

with terminal and boundary conditions

[TABLE]

In contrast to the terminal control problem, where the adjoint equation (39) is homogeneous and the adjoint dynamics are determined by the non-zero terminal condition (40), in the regulation control problem the terminal condition (50) is zero and the adjoint dynamics are determined by the time-dependent forcing on the right-hand side of (49). Finally, the gradient needed to update the control so as to decrease ${\mathcal{J}}_{r}$ reads

[TABLE]

V Optimal control results

Results of the optimization are reported in Fig. 8 for a wide range of unit cost values $\gamma\in[10^{-8},10^{-5}]$. Fig. 8(a) shows the objective function ${\mathcal{J}}_{t}$, and Fig. 8(b) the final uniformity measure $\mathcal{U}(t_{f})$. (In this section, we give numerical values of $\mathcal{U}(t_{f})$ in m4 and omit units.) Recall that, as defined in (17), small $\gamma$ values correspond to cheap controls (large angles $|\theta|$, $|\beta|$ are not penalized as they barely affect ${\mathcal{J}}_{t}$), whereas large $\gamma$ values correspond to expensive controls (large angles $|\theta|$, $|\beta|$ are penalized as they contribute to increasing ${\mathcal{J}}_{t}$).

The uncontrolled case (fixed horizontal substrate, $\theta=\beta=0$) leads at $t=t_{f}$ to a uniformity $\mathcal{U}(t_{f})=6.80\times 10^{-8}$ (solid line; see also section III.1). Regarding time-harmonic control obtained with the Monte-Carlo algorithm (dashed line with open symbols; see also section III.2), the best control for ${\mathcal{J}}_{t}$ yields $\mathcal{U}(t_{f})=4.35\times 10^{-8}$ for small $\gamma$ values. As $\gamma$ increases, the second-best control in terms of uniformity ($\mathcal{U}(t_{f})=5.52\times 10^{-8}$) becomes the best control for ${\mathcal{J}}_{t}$ thanks to its smaller angles and lower total cost. However, although this time-harmonic control yields a slightly better uniformity $\mathcal{U}(t_{f})$ than the uncontrolled case, it eventually becomes less efficient for ${\mathcal{J}}_{t}$ when $\gamma\gtrsim 5\times 10^{-7}$.

The adjoint optimization (filled symbols) yields a substantial improvement. For small $\gamma$, the uniformity is improved by a factor 6, down to $\mathcal{U}(t_{f})=1.16\times 10^{-8}$. As the unit cost increases, ${\mathcal{J}}_{t}$ and $\mathcal{U}(t_{f})$ increase but remain smaller than their uncontrolled counterparts up to $\gamma$ as large as $10^{-5}$.

Figure 9 shows the best control laws $\theta(t)$, $\beta(t)$ obtained with the adjoint optimization. Clearly, a larger unit cost leads to smaller angles: while inclinations as large as $45^{o}$ are admissible for $\gamma=10^{-8}$, they do not exceed $10^{o}$ for $\gamma=10^{-6},$ and barely differ from the no-control case $\theta=\beta=0$ for $\gamma=10^{-5}$. The kinematics displayed here is close to a damped rotation motion: $\theta(t)$ and $\beta(t)$ are approximately periodic, decreasing in amplitude, and in phase quadrature (azimuthally traveling wave). About two rotations are described for $\gamma\leq 10^{-7}$ and one rotation for $\gamma\leq 10^{-6}$, possibly because the larger angles allowed for small $\gamma$ can spread the film faster, which in turn allows for faster kinematics.

This is confirmed by the time evolution of the film thickness in Figs. 10-11 (see also videos in the supplementary material SM ). For $\gamma=10^{-8}$ (Fig. 10), first the bulk of the liquid is quickly moved to the rim of the disk, leaving behind a thinner film in the central region of the disk ($0\leq t\leq\Delta t$); the thickest portion of the film is then displaced along the whole rim, thus depositing some liquid in the yet uncovered disk areas (first rotation, $\Delta t\leq t\leq 4\Delta t$); finally, another rotation further distributes the liquid and improves the thickness uniformity, albeit less markedly because the film has already become more viscous. For $\gamma=10^{-6}$ (Fig. 11), smaller inclinations cause the fluid to move more slowly. Only about one rotation is completed by the time viscosity becomes so large that the thickness uniformity cannot be modified substantially. However, the optimized control manages to distribute the liquid just about everywhere over the disk. Remarkably, this is achieved by gradually slowing down the rotating motion as the film becomes more viscous (compare $0\leq t\leq 2\Delta t$ and $2\Delta t\leq t\leq 6\Delta t$ in Figs. 9(c) and 11).

For a better visualization of the final state, Fig. 12 shows contours of film thickness in regions where it is within 20% of $h_{\text{opt}}$. Without control, $h(\mathbf{x},t_{f})$ falls in this interval only in a narrow ring: the inner region of the substrate remains much thicker than $h_{\text{opt}}$, and the outer region remains much thinner, resulting in a highly non-uniform film. Control improves the film spreading, albeit in a non-axisymmetric fashion. The best harmonic control, expensive adjoint control, and cheap adjoint control result in an increasingly wider area where the film thickness departs from $h_{\text{opt}}$ by less than 20% (see also Figs. 6 and 10-11).

The evolution of $\mathcal{U}(t)$ in Fig. 13 shows interesting features. First, as already mentioned, the adjoint optimization improves the final thickness uniformity compared to the uncontrolled case, and more substantially as $\gamma$ decreases. Second, our terminal control formulation allows $\mathcal{U}(t)$ to increase temporarily because this yields a lower final value $\mathcal{U}(t_{f})$. In the specific case at hand ($\gamma=10^{-8}$), a sharp increase around $t\simeq 5$ s brings $\mathcal{U}$ back to its initial (and maximal) value. With less time available and a liquid that has become more viscous, the situation may seem compromised. However, from this point onward, $\mathcal{U}$ decreases quickly and steadily, finally yielding the best uniformity achieved in the present study.

One may wonder whether more uniform films are also smoother. While of course this is not true in general, it is worth investigating the smoothness of the final films obtained with and without control. A simple measure of smoothness is based on the gradient of the thickness:

[TABLE]

A uniform film yields $\mathcal{S}=0$, while large gradients $\bm{\nabla}h$ yield large values of $\mathcal{S}$. Other measures exist, putting for instance different weights on small or large spatial wavelengths Bewley2001 ; Foures2014 , but (55) suffices for present purposes. In all cases, we observe that $\mathcal{S}(t_{f})<\mathcal{S}(0)$: overall, the film becomes smoother as time evolves. This is consistent with the initial film exhibiting a localized but very sharp front at $r=R_{i}$, whereas final films have smoother fronts due to the effect of gravity. We also note that film uniformity and smoothness are indeed correlated. Table 2 summarizes final values of $\mathcal{U}$ and $\mathcal{S}$. Both measures yield similar trends in terms of control performance: cheaper adjoint control, followed by more expensive adjoint control, best harmonic controls, and no control. This trend is also illustrated by the final gradient fields in Fig. 14. This suggests that good results in terms of uniformity $\mathcal{U}$ could be obtained by formulating the optimization with the smoothness measure $\mathcal{S}$ in the objective function. In some cases, it so happens that optimization problems formulated with an objective function different from the true objective yield better result (see Bewley2001 for an example where drag reduction in a plane channel flow is best achieved by targeting kinetic turbulent energy).

Interestingly, we have obtained other controls, with qualitatively different kinematics but similar values of ${\mathcal{J}}_{t}$ and $\mathcal{U}(t_{f})$. This suggests that the objective function has several local minima, as mentioned in section IV.3. We cannot rule out the possibility that better controls exist. However, running the optimization algorithm from a set of initial guesses did allow us to decrease $\mathcal{U}(t_{f})$ substantially compared to both the uncontrolled case and the best time-harmonic Monte-Carlo controls, with a set of physical parameters constituting a challenging test case, as cooling occurs quickly and the fluid becomes more viscous.

The adjoint optimization algorithm generally converged within 10 to 30 gradient evaluations, each iteration consisting of 1 adjoint calculation (to determine the descent direction) and 5 to 10 direct calculations (to determine the descent distance in that direction). This is to be contrasted with the 2000 evaluations and poorer results of the Monte-Carlo algorithm. (To be fair, it must be recalled that the Monte-Carlo algorithm was restricted to harmonic controls, that are necessarily sub-optimal compared to arbitrary controls.) It is precisely the very large number of degrees of freedom, optimized efficiently within a reasonably small number of iterations, that makes adjoint methods powerful and attractive for optimal control in this and other settings.

VI Conclusions

The question at the center of this study is: can one improve the coverage of a surface by a gravity-driven liquid film using a suitable kinematics of the surface to distribute the liquid uniformly? This question is answered in the framework of PDE-constrained optimization, for which ones seeks to minimize an objective function (a combination of the film thickness uniformity and of the cost of moving the substrate) subject to a set of constraints (governing PDE with associated boundary and initial conditions). The PDE that governs the film dynamics is obtained in the lubrication approximation framework, using a temperature-dependent viscosity and a gravity force whose direction and magnitude depend on the orientation of the surface. This effectively expresses the surface kinematics as a time-dependent body force in the governing equations. This model has inherent limitations, namely that the flow must be free of inertia, and the free surface slope must be small. Both assumptions appear to be reasonable for the applications of interest here: pancake making and surface coating.

Two methods are proposed to optimize the film uniformity. The first one assumes a kinematics described by harmonic functions, and identifies the optimal amplitudes and periods from a set of randomly distributed Monte-Carlo realizations. This method is computationally costly as it blindly samples the large parameter space, and it is shown to offer some improvement in thickness uniformity (40% improvement compared to the uncontrolled case). The second method is a gradient-based method that allows for arbitrary controls. We have shown that the gradient of the objective function can be conveniently calculated using an adjoint formulation. This second approach is shown to be much more effective as it leads to a 83% improvement of the uniformity measure for the cheapest control ($\gamma=10^{-8}$). The computational cost of this second method is lower than that of the Monte-Carlo method since local information is available to guide the optimization in a suitable descent direction. Interestingly, one of the optimal controls identified by the adjoint method appears to replicate the motion one would naturally adopt when making crêpes, i.e. first draining all the batter to one end of the pan, and then slowly rotating it for one revolution or two in order to distribute the batter in the remaining part of the pan. Results also suggest the existence of multiple local minima since different kinematics can lead to commensurate uniformity. The framework presented here can be adapted to other optimal control problems involving thin liquid films.

VII Acknowledgments

The authors wish to express their gratitude to D. Sellier for carefully proofreading the manuscript. The authors also thank the financial support of the University of Canterbury through the College of Engineering Strategic Research grant and the Innovation Jumpstart grant.