1 Introduction
Unstable thin viscous liquid films coating a vertical fibre display complex interfacial dynamics and various flow regimes. Driven by Rayleigh–Plateau instability and gravity effects, a rich variety of dynamics can occur including the formation of droplets, regular travelling wave patterns and irregular droplet coalescence. Such dynamics has attracted a lot of attention from researchers in recent years due to its widespread applications in heat and mass exchangers [Reference Zeng, Sadeghpour and Ju28], desalination [Reference Sadeghpour, Zeng, Ji, Ebrahimi, Bertozzi and Ju22,Reference Zeng, Sadeghpour and Ju29] and particle capturing systems [Reference Sadeghpour, Oroumiyeh, Zhu, Ko, Ji, Bertozzi and Ju1].
With proper choices of flow rate, liquid, fibre radius and inlet geometry, three typical flow regimes are observed in experiments [Reference Ji, Sadeghpour, Ju and Bertozzi12,Reference Kalliadasis and Chang15,Reference Sadeghpour, Zeng and Ju23]: a convective instability regime where bead coalescence happens repeatedly, a travelling wave regime where a steady train of beads flow down the fibre at a constant speed, and an isolated droplet regime where widely spaced large droplets are separated by small wave patterns. When other system parameters are fixed, varying the flow rate from high to low can lead to a flow regime transition from the convective to the travelling wave regime, and eventually to the isolated droplet regime. A better understanding of the travelling wave pattern is expected to provide insights for many engineering applications that utilise steady trains of liquid beads.
At small flow rates, the dynamics of the axisymmetric flow on a cylinder is typically modelled by classical lubrication theory. Under the assumption that the film thickness is much smaller than the radius of the cylinder, Frenkel [Reference Frenkel9] proposed a weakly nonlinear thin-film equation for the film thickness h (or the height of the film) that captures both stabilising and destabilising effects of the surface tension in the dynamics. This evolution equation was further investigated by Kalliadasis & Chang [Reference Kalliadasis and Chang15], Chang & Demekhin [Reference Chang and Demekhin4], and Marzuola et al. [Reference Marzuola, Swygert and Taranets18]. Similar models for weakly rippled thin films were also studied in [Reference Rosenau and Oron19,Reference Shlang and Sivashinsky24], and the paper of Rosenau & Oron [Reference Rosenau and Oron19] incorporates fully nonlinear curvature terms that account for the deformation of the film interface based on asymptotic analysis. To relax the thin-film assumption, Craster & Matar [Reference Craster and Matar6] developed an asymptotic model which assumes that the film thickness is much smaller than the capillary length. In 2000, Kliakhandler et al. [Reference Kliakhandler, Davis and Bankoff17] extended the thin-film model to consider a thick layer of viscous fluid by introducing fully nonlinear curvature terms, which leads to an evolution equation (1.1) for the film thickness h(x,t),
where $\sigma > 0$ is the Bond number, $r_0 > 0$ is the dimensionless fibre radius and the mobility Q(h) takes the form
While equation (1.1) is a model equation that was not rigorously derived asymptotically, it constitutes all the necessary terms to describe the corresponding physical process. The term $[Q(h)]_x$ in (1.1) represents gravitational effects, ${h_{xx}}/{(1+ h_x^2)^{3/2}}$ and ${(h+r_0)^{-1}(1+ h_x^2)^{-1/2}}$ describe the stabilising and destabilising roles of the surface tension due to axial and azimuthal curvatures of the interface, respectively.
In 2019, Ji et al. [Reference Ji, Falcon, Sadeghpour, Zeng, Ju and Bertozzi10] studied a family of full lubrication models that incorporate slip boundary conditions, fully nonlinear curvature terms and a film stabilisation mechanism. The film stabilization term,
is motivated by the form of disjoining pressure widely used in lubrication equations [Reference de Gennes7,Reference Thiele26] to describe the wetting behaviour of a liquid on a solid substrate, and the scaling parameter $A > 0$ is typically selected based on a stable liquid layer in the coating film dynamics. Numerical investigations against experimental results in [Reference Ji, Falcon, Sadeghpour, Zeng, Ju and Bertozzi10] show that compared with previous studies, the combined physical effects better describe the propagation speed and the stability transition of the moving droplets. For higher flow rates, coupled evolution equations of both the film thickness and local flow rate are developed [Reference Ruyer-Quil, Treveleyan, Giorgiutti-Dauphiné, Duprat and Kalliadasis20,Reference Ruyer-Quil, Trevelyan, Giorgiutti-Dauphiné, Duprat and Kalliadasis21,Reference Trifonov27]. These equations incorporate inertia effects and streamwise viscous diffusion based on the integral boundary-layer approach. Recently, Ji et al. [Reference Ji, Sadeghpour, Ju and Bertozzi12] further extend a weighted-residual integral boundary-layer model to incorporate the film stabilisation mechanism to address the effects of the inlet nozzle geometry on the downstream flow dynamics.
While these previous studies primarily focus on the modelling of coating film dynamics, theoretical analysis of these models is still lacking. In this paper, we focus on the model (1.1) with an additional generalised film stabilisation term motivated by [Reference Ji, Falcon, Sadeghpour, Zeng, Ju and Bertozzi10,Reference Ji, Sadeghpour, Ju and Bertozzi12]. Substituting $u = h + r_0$ into (1.1) and including the generalised film stabilisation term $\tilde{\Pi}(u) = -A/u^m$ , we obtain a fourth-order nonlinear partial differential equation for u(t, x), namely,
where the scaling parameter $A \geqslant 0$ , the exponent $m > 0$ , and the mobility becomes
A useful observation is that
In this paper, we will show the existence of non-negative weak solutions to (1.4) with a focus on the travelling wave solutions. Under proper conditions, u(x,t) converges to a travelling wave solution in long time.
The structure of the rest of the paper is as follows. In Section 2, we present the main theorem on the existence of non-negative weak solutions to the problem, followed by the proof of the existence theorem in Section 3. Section 4 focuses on the travelling wave solutions of the problem. Numerical studies based on the analytical results are included in Section 5, followed by concluding remarks in Section 6.
2 Existence of non-negative weak solutions
In this section, we will define a generalised non-negative weak solution and formulate an existence of the solution statement for the following problem:
where $ \Omega \subset \mathbb{R}^1$ is bounded domain, $m > 0$ , $A \geqslant 0$ , $Q_T = \Omega \times (0,T)$ , $ T > 0$ , $Q(r_0) = 0$ , $f(z) := (1 + z^2)^{-\frac{1}{2}}$ , and
Assume that
Integrating (2.1) in $\Omega\times (0,t)$ , by (2.2) we have
Let us denote by
Here, $-J$ represents a dynamic pressure [Reference Ji, Falcon, Sadeghpour, Zeng, Ju and Bertozzi10,Reference Ji and Witelski13,Reference Ji and Witelski14] that incorporates both axial and azimuthal surface tension effects and the film stabilisation term.
Definition 2.1 A generalised weak solution of the problem (2.1)–(2.3) is a non-negative function $u(x,t)$ satisfying
where $\mathcal{P}_T = \overline{Q}_T \setminus ( \{u= r_0\} \cup\{t=0\})$ and $u$ satisfies (2.1) in the following weak sense:
for all $ \phi \in L^2(0,T; H^1(\Omega))$ and $\psi \in L^2(Q_T)$ ;
Let us denote by $\tilde{G}_0^{(\alpha)}(z)$ the following function
where $s_0$ is a positive constant. The function $\tilde{G}_0^{(\alpha)}(z)$ is a generalisation of the entropy introduced in [Reference Beretta, Bertsch and Dal Passo2].
Theorem 1 Assume that the initial function $u_0$ satisfies (2.4). Then problem (2.1)–(2.3) has a weak solution $u(x,t)$ defined in $Q_{T}$ for any $T > 0$ , in the sense of Definition 2.1. Assume that the initial function $u_0 \geqslant r_0 > 0$ also satisfies
Then the solution $u(x,t)$ satisfies $ u \geqslant r_0 >0 $ and
3 Proof of Theorem 1
3.1 Regularised problem
Note that (2.1) is a parabolic equation which degenerates when $u = 0$ , $ u = r_0$ and $|u_x| \to \infty$ . Therefore, we have to apply a regularisation technique to this equation to overcome the degeneracies. For given $\varepsilon > 0$ and $\delta > 0$ , we have
where
where $\gamma \in (0,1)$ , $\beta >6$ , $\theta \in (0,\frac{1}{3})$ , and $\kappa \in (0,\frac{1}{\beta -2})$ .
The positivity condition $\varepsilon > 0$ eliminates the degeneracies at $u =0,\, r_0$ in (2.1). This condition will ‘lift’ the initial data and smooth the initial data up to $C^{4+\gamma}(\Omega) $ . The $\delta > 0$ condition eliminates the degeneracy at $|u_x| \to \infty$ in (2.1). As a result, the regularised equation (3.1) is uniformly parabolic. Moreover, the boundary conditions (3.3) are of Lopatinskii–Shapiro type that implies the existence of a proper continuation of solutions to the whole line (cf. [Reference Solonnikov25], for example). Using the well-known parabolic Schauder estimates from [Reference Solonnikov25], one can generalise [8, Theorem 6.3, p. 302] and show that the regularised problem has a unique classical solution $u_{\delta \varepsilon} \in {C}_{x,t}^{4+\gamma,1+\gamma/4}( \Omega \times [0, \tau_{\delta \varepsilon}])$ for some time $\tau_{\delta \varepsilon} > 0$ . Here $\tau_{\delta \varepsilon}$ is the local existence time from [8, Theorem 6.3, p. 302].
3.2 A priori estimates
Now, we need to derive a priori estimates for energy-entropy functionals. Let us denote
Multiplying (3.1) by $- J_{\varepsilon \delta} $ and integrating over $\Omega$ , we obtain
As
then we get
where
Integrating (3.1) and taking into account periodic boundary conditions, we get
By (3.5) we deduce that
whence
Taking into account
(3.5) and (3.6), from (3.4) we find that
By (3.7) we obtain
for all $t \geqslant 0$ . Let us denote by
Multiplying (3.1) by $ (G^{(\alpha)}_{\varepsilon }(u))'$ and integrating over $\Omega$ , we obtain
where $(\tilde{G}^{(\alpha)}_{\varepsilon }(z) )' = g_\varepsilon(z) \, (G^{(\alpha)}_{\varepsilon } (z))' .$ By periodic boundary conditions, we have
where $C_{\alpha} = \alpha^2 + | \alpha | +1 $ , whence
From (3.9) we have
Next, we will use the following lemma:
Lemma 3.1 [5, Lemma 3.1, p. 806] Let $v \in H^1(\Omega)$ be a nonnegative function. Then for any $p > 2$ , there exists a constant $\tilde{C} $ depending on p and $\Omega$ such that
Using (3.11) with $v = g^{2}_\varepsilon(u)$ and $p = \frac{\beta-2}{2} > 2$ , i.e. $\beta > 6$ , we have
whence, due to (3.8), we find that
Choosing $\alpha = 1$ and $\varepsilon $ small enough in (3.10), taking into account (3.8) and (3.13), we deduce that
Integrating (3.14) in time, taking into account (3.8), we arrive at
for all $t \geqslant 0$ , where
Let $\tilde{\mathcal{G}}_{\varepsilon } (z) := \tilde{G}^{(1)}_{\varepsilon } (z)$ and $\mathcal{G}_{\varepsilon } (z) := G^{(1)}_{\varepsilon } (z)$ . Note that
provided $|z| \leqslant B $ , where $C = \frac{1}{2} (B +C_0(1+B^4 \log(\frac{B}{r_0})))$ . So,
whence, due to $|Q(z)| \sim | \log ( \frac{z}{r_0} )| $ as $z \to r_0$ , we find that
provided $\theta \in (0,\frac{1}{3})$ . As a result, we deduce that
Therefore, due to (3.16) we have
Note that, after taking limit $\varepsilon \to 0$ , we obtain the limit solution $u_\delta \geqslant r_0 >0$ (the proof is similar to [3, Theorem 4.1]), and for this reason, instead of (3.13), we use $\mathop {\sup} \limits_{\Omega} ( u_{\delta}^{-2} ) \leqslant r_0^{-2}$ for the limit process on $\delta \to 0$ .
3.3 Construction of a weak solution
We will construct a weak non-negative solution using Arzela–Ascoli theorem. By (3.8) we deduce a uniform boundedness of the following sequences
By (3.15) and (3.19) we arrive at
By (3.19) and (3.21) we have (see, e.g. [Reference Bernis and Friedman3])
From (3.17) and (3.28) we obtain that
By (3.29) and (3.27) we find that
By (3.29) and (3.24) we deduce that
By regularity theory of uniformly parabolic equations and by the uniformly Hölder continuity of the $u_{\varepsilon\delta}$ , we deduce that
It follows that
Next, multiplying (3.1) by $\phi$ and integrating in $Q_T$ , we have
for all $ \phi \in L^2(0,T; H^1(\Omega))$ and $\psi \in L^2(Q_T)$ . Due to (3.29)–(3.35), letting $\varepsilon \to 0$ in (3.37) and (3.38), we arrive at
for all $ \phi \in L^2(0,T; H^1(\Omega))$ and $\psi \in L^2(Q_T)$ .
Next, we pass to the limit $\delta \to 0$ in (3.39). Note that a solution $u_{\delta}(x,t) \geqslant r_0$ in $Q_T$ . By (3.17), (3.18) and (3.5) we get
and a. e. in $Q_T$ . By (3.21) and (3.41)
By (3.26), (3.27) and (3.42) we get
Letting $\delta \to 0$ in (3.39) and (3.40), in view of (3.41)–(3.46), we get a solution u(x,t) which satisfies (2.10) and (2.11).
4 Travelling wave solutions
We introduce a change of variables to the reference frame of the travelling wave,
Substituting the change of variables to the PDE (1.4) leads to the following PDE for $v(\xi, t)$
with L-periodic boundary conditions on $\xi$ , where $\mu =0$ corresponds to the case without gravity and $\mu = 1$ corresponds to the case with gravity. Note that the total mass of the film M satisfies the conservation-of-mass condition,
Let us denote
Here, $\mathcal{E}(v)$ is an energy functional that combines the surface tension of the liquid and a local free energy that corresponds to the film stabilisation mechanism, $\lambda$ is a Lagrange multiplier and $\mathcal{L}(v)$ incorporates both $\mathcal{E}(v)$ and the total mass of the film $M = \int_0^L{v^2 d\xi}$ . The functional $-J(v)$ represents dynamic pressure of the liquid, and $\bar{M}$ is the average film thickness.
In the Subsection 4.1–4.3, we consider the case of $\mu =0$ (without gravity), and the case $\mu =1$ (with gravity) is studied in the Subsection 4.4.
4.1 Critical points of the energy functional
Lemma 4.1 Let $m > 2$ and $\mu =0$ . Then, for every $M$ , the functional $\mathcal{E}$ attains its minimum $v_{\min}$ on $\mathcal{M}$ . Moreover, if $0 \leqslant A < r_0^{m-1}$ , then
Proof of Lemma 4.1. First of all, we will show that $ \mathcal{E}(v)$ dissipates. Multiplying (4.1) by $ - J $ and integrating on $\xi$ , we find that
As
then by (4.3)
whence
Moreover, taking into account
we deduce that
The functional $\mathcal{E}(v)$ is non-increasing and bounded from below therefore it attains its minimum $v_{\min}$ on the set $\mathcal{M}$ .
Let us denote by
Note that
whence
As
where
then
On the other hand, we have
So, from (4.5) we obtain
where
Then by (4.4) we arrive at
Compare (4.7) with a solution of the following problem
where $ y(t) := \mathcal{E}(v) > 0$ ,
Then we have
As a result, we obtain that
provided
4.2 Structure of energy functional minimisers
The Euler–Lagrange equation for $\mathcal{E}(v)$ under the constraint $ \int \limits_0^L { v^2 d\xi} = M$ is given by
where $\lambda$ is the Lagrange multiplier. We also need to incorporate the constraint $v \geqslant r_0 > 0$ . If $v \in \mathcal{M}$ , we decompose (0,L) according to the value of v into two sets, P(v) and Z(v).
We compute the first variation of $\mathcal{L}(v)$ about $v_{\min}$ and obtain
where $\phi$ is a smooth test function supported in P(v). Since the first variation of $\mathcal{L}(v)$ along $v_{\min} + \epsilon$ must vanish for every such test function $\phi$ , the Euler–Lagrange equation (4.10) holds.
Similarly, when we take $\phi$ to be a smooth test function supported in Z(v), the first variation of $\mathcal{L}(v)$ along $v_{\min} + \epsilon$ must vanish for every such $\phi$ , and we obtain
Therefore, the energy minimiser $v_{\min}$ satisfies (4.10) on (0,L) in the sense of distributions as $\lambda$ in (4.10) is a piece-wise constant function, namely,
Lemma 4.2 Let $m > 2$ , $\mu =0$ , and $0 \leqslant A < r_0^{m-1}$ . If $v_{\min}$ minimizes $\mathcal{E}$ on $\mathcal{M}$ , then it solves (4.10) with $\lambda(\xi)$ on $(0, L)$ . The Lagrange multiplier $\lambda(\xi)$ is negative and $\lambda^*$ satisfies
where
Furthermore, $v_{\min}$ is of class $C^1$ , and $v'_{\min} =0$ on $ \partial P(v_{\min})$ .
Proof of Lemma 4.2. Next, we consider (4.10) on P(v). By the substitution
we have
On the other hand, if $v' = 0$ , then
The equation (4.12) has the following general solution
where $ C_0 \in \mathbb{R}^1$ . For the rest of the proof, we will separately consider the case $A = 0$ and the case $A > 0$ , $m > 2$ .
Case $A = 0$ : In this case, by (4.13) we get
provided
where
Note that here we only consider the case $\lambda < 0$ and $C_0 > 0$ . This is because for $\lambda \geqslant 0$ or $C_0 \leqslant 0$ , equation (4.14) does not have any real-valued non-trivial periodic solution $v \geqslant 0$ . In particular, since $0 \leqslant f(z) = (1+z^2)^{-\frac{1}{2}} \leqslant 1$ , the case $\lambda \geqslant 0$ and $C_0 \leqslant 0$ cannot lead to any real-valued solution $v \geqslant 0$ to (4.14). For the case $\lambda > 0$ and $C_0 > 0$ , we have $v \leqslant v_1$ or $v \geqslant v_2$ where $v_1$ and $v_2$ are defined in (4.15) and $v_1 < 0 < v_2$ . Smooth periodic solutions that satisfy these criteria do not exist. Moreover, for $C_0 < 0$ and $\lambda < 0$ , we have $v_1 < 0 < v \leqslant v_2$ and non-negative smooth periodic solutions do not exist for this case.
From here we arrive at
This equation has a periodic solution with the period
In the limit case when the droplet touches the substrate surface, we have the minimum and the maximum of the periodic solution given by
and the corresponding $\lambda^*$ satisfies
A typical smooth periodic solution $v(\xi)$ defined on $a \leqslant \xi \leqslant a + \tau$ with $v(a)= v(a +\tau) = v_1 = r_0$ , $v'(a)= v'(a +\tau) = 0$ , and $P(v) =(a,a+\tau) \subset (0, L) $ is given by
If we set $C_0 = \tfrac{r_0}{2}$ , then from (4.18) and (4.19) we obtain $\lambda^* = -r_0^{-1}$ and the trivial solution $v \equiv r_0$ .
Case $A > 0,\, m > 2$ : Similar to the case $A = 0$ , we only need to consider the case $\lambda < 0$ and $C_0 > 0$ . In this case, by (4.13) we get
provided
So, if
where
and $ - \tfrac{(m-1)^2}{2m(m-2)} < \lambda C_0 < 0$ , then there exist $v_1^*(A) $ and $ v_2^*(A) $ such that
and (4.20) is true for $v_1^*(A) \leqslant v \leqslant v_2^* (A) $ . From here we arrive at
This equation has a smooth periodic solution with the period
Taking the minimum of the solution $v^*_1(A) = r_0$ , we have
If we set $v(a)= v(a +\tau) = v_1 = r_0$ and $P(v) =(a,a+\tau) \subset (0, L) $ , then our smooth periodic solution $v(\xi)$ is given by
and $v'(a) = v'(a+\tau) = 0$ . Furthermore, if $C_0 = \frac{1}{2}(r_0 + \frac{A\,m}{m-2}r_0^{2-m})$ , then we obtain the trivial solution $v \equiv r_0$ on (0,L).
On the other hand, if $\lambda > 0$ and $C_0 > 0$ , then
provided
Since $v_2 < 0 < v_1 $ , smooth periodic solutions do not exist.
If $\lambda =0$ and $C_0 > 0$ then
provided
Therefore, we do not have any smooth periodic solutions.
4.3 Convergence to an energy minimiser without gravity
Next, we study the long-time behaviour of solutions to (4.1). We start by considering the case without the gravitational term ( $\mu = 0$ ).
Theorem 2 Let $m > 2$ , $\mu =0$ , $v(\xi,t)$ be a weak solution to (4.1) with periodic boundary conditions and $v_{\min}(\xi)$ be a solution from Lemma 4.2. Assume that
Then there exist $B_0 > 0$ and $B_1 > 0$ such that
and
Proof of Theorem 2. Let us denote
Similar to (4.6), we deduce that
where
Then by (4.25), similar to (4.7), we arrive at
Compare (4.26) with a solution of the following problem
where $ y(t) := \mathcal{E}(v | v_{\min}) \geqslant 0$ ,
As $\beta = 0$ provided $ \mathcal{E}(v_{\min}) \geqslant \mathcal{E}^*$ then from (4.27) we have
Solving (4.28), we deduce that
As a result, we arrive at
where
provided
Since the energy $\mathcal{E}(v)$ is bounded by (4.29) and $v\geqslant r_0 > 0$ , it follows that $\| v \|_{W_1^1(0,L)} \leqslant \bar{C}$ , where $\bar{C}$ is a constant, and the conclusion (4.24) holds.
4.4 Convergence to a travelling wave with gravity
The case with the gravitational term $(\mu = 1)$ is more complicated and will require us to introduce additional conditions to compare to the case without gravity. Let us denote by
where $\tilde{\mathcal{E}}(v)$ is a modified energy functional that incorporates surface energy, the stabilisation mechanism and a gravitational potential energy represented by $\int_0^L \tfrac{1}{2}v^2 F(\xi,t)\;d\xi$ .
Lemma 4.3 Let $m > 2$ and $\mu =1$ . Assume that
and
Then, for every $M$ , the functional $\tilde{\mathcal{E}}$ attains its minimum $v_{\min}$ on $\mathcal{M}$ .
Proof of Lemma 4.3. Multiplying (4.1) with $\mu =1$ by $- (J - F(\xi,t))$ and integrating on $\xi$ , we deduce that
Then from (4.34), due to (4.33), we get
Moreover, taking into account
we deduce that
As the functional $\tilde{\mathcal{E}}(v)$ is non-increasing and bounded from bottom it attains its minimum $v_{\min}$ on the set $\mathcal{M}$ .
The Euler–Lagrange equation for $\tilde{\mathcal{E}}(v)$ under the constraint $ \int \limits_0^L { v^2 d\xi} = M$ is given by
where $\lambda$ is a Lagrange multiplier and
We also need to incorporate the constraint $v \geqslant r_0 > 0$ . If $v \in \mathcal{M}$ , we decompose (0,L) according to the value of v into the set P(v) and the set Z(v). As equation (4.36) is not autonomous, we have $meas (Z(v)) =0$ .
Lemma 4.4 Let $m > 2$ , $\mu =1$ , and $0 \leqslant A < r_0^{m-1}$ . The functional $\tilde{\mathcal{E}}$ has the minimiser $v_{\min}$ on $\mathcal{M}$ such that $v_{\min} \in C^1[0,L]$ , $v_{\min}(0) = v_{\min}(L) $ and $v'_{\min}(0) = v'_{\min}(L)=0$ and solves (4.36) with
where
provided
Proof of Lemma 4.4. Note that equation (4.36) has a periodic solution provided
Multiplying (4.36) by $-v\,v_{\xi} $ $(v_{\xi} \neq 0)$ , we obtain that
We will look for the first integral to (4.38) in the form:
where $a(0)=a(L)$ and $b(0) = b(L)$ . Substituting (4.39) into (4.38), we find that
By $b(0) = b(L) $ we have
In particular, from (4.37) and (4.40) it follows that
As a result, by (4.39) we arrive at
If $v(0) = v(L) = r_0$ , then $v'(0) = v'(L) = 0$ provided
5 Numerical studies
To simulate the fibre coating dynamics and explore beyond the analytical results presented in previous sections, we numerically investigate the problem (2.1)–(2.3). Specifically, we are interested in the structure of energy minimisers and the travelling wave solutions governed by the PDE (4.1).
Firstly, we investigate the energy minimisers $v_{\min}(\xi)$ discussed in Lemma 4.2 and their structures. For $\lambda \equiv \lambda^*$ and $A = 0$ , the profile of $v_{\min}(\xi)$ has a unique maximum $\max_{\xi}(v_{\min}) = v_1$ and minimum $\min_{\xi}(v_{\min}) = v_2$ , where $\lambda^*$ , $v_1$ and $v_2$ are defined in Lemma 4.2 by (4.18) and (4.17), respectively. For given values of $\lambda^*$ and $r_0$ , the corresponding value of $C_0 = r_0 + \lambda^* r_0^2 / 2$ is obtained based on (4.18), and the period $\tau$ of the minimiser is obtained from (4.16). Then we numerically solve the ODE (4.10) for the solution profile $v_{\min}(\xi)$ on $ -\tau/2 \leqslant \xi \leqslant \tau/2$ subject to the Neumann boundary conditions $v'(\xi) = 0$ at $\xi = \pm \tau/2$ . When other system parameters are fixed, a larger absolute value of $\lambda^*$ , $|\lambda^*| = -\lambda^*$ , leads to an energy minimiser $v_{\min}(\xi)$ with a smaller magnitude and a smaller period (see Figure 1). The comparison between the two plots with the fibre radius $r_0 = 0.2$ and $r_0 = 1$ in Figure 1 also shows that with identical $\lambda^*$ values, a larger value of $r_0$ leads to a $v_{\min}(\xi)$ profile with a larger period and smaller spatial variations. This implies that with identical dynamic pressure, droplets on a thick fibre tend to have a smaller bead size than those on a thinner fibre.
Figure 2 (left) presents the effects of the stabilisation term $A/u^m$ , showing that with other system parameters fixed, a positive $A>0$ yields smaller spatial variations in $v_{\min}$ . This is also consistent with the observation drawn in [Reference Ji, Falcon, Sadeghpour, Zeng, Ju and Bertozzi10]. Similar to the $A = 0$ case, the energy minimisers $v_{\min}(\xi)$ for $A > 0$ presented in Figure 2 (left) are numerically obtained by solving the ODE (4.10) for $-\tau/2 \leqslant \xi \leqslant \tau/2$ , where the period $\tau$ is given by (4.21) and the constant $C_0$ is obtained from (4.11) with the fibre radius $r_0 = 1$ , $\lambda=\lambda^* = -0.3$ , and $m = 3$ . With different values of A and $r_0$ , Figure 2 (right) shows a comparison of the period $\tau$ of energy minimisers based on (4.21) for a varying $\lambda = \lambda^*$ that satisfies (4.11). This plot shows that when other system parameters are fixed, larger values of $r_0$ and A typically lead to a larger period $\tau$ for the energy minimiser.
Next we explore the dynamics of the problem (2.1)–(2.3) without gravity effects, that is, without the last term $[Q(u)]_x$ . Typical dynamic simulation results of (2.1) without gravity modulation subject to periodic boundary conditions are plotted in Figure 3. The initial data for this simulation are given by $u_0(x) = v(x)$ which satisfies the ODE (4.10) with $C_0 = 0.5$ , $\lambda = -0.5$ over $0 \leqslant x \leqslant L$ for the period $L = \tau = 10.79$ obtained from (4.16). Other system parameters are $A = 0$ , $r_0 = 0.2$ , $\sigma = 0.01$ and $M = \int_0^L u^2\;dx = 68.6$ . We note that the parameter pair $(\lambda, C_0)$ do not satisfy the condition (4.11), which implies that the initial profile $u_0(x)$ is not an energy minimiser of the system. Starting from the initial profile, the dynamic solution quickly evolves into a two-hump solution with an additional local maximum at the edge $x = 0$ and $x=L$ . Moreover, a local singularity, $u = r_0$ , is approached at two global minima near $x = 0.9$ and $x = 9.9$ where the two hump solutions are connected. An inset plot is also included in Figure 3 (left) to display the dynamics near the region where the two humps are connected. The history of the energy $\mathcal{E}(t)$ in Figure 3 (right) shows that the energy is dissipating as the solution approaches to the two-hump profile. This confirms the conclusion in Theorem 2. These numerical studies are conducted using fully implicit finite differences with a uniformly spaced grid, and the numerical scheme is second-order accurate in space and first-order accurate in time. In our numerical studies, we also observed cases when a single-hump initial profile can evolve into three or more humps. We suspect that linear stability analysis for steady states of the governing PDE can lead to an estimate of the number of humps and the structure of long-time dynamics. However, a systematic investigation of this problem is beyond the scope of this paper, and we may further study this aspect of the problem in the future.
Figure 4 (right) presents the plot of $\lambda(\xi)$ corresponding to the dynamic solution in Figure 3 at time $t = 20$ . It shows that $\lambda(\xi)$ approaches a piece-wise constant function where the values of critical $\lambda^*$ are given by $\lambda^*_A = -0.523$ and $\lambda^*_B = -3.955$ , respectively. Using these critical $\lambda^*$ values, we numerically solve for the energy minimisers $v_{\min}(\xi)$ derived in Lemma 4.2 that correspond to $\lambda = \lambda^*_{A,B}$ and compare them (see Figure 4 (left)) against the long-time dynamic solution presented in Figure 3 (left). The two solutions, $v_{\min}(\lambda=\lambda^*_A)$ and $v_{\min}(\lambda=\lambda^*_B)$ , are connected at the point where the common minimum value $v = r_0$ is reached (see the inset plot in Figure 4 (left)). This comparison shows that the two-hump solution obtained from the direct numerical calculation is in good agreement with the energy minimiser characterised by $v_{\min}$ with a piece-wise constant function $\lambda(\xi) = \lambda^*_{A,B}$ .
The gravity effects incorporated in (2.1)–(2.3) can lead to a travelling wave solution at a speed V characterised in the travelling wave PDE (4.1). To better understand the influence of gravity effects in the model, we numerically solve the equation (2.1) starting from the two-hump profile obtained from the simulation shown in Figure 3 where gravity effects are excluded. The system parameters are also set to be identical to those used in Figure 3. The simulation shows that the dynamic solution quickly converges to a travelling wave at a constant speed $V = 5.12$ . A comparison of this travelling wave solution and the two-hump profile without gravity is shown in Figure 5 (left), where the spatial symmetry around $x = L/2$ in the travelling wave solution is broken due to the presence of gravitational effects. Figure 5 (right) depicts the corresponding $\lambda$ profiles for the two solutions. This plot reveals that without gravity the $\lambda$ profile for the two-hump solution approaches a piece-wise constant function, while the travelling wave solution only keeps the large hump corresponding to $\lambda=\lambda_A^*$ and the small hump with $\lambda=\lambda_B^*$ is saturated by the gravity. In addition, the global minimum of the travelling wave solution is above the fibre radius parameter $r_0 = 0.2$ , which indicates that the gravitational effects prevent the singularity observed in the no-gravity case from happening.
Finally, we investigate the convergence to travelling wave solutions with the presence of gravity using the PDE (4.1) in the moving reference frame with $\mu = 1$ . A travelling wave solution $v_{\min}(\xi)$ at a propagating speed V is a steady state of PDE (4.1) and satisfies the fourth-order ODE (5.1) for $v(\xi)$ on a periodic domain $0 \leqslant \xi \leqslant L$ ,
The profile of the travelling wave is determined by the period L and the mass $M = \int_0^L v^2\;d\xi$ . For a given (L, M) pair, we numerically solve the travelling wave ODE (5.1) for $v_{\min}(\xi)$ and the speed V simultaneously on a periodic domain $0\leqslant \xi \leqslant L$ using finite differences. A continuous family of travelling wave solutions to (5.1) parametrised by (L, M) can be identified. Using a branch continuation method, we numerically track the travelling waves as the parameter L or M changes. Similar parametric and pseudo-arclength continuation methods are commonly used for bifurcation analysis [Reference Ji and Li11,Reference Ji and Witelski14,Reference Keller16]. For $L = 9$ and $M = 57.2$ , we obtain a travelling wave solution $v_{\min}(\xi)$ with the associated velocity $V = 4.24$ (see the solid curve in Figure 6 (top left)). This velocity agrees with the analytical expression for V provided in Lemma 4.4. The other system parameters are identical to those used in Figure 3.
Starting from a perturbed initial condition
where $\bar{\epsilon}$ is the magnitude and $\bar{f}/L$ is the frequency of the perturbation, we numerically solve the travelling wave PDE (4.1) and observe that the PDE solution converges to the travelling wave $v_{\min}(\xi)$ as $t\to \infty$ . Here we select the magnitude of the perturbation
such that the initial data satisfy the condition $\int_0^L v_0^{2}\;d\xi = \int_0^L v^2_{\min}\;d\xi = M$ which is necessary for the convergence to the travelling wave $v_{\min}$ due to the conservation of mass condition. For the simulation shown in Figure 6, we set $\bar{f} = 1$ which corresponds to the magnitude $\bar{\epsilon} = 0.51$ based on the formula (5.3) for $L = 9$ . To calculate the modified energy $\tilde{\mathcal{E}}(t)$ defined in (4.31), we specify $\nu(t)$ in the definition of $F(\xi,t)$ in (4.30) based on the formula
so that condition (4.32) is satisfied. Figure 6 (top right) depicts the evolution of $F(\xi,t)$ in time as the PDE solution converges to the travelling wave, showing that the condition $F(0,t) = F(L,t) = 0$ in (4.34) is satisfied. The plot in Figure 6 (bottom) shows that as the travelling wave $v_{\min}(\xi)$ is approached, both the modified energy $\tilde{\mathcal{E}}(t)$ and the corresponding $\int_0^L v^2 F_t\;d\xi$ decay in time, which is consistent with the result shown in Lemma 4.3. We note that while the simulations in Figures 4 and 6 are conducted with $A = 0$ , similar dynamics are observed for simulations with small stabilisation parameter $A > 0$ .
6 Conclusions
The main contribution of this paper is showing the existence and long-time behaviour of non-negative weak solutions for the generalised nonlinear PDE (2.1)–(2.3) using a priori estimates for energy-entropy functionals. Typical numerical studies of the energy functional minimisers and dynamic simulations of the PDE with and without gravitational effects are presented in support of the analytical results. The travelling wave solutions of the model are investigated both analytically and numerically. As $t\to \infty$ , with proper system parameters and initial conditions, the solution to (2.1) converges to a travelling wave solution characterised by (4.1).
Acknowledgements
H. Ji was supported by the Simons Foundation Math+X investigator award number 510776. The authors are also grateful to Prof. Andrea L. Bertozzi for helpful discussions.
Conflict of interest
None.