Instabilities and pattern formation in viscoelastic fluids

Instabilities and pattern formation in viscous fluids have been a major topic of non-linear fluid dynamics for several decades. The study of pattern formation in viscoelastic thin films offers the opportunity to find new fascinating structures that cannot be observed in viscous fluids. Rayleigh–Taylor and Faraday instabilities, such as the resulting patterns in thin films of viscoelastic fluids, are investigated. We use the long-wave approximation and a Karman–Pohlhausen approach to simplify the mass and momentum equations. The viscoelastic stress tensor is calculated applying the linear Maxwell model. Conditions for the Faraday instability have been found using Floquet’s theorem. It is shown that viscoelastic films can exhibit harmonic resonance under external vibration. Moreover, a simulation of the non-linear problem in 2D and 3D is conducted with a finite difference method. Unstable oscillating Rayleigh–Taylor modes occur in the 2D numerical solution. Furthermore, we find that the wavenumber changes with the relaxation time of the fluid. Faraday patterns in viscous films emerge as regular structures of the surface, like squares or hexagons. Numerical simulations of the viscoelastic fluid also show regular structures. However, they collapse into a chaotic stripe-like pattern after a certain time.


Introduction
This paper discusses instabilities in the case of a viscoelastic fluid. A pure viscous (Newtonian) fluid exposed to external stress converts all the stress energy into heat. A viscoelastic fluid on the other hand is able to store a part of the energy. To model this behavior, we choose the linear Maxwell fluid.
One of the instabilities we want to discuss is the Rayleigh-Taylor instability (RTI). RTI appears if a more dense liquid is located above (below) a less dense one and a volume force pointing downward (upward) exists. Such a constellation is unstable and any small disturbance of the flat interface between the two liquids will grow exponentially in time. The instability is named after Taylor [1] and Strutt, 3rd Baron Rayleigh [2]. The Rayleigh-Taylor instability appears quite often, from the structures which are formed by a collapsing star to the patterns in the billows of cigarette IMA10 -Interfacial Fluid Dynamics and Processes. Guest editors: Rodica Borcia, Sebastian Popescu, Ion Dan Borcia. a e-mail: schoefra@b-tu.de (corresponding author) b e-mail: bestehorn@b-tu.de smoke. RTI of a Newtonian two-layer system is studied in [3].
RTI of viscoelastic drops is examined in [4] using the Oldroyd-B fluid model. Reference [5] discusses gravity and Rayleigh-Taylor waves with a Jeffrey fluid model. In Ref. [6], the linear stability analysis is applied on a Maxwell fluid. Two Maxwell fluid layers are investigated in Ref. [7] in context of tectonics in the upper mantel, while Ref. [8] assumes a Maxwell fluid over a Newtonian fluid to describe large-scale processes in the lithosphere. RTI of linear viscoelastic fluids is investigated by Ref. [9] through a displacement field which is applied to the stress tensor. Linear and weakly nonlinear results are presented.
We shall also consider the Faraday instability, first studied by and named after Faraday [10]. The Faraday instability is caused by an external vertical vibration. The surface of the fluid starts oscillating after some time if a critical amplitude is exceeded. Depending on the vibration amplitude and frequencies, the surface of a vibrating fluid is organized in regular structures like triangles, rectangles, or hexagons. Stability analysis of Faraday waves on viscoelastic fluids is discussed in Ref. [11]. Experimental results of vibrating viscoelastic fluids are shown in Ref. [12,13]. Numerical investigations of Newtonian Faraday instability using a phase field model one can find in [14]. Reference [15] studies the stability of a viscous fluid layer covered with an elastic polymer film under external vibration by usage of Floquet's theorem.
A combination of Faraday and Rayleigh-Taylor instability is studied in [16]. The falling of viscoelastic films down an inclined plane is investigated through long-wave approximation, linear stability analysis, and weakly non-linear approaches by [17].

Governing equations
Using the linear Maxwell fluid, we have two parameters for the description of the viscoelastic behavior. The relation between the stress tensorσ and the shear rate tensorD is given by [18] Λ∂ tσ +σ = 2ηD (1) where η and Λ are the dynamic viscosity and the relaxation time of the fluid. If Λ = 0, Eq. (1) constitutes the inner stress of the Newtonian fluid. The ratio η/Λ is equal to Young's modulus E , so we can find another special case. For η → ∞ and Λ → ∞, but η/Λ = const., the linear Maxwell fluid is equivalent to the elastic solid. We describe conservation of mass and momentum with the incompressible Navier-Stokes equations A two-dimensional sketch of our system is shown in Fig. 1. At z = 0, the fluid is in contact with the solid substrate and no fluid movement is possible On the other hand, the boundary at z = h(x, y, t) is assumed to be a free surface and the tangential stress must disappear. Thust The normal stress balance on the free surface readŝ Here, p 0 is the pressure of the environment, p the fluid pressure under the surface, γ the surface tension, and K the curvature. The free surface is defined by a function h(x , y, t). This function gives the fluid depth at every location and time and its temporal evolution is ruled by the kinematic condition In lateral directions, we assume periodic boundaries with x ∈ [0, L] and y ∈ [0, L]. Equation (1) can be simplified by assuming a viscoelastic force field which leads with the use of Eq. (3) to One characteristic time scale of the system is the diffusion time where ν is the kinematic viscosity. We introduce the dimensionless quantities (bearing tildes) Applying long-wave approximation [19] and the Karman-Pohlhausen method [20], our system (leaving all tildes) simplifies to (see Appendix A) qiqj h . A similar approximation of a thin film is also applied by Ref. [21].
The system variables are the flow rate q = (q x , q y ) T , Note that all dependent variables are now 2D instead of the original 3D fields. In addition to Λ, two new parameters appear, the dimensionless surface tension and the Galileo number G: When we apply harmonic oscillations, the Galileo number becomes time-dependent and reads

Rayleigh-Taylor instability
First, we investigate the instability of Rayleigh-Taylor waves. Due to isotropy in the horizontal plane, it is sufficient to examine the 2D case After linearizing with respect to the small amplitudes u, f , and η, the resulting eigenvalue problem yields the characteristic equation The most dangerous mode has the growth rate max( (λ)). The system is unstable if max( (λ)) > 0.
Only with a negative Galileo number, an instability can be observed. The dependency of max( (λ)) from k can be understand as dispersion relation. We notice that the growth rate has a maximum in k , which is independent of Λ. With increasing Λ, the maximum growth rate increases stagnating. A destabilizing effect of relaxation time was also found by Ref. [22,23]. An increase of the surface tension Γ leads to a decrease of the fastest growing wavenumber in the unstable RTI domain. A decrease of G shifts the unstable domain to the short-wave region of the dispersion spectrum. In  mentioned that viscoelastic solutions have higher viscosities than their solvents; therefore, it is acceptable to assume a viscosity 440 times higher than for water.

Faraday instability
Under harmonic external excitations, we utilize Floquet's theorem to estimate the stability condition of the system. We take The time-dependent components build a vector ξ(t) = (η(t), u(t), f(t)) T , which fulfils the linear system Here,Ā denotes the Jacobian whereĀ(t) =Ā(t + T ) with T = 2π/ω due to Eq. (15). A matrixM which transform the system one period ahead is introduced. Therefore We determine the columns of the monodromy matrix M by integrating for three linearly independent initial conditions i = 1, 2, 3. The vector ξ can be decomposed into where w i (t) = w i (t + T ). The λ i are the Floquet exponents. According to Floquet's theorem, it applies where σ i are the eigenvalues ofM . Due to Eqs. (22,23), max|σ| determines the system stability. We compute σ i by solving Eq. (21) applying a fourth-order Runge-Kutta method. As initial condition, we take an Euclidean base. Figure 3 shows lines of marginal stability (Arnold tongues), max|σ| = 1. We observe successive subharmonic and harmonic unstable solutions. A change of relaxation time causes also a change of the position of the marginal lines in amplitude and wavenumber.
There is a critical amplitude where the system becomes unstable. This amplitude decreases with increasing frequency. The decrease in frequency is faster for higher relaxation times. The wavenumber of this amplitude is denoted with k crit . In Fig. 4, the critical wavenumber is plotted over the frequency. The k crit (ω)-curves show a discontinuity, when an Arnold tongue replaces another one as the most dominant in the system. In the region where Λω ≈ 1, we observe a harmonic tongue as critical. Also, the size of the appointed frequency window (H 2 in (a) and H 1 in (b)) increases compared to the Newtonian fluid (c) in Fig. 4. Reference [11] also found a harmonic resonance in that region and a minimum for the critical amplitude.

Numerical results
Applying a standard finite difference method [24], solutions in two (2D) and in three dimensions (3D) have been determined. The code has been developed in C++. In Fig. 5, a 2D Rayleigh-Taylor wave is shown for different time steps. There are two resonances, one at x ≈ 27 and one at x ≈ 66. Two perturbations with nearly the same amplitudes merge at t = 40 and form a new bigger RT wave crest at t = 50. After its formation, it begins to shrink as quickly as it has been formed, and at t = 54, the drop is not much bigger than the others. This swapping back behavior is a hint to the elastic property of the film. The second resonance, at x ≈ 66, grows to such an extent that it finally causes rupture. Figure 6 shows an RTI in 3D. Surface perturbations grow over time. The flow field has sinks where the surface height increases and sources where it decreases. Average flow rates and surface perturbations increase in time.
The 2D Faraday waves of a Newtonian fluid are shown in Fig. 7. There is a π-phase shift after one period; the structures are subharmonic. The wavenumber fits to the results from Floquet analysis, where we expect to be in the 1 S and 1 H Arnold tongue region. Therefore, the harmonic solution is not visible in the numerics.
In Fig. 8, the height for a viscoelastic fluid with Λ = 0.5 is illustrated. We observe a fundamental difference between the patterns of the Newtonian fluid and the linear Maxwell fluid. In the beginning, the Maxwell fluid acts quite similar to the Newtonian fluid. It takes some time before the structures of the surface appear, and then, they look quite similar to their viscous equivalent. After some time, the repeating structures begin to change. The surface oscillates almost chaotically, a full phase repeats only almost the same pattern. The explanation to this can be that the linear Maxwell fluid has a memory of every stress to which it was exposed. This phenomenon is called stress relaxation, it is the major difference between the Newtonian fluid and the linear Maxwell fluid and the best explanation for the appearance of these strange surface patterns. A more detailed presentation of stress relaxation can be found in [18].
The Faraday instability of Newtonian fluids in 3D produces surface patterns like squares or hexagons [21,25]. In Fig. 9, we can see how the 3D viscoelastic film develops from the random perturbation in the beginning to a subharmonic stripe structure. At t = 10, we see a random surface, which becomes a subharmonic square structure at t = 124. The squares repeat periodically for some time with small changes. At t = 208, we observe that some of the former squares have developed into diamonds and connect with their neighbors. Then, at t = 228, the regular structure begins to collapse. It follows a kind of transition phase where irregularly changing structures occur, similar to the almost chaotic surface patterns we also observe in the 2D calculation (t = 252). After the transition phase (t = 320), the surface ends up in a stripe pattern. This is a similar behavior as the corresponding 2D solution. The applied scaling is G = Γ = 40. This can be associated with ρ 0 = 789 kg/m 3 , γ = 0.022 N/m, ν = 3.4 · 10 −5 m 2 /s and d = 1.67 mm (e.g., ethanol-based polymer solution). Vibrating with 12 Hz (ω = π) or 24 Hz (ω = 2π).
Reference [13] studied Faraday waves of viscoelastic fluids experimentally. The authors use a 1.5% polyacrylamid-co-acid solvent. A similar transition of hexagonal patterns to a stripe-like structure is observed.

Conclusions
The Rayleigh-Taylor and the Faraday instability for a viscoelastic Maxwell fluid are examined applying a linear stability analysis as well as direct numerical simulations of a recently derived long-wave model. The results of the numerical solution and the linear analysis confirm each other. The wave numbers and frequencies predicted by the linear model are reproduced by the numerical approach.
The linear stability analysis has shown that the highest growth rate λ max of RTI increases with Λ. For high values of Λ, the increase stagnates. Here, it should be mentioned that in practice, it is very difficult to increase a fluids' relaxation time without increasing its viscosity too. The corresponding k max is not influenced by the relaxation time. Numerical solutions of 2D RTI waves show elastic behavior and resonance.
The Faraday instability in 2D and 3D creates erratic surface patterns. The 3D simulation shows a subharmonic stripe structure which does not occur in Newtonian fluids.
Further research can focus on the implementation of various viscoelastic models and can clarify how pattern formation in these materials distinguish from each other. Also, the Faraday instability of viscoelastics offers many possibilities for further research. For example, the characterization of the patterns which are not observable in Newtonian fluids or a detailed study of the resonance behavior is conceivable. High-frequency regions are not considered by our model due to their short-wave nature.
The long-wave model we developed here can be extended in several directions. One can add horizontal vibrations, resulting in another kind of Faraday instability [27].
Certain combinations of instabilities are possible too. One can combine RTI and Faraday instability or consider the Faraday instability on an inclined surface. The study of instabilities and pattern formation in viscoelastic fluids is a complex topic which offers many possibilities for future research.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by Projekt DEAL.

Data availability
The datasets generated during and analysed during the current study are available from the corresponding author on reasonable request.