STABILITY OF SOLUTIONS TO THE RIEMANN PROBLEM FOR A THIN FILM MODEL OF A PERFECTLY SOLUBLE ANTI-SURFACTANT SOLUTION

. In this article, we consider a quasilinear hyperbolic system of partial diﬀerential equations governing the dynamics of a thin ﬁlm of a perfectly soluble anti-surfactant liquid. We construct elementary waves of the corresponding Riemann problem and study their interactions. Further, we provide exact solution of the Riemann problem along with numerical examples. Fi-nally, we show that the solution of the Riemann problem is stable under small perturbation of the initial data.


1.
Introduction. In this article, we are concerned with the solution to a Riemann problem and elementary wave interactions for a quasilinear hyperbolic system of partial differential equations (PDEs) governing the dynamics of a thin film of a perfectly soluble anti-surfactant. Anti-surfactants arise in many practical situations, including many salts, i.e., when common salt is added to water [15,12], when water is added to short-chain alcohol [7]. Conn et al. [3] proposed a mathematical model for an anti-surfactant solution. Subsequently Conn et al. [4] obtained exact solutions to a family of Riemann problems for a reduced version of their model which describes a thin film of a perfectly wetting anti-surfactant solution (i.e., the case in which the surface concentration of anti-surfactant is identically zero). Specifically, Conn et al. [4] studied a governing system of PDEs given in Cartesian coordinates by ∂h ∂t where x and t are independent space and time coordinates and the dependent variables h and c, respectively, denote the dimensionless film thickness and concentration of the solute. Here, P e is the Péclet number that indicates the ratio of advection to diffusion, while the Capillary number Ca indicates the relative effect of viscous drag forces versus surface tension forces. If diffusion effects are neglected we have P e → ∞ and further in case of large Capillary number Ca → ∞, the governing system of PDEs (1) and (2) reduces to ∂h ∂t + ∂ ∂x ∂c ∂t Note that (3) represents the conservation of mass of the fluid whereas (4) represents the conservation of mass of the solute. Let us denote the gradient of concentration of the solute by b(x, t), i.e., ∂c/∂x = b(x, t) and differentiate (4) with respect to x to obtain a system of PDEs in a conservative form [4] ∂h ∂t + ∂ ∂x In case of non-uniform surface tension, the surface tension is a function of temperature or concentration or both. This phenomena is typically known as Marangoni effect [23]. It is known from literature that the surface tension, denoted by Γ, varies with concentration as Γ = k 1 − k 2 c, k 1 , k 2 ≥ 0.
In case of balancing the shear stress with the Marangoni effects, we typically have where ∇ s denotes surface gradient operator and [ς] denotes jump in shear stress. Using (6) in (7), we have [ς] = k 2 b. Therefore, the shear stress at the free surface of the film is positive if the value of b is positive. Also due to the Marangoni effect, the fluid will flow from left to the right for positive values of b and vice versa for negative values of b. One may see a similar type of anti-surfactant properties that arise in certain resins that are included in solvent-based paints [17,29,8].
Elementary wave interactions for Riemann problem is not only key mechanism to determine the qualitative properties of physical variables but it exhibits certain essential features of the solution. Further, understanding relevant properties based on wave interactions play a significant role in several practical applications. In particular, in the field of two-phase flows [16,11], magnetogasdynamics [20,13], nonlinear elasticity [22], traffic flow [1] etc. Hence, it is a subject of great interest both from the mathematical and physical point of view. Various methods to deal with such system of equations very much relate to the mathematical theory of quasilinear hyperbolic system of PDEs. Raja Sekhar et al. [19] discussed the interaction of elementary waves for shallow water equations. A similar analysis for blood flow equations in one dimension is done by Raja Sekhar et al. [18]. Using the method of characteristics, the interaction of elementary waves and stability of Riemann solution for Aw-Rascle model and chromatography equations are investigated in [27,28,25]. The stability of rarefaction waves for Navier-Stokes-Poisson equations may be found in [5]. Recently, the structural stability of Riemann solution for nonlinear elasticity has been discussed in [21]. More details on Riemann problem and interactions of elementary waves for various practical problems may be found in [24,26,10,2,9,14] and the references cited therein.
In a fundamental work, Conn et al. [3] considered linear stability of an infinitely deep layer of initially quiescent fluid and discussed the occurrence of an instability driven by surface-tension gradients, which occurs for anti-surfactant solutions. Further, Conn et al. [4] derived a set of Riemann solutions to the system (5) for different Riemann initial data and discussed their properties.
In this work, our main purpose is to investigate all possible elementary wave interactions and discuss the stability of Riemann solution for the system of PDEs (5). To the best of our knowledge, this kind of interactions and stability of Riemann solutions for this system (5) have not been studied previously. In order to determine all possible cases of wave interactions for the system (5), we consider the following initial data with three piecewise constant states where > 0 is arbitrarily small. Note that, the initial data (8) is nothing but the local perturbation of the Riemann initial data Our objective is to study whether the solution of Riemann problem (5) together with the data (9) is the limit of the solution of (5) and (8) when the parameter → 0. We construct the elementary wave interactions case by case to resolve this problem.
This paper is structured as follows. In section 2, we briefly express the elementary waves and discuss their properties. We study the Riemann problem (5) and (9) along with numerical examples to see various wave patterns with different possible initial data in section 3. In section 4, we mainly focus on the interaction of elementary waves for all possible cases and investigate the global stability of Riemann solution under the small perturbation of initial data. Finally, we list brief conclusions in section 5.
2. Elementary waves. It may be noted that Conn et al. [4] considered the system of PDEs (5) subject to the Riemann data (9) and derived the corresponding Riemann solution. However, we emphasize here that they have followed the method of characteristics to obtain the solution with a restriction on the initial data, i.e., b l < 0 and b r < 0. In order to study the stability of the perturbed Riemann problem (5) and (8), we need to study the properties of the elementary waves. Hence, we investigate the behavior of the elementary waves in phase plane. Accordingly, we consider the quasilinear form of the system (5) where U and A(U ) denote the primitive variable and the Jacobian matrix, respectively, which are given by The eigenvalues of the Jacobian matrix are λ 1 = bh/2 and λ 2 = 3bh/2 and the corresponding eigenvectors are given by r 1 = (h, −b) T and r 2 = (h, b) T , where T denotes the transposition. It can easily be seen that ∇λ 1 ·r 1 = 0 and ∇λ 2 ·r 2 = 3bh, i.e., 1-characteristic field is linearly degenerate and 2-characteristic field is genuinely nonlinear. Therefore, the solution corresponding to 1-family is always a contact discontinuity and the solution corresponding to the 2-family is either a shock wave or a rarefaction wave. Now, we briefly express these three elementary waves and discuss their properties (for details, we refer [24,26,6]).
2.1. Contact discontinuity. Let us denote the left hand and right hand states of an elementary wave by U l and U , respectively. Since the characteristics on either side of the contact discontinuity are parallel, we have which implies that b l h l /2 = τ = bh/2, where τ denotes the speed of the contact discontinuity. Therefore, for the contact discontinuity we have and the corresponding J curve is shown in Figure 1.

Shock wave.
The Rankine-Hugoniot jump conditions across the shock wave are given by where the speed of shock wave is denoted by σ and we use the notation [h] = h l − h which denotes the jump in h across the shock. Since, (14) consist of two equations in three unknowns σ, h and b, we express σ and b in terms of h as Lax entropy conditions across the shock wave are given by which follow that b l h l /2 < σ < 3b l h l /2 and 3bh/2 < σ. One can observe that, these inequalities hold only when b l > 0 as h is always positive. Therefore, without loss of generality, we assume b l > 0 throughout the article. Since, 3bh/2 < σ < 3b l h l /2 holds across the shock wave, using (15) we obtain h 2 < h 2 l , i.e., h < h l and consequently we obtain b < b l as both h and h l are always positive. Therefore, the shock curve S is given by (see, Figure 1) 2.3. Rarefaction wave. The Riemann invariants corresponding to the 1-and 2characteristic fields are Π 1 = bh and Π 2 = b/h, respectively. Since the Riemann invariants are constant across the rarefaction wave, we have b/h = b l /h l . Further, the characteristic speed increases from left to right across the rarefaction wave, i.e., λ 2 (U l ) ≤ λ 2 (U ). Therefore, the rarefaction wave curve (R) can be expressed as (see, Figure 1).
We summarize the properties of elementary waves as follows: Theorem 2.2. Contact discontinuity curve corresponding to the system (5) is monotonically decreasing and convex while the curve of shock and rarefaction waves are monotonically increasing straight lines, respectively.
3. Riemann problem. In this section, we consider the Riemann problem (5) and (9) and present the corresponding solution. We have already seen from the characteristic analysis that the solution of Riemann problem associated with 1-wave is always a contact discontinuity and the solution corresponding to 2-wave is either a shock or a rarefaction. Therefore, these elementary waves separate the (x, t)-plane into at most three constant states (h l , b l ), (h * , b * ) and (h r , b r ). The choice of the initial data generates two types of wave patterns as shown in Figure 2 From (17) and (18), one can observe that the shock curve and the rarefaction curve coincide with each other and in either of the cases the unknown state Combining (19) and (20), we obtain the unknown state (h * , b * ) given by We can determine the solution inside the rarefaction wave fan by using the slope of the characteristic joining (0, 0) to (x, t). The corresponding characteristic speed is given by Therefore, from (20) and (21), we obtain the solution inside the rarefaction fan as One can observe that, if the solution of the Riemann problem consists of contact discontinuity followed by a shock wave then , solution consists of a shock wave. Similarly, the solution of the Riemann problem consists of contact discontinuity followed by a rarefaction wave if and only if b l h l ≤ b r h r . Therefore, depending upon the choice of the initial data (9), there are two possible wave patterns of the solution to the Riemann problem (5) and (9) which are described as follows:   Table 1. Initial data and solution for the Riemann problem contact discontinuity J followed by a shock wave S and the corresponding solution at (x, t) is given by where τ = b l h l /2 and σ = (b * (h 2 * + h * h r + h 2 r ))/2h * are, respectively, the speeds of propagation of contact discontinuity and shock wave. Case B: If b l h l ≤ b r h r , then the solution of Riemann problem consists of contact discontinuity J followed by a rarefaction wave R and the solution at (x, t) is given by In order to analyze different types of wave patterns for the Riemann problem we consider two test cases. The initial data and the analytical solution for the Riemann problem in the unknown region for the two test cases are listed in Table 1. The analytical solution of the Riemann problem consists of a left contact discontinuity and a right rarefaction wave for test 1. The solution profiles at time t = 0.95 are shown in Figure 3. For test 2, the solution consists of a left contact discontinuity and right shock wave, the solution profiles at time t = 0.95 are depicted in Figure  4.  We have observed that depending on the choice of the initial data, there are two kind of solutions to the Riemann problem. Further, we have presented the analytical expression of the solution for the Riemann problem in (22) and (23). In the next section, we use these solutions to study the wave interactions and stability of solution to the Riemann problem.

4.
Interaction of elementary waves. Here, we consider the system of PDEs (5) with the initial data (8) to study the elementary wave interactions. Let us assume that (h , b ) be the solution of perturbed Riemann problem (5) and (8). Now, we wish to investigate the stability of the solution of the Riemann problem, i.e., whether the solution of the Riemann problem (5) and (9) is the limit of (h , b ) as the perturbation parameter → 0.
Since, we have two local Riemann problems for the initial data (8) at (− , 0) and ( , 0), we realize that four different combinations of the elementary waves are possible. These are given by J + S and J + S, J + R and J + R, J + R and J + S, J + S and J + R. We now discuss these interactions.

Case 1: J + S and J + S
In this case, we consider the interactions between a contact discontinuity followed by a shock wave originating from (− , 0) and a contact discontinuity followed by a shock wave originating from ( , 0). This case occurs only when the initial data (8) satisfies the inequality b l h l > b m h m > b r h r . For a small enough time t, the solution of the initial value problem (5) and (8) can be expressed as follows: (h l , b l ) + J 1 + (h 1 , b 1 ) + S 1 + (h m , b m ) + J 2 + (h 2 , b 2 ) + S 2 + (h r , b r ), where "+" denotes "followed by", and (h 1 , b 1 ), (h 2 , b 2 ) are solutions of the unknown states of first and second Riemann problems, respectively. One can easily compute that ))/2h 1 and τ 2 = b m h m /2 are propagation speeds of S 1 and J 2 , respectively. Hence, S 1 and J 2 interact after a finite time, say, t 1 (see, which yields Therefore, at time t = t 1 we have a new Riemann problem with initial data (h l , b l ) = (h 1 , b 1 ) and (h r , b r ) = (h 2 , b 2 ) and the solution of this new Riemann problem consists of a new contact discontinuity J 3 and a new shock S 3 . We see that the propagation speeds of the contact discontinuities J 1 and J 3 are equal (i. e., τ 1 = τ 3 ). Hence, J 1 and J 3 are parallel and never interact. Now, for t > t 1 the propagation speeds of S 2 and S 3 are respectively given by Therefore, S 3 overtakes S 2 after a finite time, say, t 2 and the interaction point (x 2 , t 2 ) satisfies the following equations Solving (26), we obtain Therefore, at t = t 2 again a new Riemann problem forms with the initial data and (h r , b r ) can be connected only by a shock wave (see, Figure 5(b)) and hence, for t > t 2 , a new shock S 4 occurs (see, Figure 5(a)) with a speed of propaga- ))/2h 3 < 0 implies that contact discontinuity J 3 can not overtake S 4 . Therefore, for t > t 2 , the solution can be expressed as (see, Figure 5 One can observe that, when the perturbation parameter → 0 then both (x 1 , t 1 ) and (x 2 , t 2 ) tend to (0, 0) and the curves J 1 and J 3 coincide. Hence, we can conclude that for a sufficiently large time t, the results of interaction can be expressed as (h l , b l ) + J + (h 3 , b 3 ) + S 4 + (h r , b r ). This result indicates that the limit of the solution of initial value problem (5) and (8) is exactly the corresponding Riemann solution of (5) and (9). Therefore, the solution of Riemann problem (5) and (9) is globally stable under small perturbation.

Case 2: J+R and J+R
Here, we consider the initial data (8) such that the solution of both the local Riemann problems consist of contact discontinuity and rarefaction wave only. Therefore, for sufficiently small time t, the solution of (5) and (8) can be expressed . This situation occurs only when the initial data (8) The propagation speeds of head of rarefaction wave R 1 and contact discontinuity J 2 are respectively given by ω 1 = 3b m h m /2 and τ 2 = b m h m /2. It is clear that ω 1 > τ 2 and they interact at some point, say, at (x 1 , t 1 ) (see, Figure 6(a)) which is determined by A direct computation yields For t > t 1 , the contact discontinuity J 2 is at the onset of interaction with the  rarefaction wave R 1 . Accordingly, during the process of penetration, the curve of contact discontinuity J 2 is determined by From (30a), we obtain Differentiating (30b) with respect to t, we get Using (30a) and (31) in (32), we obtain which implies that the contact discontinuity J 2 will decelerate during the process of penetration. Further, combining (30a) and (30b), we obtain The solution of the initial value problem (34) is given by It may be noted that the curve in equation (35) denotes contact discontinuity during the process of penetration. Also, the equation of the tail of the rarefaction wave R 1 curve is given by The intersection of these two curves will give the end point of penetration. Accordingly, we solve (35) and (36) to get the end point of penetration , say, (x 2 , t 2 ) given by Therefore, the contact discontinuity J 2 sweeps the entire rarefaction fan completely after a finite time t 2 . As b 1 h 1 ≥ b 2 h 2 , for t > t 2 , a new contact discontinuity J 3 and a new rarefaction wave R 3 can occur (see, Figure 6) as a solution of new Riemann problem. Since 3b 3 h 3 /2 = 3b 1 h 1 /2 and 3b m h m /2 = 3b 2 h 2 /2, which implies that the head and tail of the rarefaction wave R 3 coincides with the head and tail of the rarefaction wave R 1 . Further, with τ 1 = τ 3 , the contact discontinuities J 1 and J 3 are parallel.
In the limit → 0, the points (x 1 , t 1 ) and (x 2 , t 2 ) approach (0, 0) and J 1 and J 3 coincide each other. Furthermore, in this limit the speed of propagation of the head of rarefaction wave R 3 and the propagation speed of the tail of rarefaction wave of R 2 are equal and hence R 3 and R 2 coincide each other in this limit. Therefore, for sufficiently large t, the result of interaction can be expressed as (h l , b l ) + J + (h 3 , b 3 ) + R + (h r , b r ). Hence, the solution is stable under small perturbation.

Case 3: J+R and J+S
This case occurs when the initial data (8) satisfies the inequalities b l h l ≤ b m h m and b r h r < b m h m . For small enough time t, the solution of the initial value problem (5) and (8) can be represented as follows: (h l , b l ) + J 1 + (h 1 , b 1 ) + R 1 + (h m , b m ) + J 2 + (h 2 , b 2 ) + S 1 + (h r , b r ). As discussed in Case 2, the interaction of R 1 and J 2 produce a contact discontinuity J 3 and a rarefaction wave R 2 . Moreover, the propagation speeds of J 1 and J 3 are equal as τ 1 = τ 3 and hence J 3 is parallel to J 1 .
(a) Interaction of elementary waves in (x, t)-plane when solution of first local Riemann problem consists of a contact discontinuity followed by a rarefaction wave and the second local Riemann problem consists of a contact discontinuity followed by a shock wave Also the interaction point (x 1 , t 1 ) of R 1 and J 2 is given in (29) and the end point of penetration (x 2 , t 2 ) of R 1 and J 2 satisfies (37). The propagation speed of the head of rarefaction wave fan R 2 and S 1 are respectively given by ω 2 = 3b 2 h 2 /2 and σ 2 = (b 2 (h 2 2 + h r h 2 + h 2 r ))/2h 2 , where h r < h 2 . Straight forward calculations yield (a) Interaction of elementary waves in (x, t)-plane when solution of first local Riemann problem consists of a contact discontinuity followed by a rarefaction wave and the second local Riemann problem consists of a contact discontinuity followed by a shock wave Therefore, R 2 and S 1 interact at a point (x 3 , t 3 ) which can be determined by From (39), we obtain Therefore, for t > t 3 , S 1 begins to penetrate R 2 and during the process of penetration, S 1 propagates with a varying speed which can be determined by From (40a) and (40b), we obtain Usage of (40c) in (41) yeilds Using the properties of the rarefaction wave (please refer Theorem 2.1), it may be noted that dh/dt < 0. Now differentiating (40a) and combining with (40c), we get The estimate (43) indicates that S 1 decelerates during the process of penetration. Integrating (42), we obtain Hence, for t > t 3 , S 1 sweeps the entire rarefaction fan when h r < h 3 (see, Figure 7) after a finite time t 4 given by This indicates that S 1 can never overtake the entire R 2 when h 3 ≤ h r (see, Figure  8) since t → ∞ as h → h r and the propagation speed of the shock is 3b r h r /2 as t → ∞. Therefore, in the limit → 0, for a sufficiently large time t, when h r b r < h l b l , the solution can be expressed as follows (h l , b l ) + J + (h 3 , b 3 ) + S + (h r , b r ) (see, Figure 7), whereas when h r b r ≥ h l b l , the solution can be expressed as (h l , b l ) + J + (h 3 , b 3 ) + R + (h r , b r ) (see, Figure 8). Further, for a sufficiently large time, the shock wave has x + = 3b r h r t/2 as its asymptote in the latter case (i.e., when h r b r ≥ h l b l ). Hence, the solution of the Riemann problem (5) and (9) is stable under small perturbation.

Case 4: J+S and J+R
For a sufficiently small time t, the solution of (5) and (8) can be expressed as (h l , b l ) + J 1 + (h 1 , b 1 ) + S 1 + (h m , b m ) + J 2 + (h 2 , b 2 ) + R 1 + (h r , b r ) under the assumptions b l h l > b m h m and b m h m ≤ b r h r . As discussed in Case 1, S 1 and J 2 interact after a finite time, say, at (x 1 , t 1 ) (see, Figure 10(a)) which satisfies (25). This interaction leads to a new contact discontinuity J 3 and a new shock wave S 2 . The propagation speeds of shock wave S 2 and the tail of rarefaction wave R 1 For t > t 2 , the shock wave S 2 begins to penetrate R 1 with a varying speed of propagation given by Further, from (46a) and (46c) we get Using (46a) and (46b), we get Hence, we see that d 2 x/dt 2 > 0 indicating that the shock wave accelerates during the penetration process. Therefore, as done in Case 3, for a sufficiently large time t, when h r b r ≥ h l b l (i.e., when h r ≥ h 3 ), the solution can be expressed as follows (h l , b l ) + J + (h 3 , b 3 ) + R + (h r , b r ) (see, Figure 9), whereas when h r b r ≥ h l b l (i.e., when h r < h 3 ), the solution can be expressed as (h l , b l ) + J + (h 3 , b 3 ) + S + (h r , b r ) (see, Figure 10) as the limit → 0. Hence, the solution of Riemann problem (5) and (9) is stable under such a small perturbation. We summarize our results as in Theorem 4.1: Theorem 4.1. If the perturbation parameter → 0, then the limiting solution of the perturbed Riemann problem (5) and (8) exactly matches with the solution of the corresponding Riemann problem (5) and (9). Moreover, the initial states (h l , b l ) and (h r , b r ) completely govern the asymptotic behavior of the perturbed Riemann solution which implies that the solution of Riemann problem (5) and (9) is stable under the small local perturbation in the initial data (9) of the Riemann problem.

Conclusions.
We discussed all possible elementary wave interactions using characteristic analysis method and the solution of the perturbed Riemann problem is constructed globally. By letting the limit → 0, it is observed that the solution of perturbed Riemann problem (5) and (8) exactly matches with the corresponding Riemann solution of (5) and (9). Therefore, the limiting solution of perturbed Riemann problem consists of contact discontinuity followed by a rarefaction wave or a shock wave. In the first case, the solutions for h and b are continuous everywhere except for contact discontinuity between regions U l and U * , which propagate towards right into the region x > 0 at a constant speed τ . However, in the latter case the solutions for both h and b are uniform everywhere except across contact discontinuity and shocks, which propagate towards right into the region x > 0 with constant speeds τ and σ, respectively. Moreover, in this case the concentration gradient (b * ) always lies between b l and b r and the value of h * is always greater than both h l and h r , i.e. the film is always thickest in the middle region between the contact discontinuity and the shock. Also, it is noticed that due to positive initial values of b, concentration gradient is always positive which drives the fluid to right. Hence, the asymptotic behavior of solution is completely governed by the initial states (h l , b l ) and (h r , b r ).