ON THE BENILOV-VYNNYCKY BLOW-UP PROBLEM

We study an initial-boundary value problem for a fourth-order parabolic partial differential equation with an unknown velocity. The equation originated from the linearization of a two-dimensional Couette flow model, that was recently proposed by Benilov and Vynnycky. In the case of a 180◦– contact angle between liquid and a moving plate Benilov and Vynnycky conjectured that the speed of the contact line blows up to infinity in finite time. In this paper we present numerical simulations and qualitative analysis of the model. We show that depending on the initial data and parameter values different long time behaviors of velocity can be observed. The speed of the contact line may blow up to infinity or converge to a constant.

1. Introduction. The contact line is the triple junction between solid, air and liquid flow. It is well known that Navier-Stokes equations with classical boundary conditions are not applicable if the free boundary of the flow intersects with a rigid boundary, resulting in the contact line. In 1980 Benney and Timson [3] analyzed the viscous liquid flow near the contact line and showed that, if the contact angle is 180 • (the angle between the solid and the liquid interfaces), the contact line singularity, that is well known for the zero contact angle, does not arise and as a result the interface propagation is well-defined. Their local analysis did not include an asymptotic behavior of the contact line velocity. In a recently published paper, Benilov and Vynnycky [4], under an assumption of lubrication approximation regime, complemented the Benney and Timson results via asymptotic analysis of the global flow. Among other characteristics of the global flow they also determined the contact line velocity.
The Couette flow is represented schematically in Figure 1, where two parallel horizontal rigid plates are separated by the distance H. The upper plate is moving to the left with a velocity U 1 and the lower plate is moving to the right with a velocity U 2 . The volume between these plates is filled with an incompressible viscous fluid on the left and with vacuum on the right, with the free boundary separation. The contact line is located on the upper plane and the position of the contact line in the reference frame that is co-moving with a contact line is fixed at the point x = 0. Under an assumption that a velocity of the upper plane matches a velocity of the contact line we have: U 1 = −V (t) and U 2 = U −V (t), where U is a constant velocity of the lower plate relative to the upper plate. The time evolution of the profile of the liquid/vacuum free boundary is described by the graph of the function h(x, t) for x > 0, and where h is the thickness of the liquid film.
We study the following initial-boundary value problem that was derived in [4, Eqs.
with two unknowns h(X, T ) and v(T ). With the scaling In the special case when the contact line velocity was assumed to be a constant V (t) = V 0 the problem above was studied in [8]. The class of self-similar solutions for this partial differential equation (self-similar solutions do not satisfy the boundary conditions of the original problem) was constructed in [7].
Benilov, in personal communications, pointed out that a classical solution of the problem above (if it exists) should have an infinite number of constraints. It follows from the boundary conditions and from the partial differential equation that if a solution h(x, t) is a classical one it should satisfy the condition h xxxx (0, t) = 0. By differentiating this new condition with respect to t and by using the partial differential equation above as a substitution for h t one obtains an infinite series of constraints (under an assumption of an infinite smoothness of the solution).
The existence of an infinite series of constraints for a solution of a partial differential equation is unusual, but this is not a unique case. A similar property has been noted for the so-called Ostrovsky equation for waves in a rotating ocean or in a channel with bottom topography [6,2], as well as for Kadomtsev-Petviashvili equation [5,1]. Finally, there are numerous models with a finite number of constraints, associated mostly with a scale of oceanic dynamics (see [10] and references therein). In all such cases, the constraints reflected an adjustment of the solution by fast dynamics, which were present in the original (exact) problem, but have been scaled out while deriving a slow-time asymptotic model. If the initial data taken does not comply with all the boundary constraints, it instantaneously evolves into a state satisfying all of them. However, since numerical methods cannot, generally, handle infinitely fast evolution, the adjusted state is not computed accurately. Once the adjustment is complete, the numerical solution begins to satisfy the equation accurately enough, but this adjustment makes the adjusted state different from the one originating from the initial condition given.
The structure of the article is as follows. In section two we prove non-existence of physically relevant stationary solutions in the Benilov-Vynnycky model and present analytical stationary solutions in non-physical regime. In the third section using energy method and functional inequalities we analyze a finite-interval approximation of the original problem. We derive estimates for V (t) and for the existence times of solutions for different types of initial data. In section four we propose a simple numerical method, based on the finite difference approach, to solve the initialboundary value problem of the Benilov-Vynnycky model. We show that depending on the initial data and parameter values, the magnitude of speed of the contact line may blow up to infinity or converge to a constant. The short discussion is presented in section five with a comparison of blow-up rates of numerical V (t) with the logarithmic rate V (t) c 1 ln(t * − t) + c 2 predicted in [4] and with a power law rate V (t) c 1 (t * − t) −1/2 + c 2 predicted recently in [9].
2. Formulation of the problem and analysis of stationary solutions. We consider the following initial-boundary value problem on R + : with two unknowns h(x, t) and V (t).
The solution h(x, t) has a physical meaning only if 0 ≤ h(x, t) ≤ 1 for t ≥ 0 that implies lim The stationary solution of the problem above must satisfy The general solution can be written as From the condition u(0) = 1 it follows that The condition u x (0) = 0 provides that If V > 0, the conditions lim Notice that two boundary conditions at x → ∞ lead to the same constraint and do not identify V uniquely. Furthermore, this solution has lim It is worth mentioning that if V < 0, the conditions lim This is in contradiction with the conditions (4,5). Thus there is no physical stationary solution for this problem.
In the following section, we will study the problem above restricted to a finite interval. We will explore numerically how properties of solutions depend on the length of the truncated interval. 3. Formulation of the finite-interval approximation and some qualitative analysis. We consider the following initial-boundary value problem on the finite interval [0, L]: with two unknowns h(x, t) and V (t). Let us show that V (t) can be positive only during a finite time interval provided that 0 ≤ h ≤ 1.
Due to non-negativity of V (t) and 1 − h 2 (L, t) this implies that and integrating over [0, T ] we obtain the claimed upper bound for the time T .
Proof. Integration of the equation over domain (0, L), yields then the total mass is increasing.
of the problem (P fin ) then the following inequality holds true Proof. Multiplying Equation (6a) by h xx and integrating over finite domain [0, L], As x = 0 is the attachment point of the liquid film to the upper moving plate that leads to the assumption that h xx (0, t) ≤ 0 and in this case we will obtain the energy dissipation If we apply the Poincare' inequality to the h It follows from (10, 11) that This implies that L we derive the upper and the low bounds for the thickness of the liquid film If the solution was global in time the bounds (9) in the theorem above would imply the uniform in time convergence toward h = 1. We can also use (9) to obtain the estimation of total mass: The approach we used above is not straightforwardly applicable to the semiinfinite domain (P) due to the integrability problem. If we assume that a short-time initial dynamic is local in space and lim x→∞ h(x, t) = h ∞ (does not depend on time), then we can show that the contact line velocity V (t) can be positive only during a finite time interval. Indeed, for any constant value of h ∞ we can introduce a new variableȟ = h − h ∞ . Due to the invariance of the partial differential equation, with respect to the shift and scaling transformations, the boundary conditions can be preserved. We omitȟ-notations below.
and integrating over [0, T ] we obtain the claimed upper bound for the time T .
Note that because a blow-up of V (t) is possible only if there is a blow-up of the rate of change of the total mass or if h ∞ = 1.
Theorem 3.6. There are countably many stationary solutions with V > 0 for the finite interval problem (P fin ).
Proof. Denote a = V 1 3 , the solution (2.8) can be written as The boundary conditions u x (0) = 0, u x (L) = 0, and u xxx (L) = 0 lead to and θ = aL. The determinant of the matrix has to be zero in order to have nontrivial solutions, i.e., Notice that θ = 0 is a solution. Furthermore, since f 1 (θ) = sin √ Theorem 3.7. If u is the stationary solution of (P fin ), h is the time-dependent solution of (P fin ), and lim Proof. In the finite interval [0, L],    Multiplying both sides by w xx and then integrating from 0 to L, we obtain By using Cauchy-Schwarz for the right hand side term and Poincare inequality for

MARINA CHUGUNOVA, CHIU-YEN KAO AND SARUN SEEPUN
Since Thus We then have Taking ε < π L , we have lim 4. Numerical method. In this section, we apply a semi-implicit method to study the solutions of the problem (P fin ). The computational domain is (0, L). A uniform grid of points x j = j∆x where 0 ≤ j ≤ N and N = L ∆x is used. The step size in time is ∆t. The discretization of Eq. (6a) in time yields where h n j is the numerical approximation of h(x j , t n ). The implicit discretization of h xxxx is chosen to ensure that the step size can be chosen reasonably, i.e., ∆t ∼ O(∆x) instead of ∆t ∼ O(∆x) 4 . The explicit discretization of h x is chosen to avoid nonlinearity in unknown variables h n+1 j and V n+1 . Thus the updating formula is where I is the identity matrix and D 4 , D 1 are numerical operators which approximate fourth-order and first-order differential operators respectively. We simply use the five-point central scheme for D 4 and the first order upwind scheme for D 1 which uses different discretization depends on the wind direction −V , i.e., and For the boundary conditions, we use the following discretization.
The time step is chosen to satisfy the CFL condition ∆t ≤ ∆x/|V | and the number of interval N used is 128.
Notice that the initial velocity V 0 is not required when we solve (15) with the aforementioned boundary conditions. The discretized equation (15) and boundary conditions (16) form a linear system with N + 2 unknown h n+1 j (j = 0, ..., N), V n+1 , and N + 2 equations, which can be solved easily. The first step of calculation sometimes yields a large change in thickness h when the given initial data does not comply with all the boundary constraints. This is because the algorithm seeks the solution h and V to satisfy the equation and boundary conditions in one step. The contradiction between initial values and a partial differential equation was previously analyzed in [11,12] and it is also imposed an additional difficulty on construction of a numerical solution of two-phase Stefan problem.
We first demonstrate the first order convergence of our numerical scheme with two initial conditions. The first initial condition is a 5-th order polynomial function which satisfies all the boundary conditions and h(L) = 0.2, i.e., where with L = 1 while the second initial condition is 0.8 cos 10 ( π 2 x) + 0.2 on the domain [0, 1]. We choose the second initial condition because it generates a wave front which looks like the one shown in [4]. Since the exact solution is not available, numerical tests were conducted on four different sizes of meshes: ∆x, 2∆x, 4∆x and 8∆x with N = 256. The step size ∆t = c∆x is chosen to satisfy the CFL constrain, i.e. c < 1/|V |. We choose ∆t = 5 × 10 −6 and ∆t = 1 × 10 −5 which are small enough for these two initial conditions, respectively. Denote the solution with mesh size ∆x as h ∆x . We compute the differences between solutions in L 2 -norm for various times and list them in Table 1 and Table 2. For the polynomial initial condition, the blow-up time for the velocity V happened before t = 0.02 thus the value is not listed (In Table 1, we denote it as N/A). We use the notation order 1 (order 2 ) for the base-2 logarithm of ratio of consecutive differences between three solutions obtained by grid sizes ∆x, 2∆x, and 4∆x (by grid sizes 2∆x, 4∆x and 8∆x ). It is clear that our scheme achieves the first order accuracy for both initial conditions at any given time which is away from the blow up time.   Table 1. Accuracy test for the algorithm with the 5-th order polynomial initial condition which satisfies five boundary conditions and h(1) = 0.2.  Table 2. Accuracy test for the algorithm with the initial condition 0.8 cos 10 ( π 2 x) + 0.2.
In Figure 6, we demonstrate how the evolution of h(x, t) and V (t) varies with respect to L with the initial condition as the 5-th order polynomial which satisfies the boundary conditions and h(L) = 0.2 for L = 2, and 6 . We observe the finitetime blow up for V (t) in both cases. For L = 2, V (t) approaches −∞ while, for L = 6, V (t) approaches ∞. Note that the initial condition for the case L = 6 is nonphysical. In Figure 7, we change the initial condition to 0.8 cos 10 ( π  Figure 10, we provide a comparison for solutions in the short (L = 4), medium (L = 16) and long (L = 64) intervals. As we discussed before, the first step of calculation sometimes yields a large change in the thickness h when the given initial data does not comply with all the boundary constraints. We thus compute the solution at t = 10 −3 with the initial condition 0.8 cos 10 ( π 2 x) + 0.2 in the interval [0, 1] and 0.2 for the rest of the interval [1,128] with the mesh size ∆x = 1/32 (4097 grid points in total) and the time step ∆t = 10 −5 first. Then we use part of this solution (the first 129, 513, and 2049 grid points) as the initial conditions for the intervals L = 4, 16, and 64. The evolution of h(x, t) and V (t) are shown in Figure 10 in the interval [0, 6] since the solutions stay close to a constant after 6.
Notice that in the infinite interval, the boundary condition lim x→∞ h x = 0 implies lim x→∞ h x x = 0. However, this is not true for a finite interval truncated problem. Nevertheless, h x x(L, t) becomes closer to zero as L becomes large. The second order finite difference estimations give 3.27 × 10 −4 , −1.05 × 10 −11 and 2.18 × 10 −14 for L = 4, 16 and 64 at t = 0.0625, respectively. We can see that the solutions h are very close to each other for different choices of L and the difference of V only becomes noticeable when the time is close to the blow up time. In Figure  11, we use the numerical results for L = 64 and do least-square fitting by using V (t) c 1 ln(t * − t) + c 2 [4] and V (t) c 3 (t * − t) −1/2 + c 4 [9], respectively, for 800 grid points before the blow up time. With the choice t * = 0.064563, the best fittings are V (t) 15.14 ln(t * − t) + 61.54 denoted by blue line and V (t) −0.5459 ln(t * − t) − 14.51 denoted by green line, respectively. We see that the blow up rate of the contact line velocity V (t) is between the logarithmic one proposed in [4] and the power law one in [9]. From the above simulations, we observe that |V | blows up in finite time for physical initial conditions while V may blow up in finite time or stay bounded for nonphysical initial conditions. 6. Discussion. We have thus explored the behavior of viscous liquid film between two parallel horizontal moving plates. In particular, we showed analytically and numerically that the contact line velocity V (t) exhibits a finite time singularity if the contact angle is 180 • . We have also estimated the local existence times of solutions for the related model (1a) and studied stationary solutions.
First, it is suggested by the results shown in Figure 11 that the blow up rate of the contact line velocity V (t) is between the logarithmic one proposed in [4] and the power law one in [9]. Rigorous asymptotical analysis of the blow-up rate of V (t) is not straightforward and is open for future study. Second, in non-physical regime the steady state seems to be an attractor but we were able to prove this only under additional assumptions. Finally, one can develop more accurate numerical computations to better resolve the singular behavior of V (t).  Figure 11. Left: the red line corresponds to the numerical contact line velocity V (t), the blue line is the least-square fitting of the numerical data for V (t) using V (t) c 1 ln(t * − t) + c 2 and the green line is the least-square fitting using V (t) c 3 (t * −t) −1/2 +c 4 . Right: the semilogy plot to see different fittings more clearly.