Viscoelasticity-Induced Instability in Plane Couette Flow at Very Low Reynolds Number

Abstract

Elasto-inertial turbulence (EIT), a new turbulent state found in polymer solutions with viscoelastic properties, is associated with drag-reduced turbulence. However, the relationship between EIT and drag-reduced turbulence is not currently well-understood, and it is important to elucidate the mechanism of the transition to EIT. The instability of viscoelastic fluids has been studied in a canonical wall-bounded shear flow to investigate the transition process of EIT. In this study, we numerically deduced that an instability occurs in the linearly stable viscoelastic plane Couette flow for lower Reynolds numbers, at which a non-linear unstable solution exists. Under instability, the flow structure is elongated in the spanwise direction and regularly arranged in the streamwise direction, which is a characteristic structure of EIT. The regularity of the flow structure depends on the Weissenberg number, which represents the strength of elasticity; the structure becomes disordered under high Weissenberg numbers. In the energy spectrum of velocity fluctuations, a steep decay law of the structure’s scale towards a small scale is observed, and this can be recognized as a ubiquitous feature of EIT. The existence of instability in viscoelastic plane Couette flow supports the idea that the transitional path toward EIT may be mediated by subcritical instability.

1. Introduction

Viscoelastic fluids are easily produced by adding a small amount of long-chain polymers to a Newtonian solvent, such as water. These fluids exhibit a property by which the skin friction in wall-bounded turbulence is significantly reduced [1]. Many researchers have experimentally and numerically studied the mechanism of drag reduction [2,3,4,5].

Certain aspects of the drag-reduction mechanism are generally accepted [6,7]. Polymer molecules move randomly under Brownian motion caused by collisions among the molecules, resulting in an equilibrium state called a random coil. In wall-bounded turbulent flow, polymers are stretched along the flow direction due to the strong shear stress near the wall. The force to be relaxed when tending towards the equilibrium state is associated with the elasticity of the polymeric fluid. The time scale of polymer relaxation is called the relaxation time, and the viscoelastic flow state is characterized by the Weissenberg number $Wi$, which is the ratio of the relaxation time to a time scale of shear flow. The elastic forces experienced by the stretched polymers near the wall suppress streamwise vortices [8,9,10]. The effect of viscoelasticity that modifies a sustaining process of Newtonian turbulence has been observed in a relatively low-$Wi$ flow, which leads to drag-reduced turbulence. In this flow state, the turbulent drag reduction rate increases with the polymer concentration. It is well-established that an increasing polymer concentration would manage to approach the maximum drag reduction (MDR) asymptote, as predicted by Virk [3]. The MDR is the limiting state in which the dynamics of turbulence are weakened and turbulence can be barely sustained. In higher concentrations of a solution, viscoelasticity is known to cause a flow instability that is significantly different from the Newtonian turbulence. Advanced numerical simulations in recent years have succeeded in capturing such instabilities, and have found a flow dominated by structures arranged in the spanwise direction, unlike the familiar turbulence dominated by streamwise structures, such as streak structures near the wall [11,12]. This flow state has been termed elasto-inertia turbulence (EIT), as it is not only triggered by a purely elastic instability, but additionally requires inertial instability. The EIT state, the existence of which has been established predominantly by numerical simulations, has been experimentally confirmed in the past few years, and it has been agreed to be an essential feature of viscoelasticity.

The pipe flow of polymer solutions may undergo a transition that is different to that of a Newtonian fluid [13,14,15]. However, the cause for this remains unclear, even when considering both the resultsthat the transition is triggered earlier than in a Newtonian fluid and that the transition is delayed. Samanta et al. [11] first demonstrated the flow regimes of viscoelastic pipe flows at high shear rates using tubes with smaller diameters than those used in previous pipe experiments. This achievement revealed that viscoelasticity triggers turbulence below the critical Reynolds number, Re, of the Newtonian fluid flow, and modifies the conventional Newtonian laminar-turbulent transition at high $Wi$ values. They also found that the friction factor of the EIT state follows the friction law of the MDR asymptote. Since then it has been hypothesized that EIT is associated with the MDR asymptote, and results supporting this association have been demonstrated. To investigate the association between the MDR asymptote and EIT, Choueiri et al. [16] focused on the path to EIT instability using the modified laminar-turbulent transition in a viscoelastic fluid. Their flow state diagram of the viscoelastic fluid, drawn based on the $Re$ and the concentration of the polymer plane, systematically shows the transition phenomena obtained from previous studies. They found that increasing the polymer concentration at a fixed, relatively low Reynolds number resulted in a transition from Newtonian turbulence to laminar flow. By further increasing the concentration while keeping the Reynolds number fixed, the transition to the MDR regime occurred regardless of the initial state of the Newtonian fluid. Some studies have focused on drag-reduced flow and have considered this separately from EIT. In other words, they regarded drag-reduced flow as a flow state that originates from Newtonian turbulence and that is modified according to the increase in either concentration (in experiments) or $Wi$ (in simulations). However, the pipe experiments conducted by Choueiri et al. [16] and the numerical simulations conducted by Lopez et al. [17], in support of their experimental work, confirmed that the MDR asymptote, as the final state of the drag-reduced flow, undergoes a re-entrant transition from the laminar flow. According to a series of studies by Hof’s group [11,16,17], the MDR asymptote originates from an instability that is different to that observed in Newtonian turbulence, and that the dynamics of the MDR asymptote have characteristic features in common with EIT. As summarized here, these results, in the context of the transition to EIT, demonstrate that drag-reduced flow is a phenomenon wherein Newtonian turbulence and EIT coexist, and in the transition under a constant Reynolds number, Newtonian turbulence is weakened as the concentration increases, and eventually is replaced by EIT. The MDR asymptote exists as a boundary of the marginal flow state, beyond which Newtonian turbulence cannot be maintained during this scenario.

Garg et al. [18] investigated the linear stability of viscoelastic pipe flow, finding it to be linearly unstable. This instability gives rise to the early transition of pipe flow at the significantly lower critical Reynolds number of viscoelastic pipe flow than that of Newtonian pipe flow, and this is considered as the first stage of EIT. This instability has additionally been confirmed in plane Poiseuille flow. Some advanced stability analyses of viscoelastic flow have revealed several unstable solutions and demonstrated that elastic instability can occur without inertia in the EIT, and appears with an infinite elasto-inertial hierarchy [19,20,21,22]. This suggests that viscoelastic instability could develop through a tangled nonlinear process. The transition process of viscoelastic flow has not yet been clarified [19]. A nonlinear mechanism of EIT has been proposed by Morozov and van Saarloos [23] using a simple model. Their model was based on elastic instability caused by curved streamlines, associated with purely elastic instability [24,25]. Using nonlinear stability analysis, they predicted that Couette and Poiseuille flows may undergo a subcritical transition to EIT [23,26]. They also found that the subcritical instability of viscoelastic flow has a structure that is very similar to that of EIT, as confirmed in other studies [12,20,27].

The stability of viscoelastic Couette flow has been investigated for models of various constitutive equations using several methods [28,29,30,31,32]. There is a consensus that the widely-used models (upper convective Maxwell, Oldroyd-B, Giesekus, and the FENE-P model) demonstrate the linear stability of viscoelastic Couette flow. Therefore, if the transition to EIT occurs in Couette flow, the instability grows through a nonlinear process. Obtaining a nonlinear solution of viscoelastic Couette flow is important for understanding the subcritical transition of viscoelastic fluids. This nonlinear solution provides an initial template to obtain unstable eigenfunctions that are difficult to predict. Furthermore, unlike pipe and Poiseuille flow, the effect of the absence of a linear instability solution on the transitional process of viscoelastic fluids can be inferred. In viscoelastic flow, especially in Couette flow, there are few studies focusing on the subcritical transition, and the full picture of the transition is as yet unknown. In this paper, we investigate the mechanism of the transition of the viscoelastic Couette flow to EIT, and the characteristics of EIT using three-dimensional direct numerical simulation. The current numerical simulation method for EIT still has many problems and needs to be improved.

2. Numerical Method and Procedure

Our three-dimensional direct numerical simulation was performed under the incompressibility condition, coupled with a constitutive equation. The flow system used is a wall-bounded shear flow, i.e., a plane Couette flow driven by moving the upper and lower walls, as shown in Figure 1. The streamwise, wall-normal, and spanwise directions are denoted as x, y, and z, respectively. The non-slip boundary condition at the walls and the periodic boundary conditions in x and z are applied. The dimensional continuity and momentum equations are, respectively, given by

$$\nabla ·\mathbf{u}=0,$$

(1)

$$\rho \{\frac{\partial \mathbf{u}}{\partial t}+\left(\mathbf{u}·\nabla \right)\mathbf{u}\}=-\nabla p+{\mu }_0{\nabla }^2\mathbf{u}+\nabla ·\tau ,$$

(2)

where $\rho$, $\mathbf{u}$, p, and $\tau$ represent the density of the fluid, a velocity vector, pressure, and the polymeric stress tensor, respectively. These equations may be scaled with the wall speed $U_w$ and the half-channel height $\delta$. Regarding the polymeric stress tensor, scaling with ${\mu }_p{U}_w/\delta$ allows us to obtain a non-dimensional momentum equation as

$$\frac{\partial {\mathbf{u}}^{*}}{\partial t^{*}}+\left({\mathbf{u}}^{*}·\nabla \right){\mathbf{u}}^{*}=-\nabla p^{*}+\frac{\beta }{Re}{\nabla }^2{\mathbf{u}}^{*}+\frac{1-\beta }{Re}\nabla ·{\tau }^{*}.$$

(3)

Here, the viscosity ratio $\beta$ is defined as the ratio of the solvent viscosity to that of the solution $\beta ={\mu }_s/{\mu }_0={\mu }_s/({\mu }_s+{\mu }_p)$, where ${\mu }_s$ and ${\mu }_p$ are the viscosity of the solvent and polymer contributions, respectively, and ${\mu }_0$ represents the zero-shear viscosity of the solution. The superscript ‘*’ indicates the scaled value within the equations. We define two dimensionless parameters, the Reynolds number $Re=U_w\delta /\nu$, and the Weissenberg number $Wi=\lambda U_w/\delta$, where $\nu$ is the kinematic viscosity and $\lambda$ is the elastic relaxation time. The polymeric stress tensor $\tau$ is solved by means of the constitutive equation, coupled with the continuity and momentum equations. The Giesekus model [33] for the constitutive equation is employed and given by

$$\tau +\lambda \stackrel{▽}{\tau }+\frac{\alpha }{{\mu }_p}\tau ·\tau =2{\mu }_{p}D,$$

(4)

where D is the deformation velocity tensor, $\alpha$ is the non-dimensional mobility factor, and $\stackrel{▽}{\tau }$ represents the upper convective derivative of the polymeric stress, which can be written as

$$\stackrel{▽}{\tau }=\frac{\partial \tau }{\partial t}+\left(\mathbf{u}·\nabla \right)\tau -(\nabla \mathbf{u})·\tau -{(\tau ·\nabla \mathbf{u})}^{\mathrm{T}}.$$

(5)

Moreover, we introduce the dimensionless conformation tensor $c_{ij}$, as follows

$${\tau }_{ij}=\frac{{\mu }_p}{\lambda }(c_{ij}-{\delta }_{ij}),$$

(6)

where ${\delta }_{ij}$ is the Kronecker delta. In the actual process of DNS, the following governing equation of $c_{ij}$ is solved numerically:

$$\frac{\partial c_{ij}}{\partial t^{*}}+u_k^{*}\frac{\partial c_{ij}}{\partial x_k^{*}}=c_{ik}\frac{\partial u_j^{*}}{\partial x_k^{*}}+c_{kj}\frac{\partial u_i^{*}}{\partial x_k^{*}}-\frac{1}{Wi}\{c_{ij}-{\delta }_{ij}+\alpha (c_{ik}-{\delta }_{ik})(c_{kj}-{\delta }_{kj})\}.$$

(7)

We follow the same numerical method used by Nimura et al. [34,35]. The fourth-order central finite difference scheme is used in the x and z directions, whereas the second-order central scheme in y. For the time integration, the second-order Crank–Nicolson scheme and the second-order Adams–Bashforth scheme are employed for the wall-normal viscous term and the other terms, respectively. The fractional step method is adopted to couple Equations (1) and (3). The pressure Poisson equation is solved using the tridiagonal matrix algorithm (TDMA) in y and the fast Fourier transform (FFT) in x and z, without the iterative method. This numerical procedure has been fully established for the Newtonian turbulent channel flow by Abe et al. [36]. As for the constitutive equations with the Giesekus model, a flux limiter of the MINMOD scheme is adopted to approximate the spatial derivatives in the advective terms; thus, an artificial diffusive term is not included, as first introduced by Yu and Kawaguchi [37] and commonly used even in recent works [38,39,40,41]. The periodic boundary conditions are applied in x and z: $\mathbf{u}(L_x,y,z)=\mathbf{u}(0,y,z)$ and $\mathbf{u}(x,y,L_z)=\mathbf{u}(x,y,0)$. The moving wall surfaces are non-slip: $\mathbf{u}=(-U_w,0,0)$ at $y=0$, and $\mathbf{u}=(U_w,0,0)$ at $y=2\delta$.

For the initial condition of the velocity field, we induce a moderately weak random disturbance, applied to the entire channel, to the linear solution of the plane Couette flow; thus, $u=U_w(y/\delta -1)$. The magnitude of the disturbance amplitude is set sufficiently small for only the main mode to grow after the initial disturbance decays rapidly. We preliminarily confirmed that the growing mode is consistent in such a range of sufficiently small magnitudes of amplitude. It was also confirmed that the mode of growth after the decay of the disturbance is consistent within a sufficiently small range of the initial amplitude. In the present study, the Reynolds number is fixed at either 10 or 100, and some rheological parameters are set as $\alpha =0.001$ and $\beta =0.9$. Regarding $\alpha$, we referred to the values in the literature [37,42], which were confirmed to be in good agreement with the experimental counterparts. The additive concentration is assumed to be diluted.

Figure 2 shows the configuration of numerical grids used consistently for all cases in this study. The size of the domain, $(L_x\times L_y\times L_z)/{\delta }^3=12.8\times 2.0\times 6.4$, and the number of grid points $(N_x\times N_y\times N_z)=256\times 96\times 64$ are kept fixed for all cases, which lead to streamwise/spanwise grid resolutions of $\mathrm{\Delta }{x}^{*}=0.05$ and $\mathrm{\Delta }{z}^{*}=0.1$, respectively. The non-uniform grids with $\mathrm{\Delta }{y}^{*}=0.017$–0.023 are used in y. For a gridconvergence test, see Appendix A. Within our simulations, we simulated the viscoelastic flow field where the elasticity was increased with constant inertia, or more specifically, the Reynolds number was fixed and the Weissenberg number was increased to systematically investigate the $Wi$-dependence of elasto-inertial instability. At high Weissenberg numbers, it is difficult to numerically simulate the viscoelastic flow due to the numerical instability of the constitutive equation caused by the absence of a diffusion term. It is often extremely difficult to accomplish a simulation with the desired high-$Wi$ regimes in a steady state. Therefore, we gradually increased the Weissenberg number, and demonstrated the results up to the Weissenberg number that developed the flow into a steady state without numerically diverging when sufficient time had elapsed.

The flow states presented in Section 3 were statistically steady in terms of the time series of a spatially averaged kinetic energy of the velocity fluctuations. The kinetic energy K is a value integrated in the whole domain, as follows:

$$K=\frac{1}{L_x{L}_y{L}_z}\int \int \int \frac{u_i^{\prime }}{2}\mathrm{d}x\mathrm{d}y\mathrm{d}z,$$

(8)

where $u^{\prime }$ represents the fluctuation value from the mean velocity, described below in Section 3.1. Figure 3 shows the time series of K for cases that exhibit the onsets of instabilities. The growth of an instability can be seen in the time evolution from $t=0$, at which a disturbance is induced in the laminar flow.In all cases, K finally converges and the flow state reaches a statistically-steady state after $t{U}_w/\delta =2000$. The flow fields shown in Section 3 are the typical ones in the steady state, obtained after the artificial disturbance was attenuated and the unstable mode was developed well.

We varied the Weissenberg number while fixing two Reynolds numbers, $Re=10$ and 100, and investigated whether viscoelastic instability appeared. These Reynolds numbers are below the $R{e}_c\approx 324$ obtained by Duguet et al. [43] as the threshold to sustain turbulence of the Newtonian Couette flow. These are also lower than the region of existence of exact coherent structures (ECSs) for the plane Couette flow ($Re\ge 127.7$), which is known to be a non-linear unstable solution [44,45,46]. Thus, if any instability appears in a flow field below these Reynolds numbers, it is considered to be due to viscoelasticity.

3. Results and Discussion

3.1. Onset of a Viscoelasticity-Induced Instability

Figure 4 shows the profiles of the mean velocity $\langle \overline{u}\rangle$, which change as the Weissenberg number increases for each Reynolds number; enlarged views of the center of the channel are inserted. The operation $\langle ·\rangle$ represents the time average and $\overline{·}$ is the spatial average on the x–z plane parallel to the walls. These profiles are almost the same as the laminar solution of Couette flow throughout the channel. The velocity gradient, i.e., the shear rate, decreases locally at the channel center for $Wi\ge 20$, as shown in Figure 5. It is clear that the velocity profile forms an inflection point at the channel center or at three points around it. For $Wi=10$ at both Reynolds numbers, the mean velocity profile has no inflection point, meaning instabilities cannot occur in the flow fields and the laminar flow is maintained. An instability can be observed in the flow fields for $Wi\ge 20$ at both Reynolds numbers, where inflection points are formed in the velocity profile. As will be shown later, spanwise-oriented structures are formed in the presence of this instability. At the high Weissenberg number $Wi=1000$ for $Re=100$ in Figure 5b, four inflection points are present around the center of channel. However, there is hardly any difference in the flow structure caused by the additional inflection points in the velocity profile. In the case of $Re=10$, for high Weissenberg numbers of $Wi>30$, the growth of the spanwise-oriented structures was confirmed, but this simulation diverged numerically and the steady state could not be obtained.

Figure 6 shows the spatio-temporal averaged rms (root-mean-square) of velocity fluctuations. It is clear that the intensity of velocity fluctuations in either direction increases with an increasing Weissenberg number. This tendency was confirmed for both $Re=10$ and $Re=100$, and the profiles of both numbers are similar. It can be estimated that there is no dependence on the Reynolds number within this range. These RMS velocity profiles are almost symmetric profiles, centered on $y/\delta =1$, and the symmetry of the flow does not break above and below the channel center. There is a sharp peak in the intensities of velocity fluctuations in the streamwise direction $u_{\mathrm{rms}}$ at the center. The intensities in other directions, $v_{\mathrm{rms}}$ and $w_{\mathrm{rms}}$, have a peak around the center, but weaken locally at $y/\delta =1$. These profiles deviate significantly from inertia-dominated turbulence.

Figure 7 shows the flow structure observed as a fully developed state at $Re=10$ in the x–z plane in the channel center $y/\delta =0$. The contour represents the fluctuating velocity $v^{\prime }$ in the wall-normal direction, and the vectors represent the fluctuating velocity $u^{\prime }$ and $w^{\prime }$ in the plane. Here, the fluctuating velocity $u_i^{\prime }$ is defined as $u_i^{\prime }=u_i-\overline{u_i}$; however, $\overline{v}$ and $\overline{w}$ are nominally zero. Hence, $v^{\prime }\approx v$ and $w^{\prime }\approx w$. It can be seen that positive and negative regions of $v^{\prime }$ appear alternately, as shown in Figure 7, and elongated spanwise-oriented structures are shown to locally align in the streamwise direction. Such structures have been observed in other systems, such as Poiseuille flow [11,12] and pipe flow [17,18]. In pressure-driven flows, spanwise-oriented structures are formed near the wall, whereas in the current Couette flow, the flow structures visualized in Figure 7 are concentrated in the channel center. These structures are simpler for lower Weissenberg numbers at the same Reynolds number. As the Weissenberg number increases, the intensity of velocity fluctuations increases, the regularity of the flow structure is lost, and three-dimensional turbulent structures are formed. Such Weissenberg-number-dependence of structures shows the same tendency for a different Reynolds number, specifically, $Re=100$. The dependence on the Reynolds number can be observed in the distribution of $v^{\prime }$. At a high Reynolds number of $Re=100$, a large-scale flow with a spatial size of >2$\delta$ in x and z is formed, as shown in Figure 8. The large-scale flow is not obvious at a low Reynolds number $Re=10$ within the Weissenberg number range of $10\le Wi\le 30$ that was simulated without the numerical divergence. It can be observed that the formation of the flow is not a feature of the high Weissenberg flow, but rather a feature of the flow at a moderately high Reynolds number, because large-scale flows also appear in the flows of $Re=100$ and $Wi=20$ (see Figure 8). The large-scale flow is described in detail in the next section.

3.2. Scale of EIT Structures

To investigate a statistical property of the spanwise-oriented structures, the two-point correlation coefficient of velocities were calculated. The two-point correlation coefficient is defined as the correlation of the velocity fluctuations for each direction at a fixed y, and is normalized with the RMS values of velocity fluctuations. It can be written as follows:

$$\begin{array}{c}\hfill R_{u_i{u}_j}(\mathrm{\Delta }x,y)=〈{\displaystyle \frac{\overline{u_i^{\prime }(x,y,z)·u_j^{\prime }(x+\mathrm{\Delta }x,y,z)}}{u_{i,\mathrm{rms}}^{\prime }·u_{j,\mathrm{rms}}^{\prime }}}〉,\end{array}$$

(9)

$$\begin{array}{c}\hfill R_{u_i{u}_j}(\mathrm{\Delta }z,y)=〈{\displaystyle \frac{\overline{u_i^{\prime }(x,y,z)·u_j^{\prime }(x,y,z+\mathrm{\Delta }z)}}{u_{i,\mathrm{rms}}^{\prime }·u_{j,\mathrm{rms}}^{\prime }}}〉.\end{array}$$

(10)

The two-point correlation coefficients for $Re=10$, as functions of the streamwise or spanwise distances where $\mathrm{\Delta }x$ and $\mathrm{\Delta }z$ form a reference point, are plotted in Figure 9. These correlation coefficients are averaged over time and space in the channel’s center plane. The streamwise wavelength of the spanwise-oriented structure coincides with the first positive peak that appears after $R_{vv}$ drops negatively from unity. According to Figure 9a, the wavelengths in the streamwise direction are ${\lambda }_x=0.65\delta$ for $Wi=20$ and ${\lambda }_x=0.70\delta$ for $Wi=30$. The wavelength of the spanwise-oriented structure ${\lambda }_x$ tends to increase as the Weissenberg number increases. The correlation coefficient of the streamwise velocity $R_{uu}$ also has a peak at the same position; however, the value of $R_{uu}$ does not drop to negative values. It is presumed that this is because the velocity does not have sufficient strength to form a vortical motion. The structure of EIT found by Dubief et al. [12] also highlighted that spanwise-oriented structures are not sufficiently strong to form vortices. A negative peak is observed at $\mathrm{\Delta }x\approx 3$–4$\delta$ for $R_{uu}$ and $R_{ww}$, which indicates a wavelength larger than the wavelength of the spanwise-oriented structure. This corresponds to the wavelength of the large-scale flow that is slightly evident in the flow field plotted in Figure 7. The wavelength of a large-scale flow cannot be estimated using the present simulations, because only the approximate two wavelengths are captured in the computational domain in the streamwise direction. In Figure 9b, the two-point correlation coefficient as a function of the spanwise distance $\mathrm{\Delta }z$ indicates that there is no small-scale structure in the spanwise direction, and demonstrates that the large-scale flow is affected by the periodical condition, as well as the streamwise direction. Therefore, in order to capture the large-scale flow induced by viscoelastic instability, simulations are required in the calculation domain, expanded in the streamwise and spanwise directions. The Reynolds-number-dependence of flow structures can be observed in $R_{vv}$ for $Re=100$, as shown in Figure 9 and Figure 10. At $Re=10$, a small-scale structure causes $R_{vv}$ to negatively drop away from unity, whereas at $Re=100$, a large-scale flow causes $R_{vv}$ to negatively drop at a large distance. It is therefore considered that the scale of the flow that dominates the flow field has changed from a small-scale structure to a large-scale flow, and this large-scale flow can be seen to clearly appear within Figure 8. This does not necessarily mean that a high Reynolds number does not give rise to the formation of spanwise-oriented structures. Considering the correlation of the flow fields separated in each scale, drops from positive to negative occur in each case. The wavelength of spanwise-oriented structures is defined as the distance of the first positive peak of $R_{vv}$ after dropping from unity, even if it does not change from positive to negative. The wavelength of spanwise-oriented structures depends on the Reynolds number. According to Figure 10, the wavelength in the streamwise direction is ${\lambda }_x=0.50$ for both $Wi=20$ and $Wi=100$ cases at $Re=100$. For $Re=10$ and $Re=100$, at the same Weissenberg number $Wi=20$, the wavelength when $Re=100$ is smaller than that when $Re=10$. It is presumed that the dependence on the Reynolds number is larger than the dependence on the Weissenberg number. The wavelength of a large-scale flow with a higher Reynolds number $Re=100$ cannot be estimated because it is also affected by the size of the computational domain, and the domain is insufficient to capture a large-scale flow that is free from the periodic boundary.

By estimating the wavelength of spanwise-oriented structures using the correlation coefficient as the threshold value, the filter operation is applied to the fluctuating velocity field, and spanwise-oriented structures and the large-scale flow are separated. The small- and large-scale velocity fields are described as $u_i^s$ and $u_i^l$, respectively, and the fluctuating velocity decomposition is written as $u_i^{\prime }=u_i^s+u_i^l$, where $u_i^{\prime }$ is the total velocity field before decomposition. They are strictly separated at the threshold in Fourier space and the cross-correlation between the small- and large-scale velocity fields is zero. Note that the large-scale component of v is almost zero, that is, $v^{\prime }=v^s$. A three-dimensional visualization is used to show the flow structures of each scale based on the decomposed velocity fields in Figure 11 and Figure 12. In these figures, the threshold of the isosurface is the same in all cases. The width of the visualized grid mesh is the same in all directions and is equal to $\delta$, but does not correspond to the computational grid width of DNS. Figure 11 shows snapshots of a typical flow field for $Wi=20$ and 30 at $Re=10$. The red and blue isosurfaces represent the positive and negative values of the second invariant of the velocity gradient tensor

$$Q=\frac{1}{2}\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i},$$

(11)

i.e., the regions of rotational motion and extensional/compressional motion, respectively. Each region is arranged alternately in the streamwise direction, and the sizes are approximately equal to each other. Their size matches the wavelength of the small-scale structure, confirmed by the two-point correlation coefficients. The regions of the second invariant represented by these isosurfaces are almost unchanged when drawn as the second invariant of the small-scale velocity field, and thus most of the displayed region of topological structures is contributed by the small-scale velocity field. These spanwise-oriented structures displayed by the isosurfaces of Q represent the characteristic structure of the EIT demonstrated by Samanta et al. [11], Dubief et al. [12] and Lopez et al. [17] in other flow systems. Moreover, the large-scale flows $u^l$ are shown by superimposing spanwise-oriented structures with black and white translucent regions, and black and white regions indicate positive and negative values of $u^l$, respectively. The positive and negative structures of $u^{\prime }$, elongated in the streamwise direction, are similar to the ECS of viscoelastic Couette flow, the existence of which was found by Stone et al. [10]. These three-dimensional visualizations provide confirmation that the dynamics of the large-scale flow and spanwise-oriented structure are closely related. The spanwise-oriented structure is a structure in which a straight vortex along the spanwise axis is bent locally. These vortex dynamics are reminiscent of the hairpin vortex found in wall turbulence. According to this dynamic mechanism, a high-speed region is generated along the bent direction, and a large-scale flow is considered to be formed. A large-scale flow with a positive value is observed at a position bent in the positive direction in the streamwise direction, and the other is observed with a negative value at a position bent in the opposite direction. Considering that a large-scale flow is rarely formed in two-dimensional flow, the large-scale flow is regarded to be a secondary flow induced by the deformation of the spanwise structure. As the spanwise structure is deformed at the position where a large-scale flow exists, the three-dimensionality of EIT and the development of the large-scale flow are closely correlated.

As mentioned earlier and as shown in Figure 9 and Figure 10, the streamwise wavelength of spanwise-oriented structures is dependent on the Reynolds number; the wavelength decreases as the Reynolds number increases. Figure 12 shows snapshots of the flow field of $Re=100$ with a varying Weissenberg number, and it can be seen that slightly smaller spanwise-oriented structures than those in in Figure 11 are observed. Large-scale flows develop as the Weissenberg number increases, and this leads to the formation of a three-dimensional turbulence. At $Re=100$ and $Wi=20$, nearly two-dimensional structures are arranged regularly and the intensity of the large-scale flow is weaker. There is no qualitative difference in flow structure for high Weissenberg numbers of $Wi\ge 100$.

Figure 13 shows a magnified visualization of the small-scale structure located in the center of the channel gap. Figure 13a visualizes a typical velocity field with the vectors of in-plane flow $(u^s,v^s)$ at an arbitrary z, which are superimposed on a color contour of the spanwise vorticity ${\omega }_z^s$, based on the small-scale velocity $u^s$. A positive ${\omega }_z^s$ corresponds to a counter-clockwise vortical motion and vice versa. The region where the small-scale velocity exhibitsa vortical motion rather than a shear motion can be also indicated by a negative value of the $Q^s$ criterion, based on the small-scale velocity components, as in Figure 12. Those vortices appear in staggered rows above and below the center of the channel. The region of high ${\omega }_z^s$, except for the layer ($0.98<y<1.02$) appearing in the center of the channel shown in (a), is consistent with the regions with negative $Q^s$ as plotted in (b). It can be noted that this layer has a sufficient possibility of being a numerical phenomenon, because the thickness of this layer corresponds to the length of the two grids in the y direction. However, the row of vortices formed across the layers closely resemble the thin layered structure observed in the EIT of plane Poiseuille flow, as demonstrated by [47]. The difference is that the EIT in plane Poiseuille flow has a thin sheet structure that is elongated in the top-right and bottom-left directions in the x–y plane and is formed near the wall, whereas in plane Couette flow, the row of vortices is unchanged with respect to the mean flow and appears in the center of the channel. Figure 13c visualizes the viscoelastic forces in the same plane and instance in time with (a) and (b), and the contour shows the trace of the conformation tensor $\mathrm{tr}\left(c_{ij}\right)$, which indicates the degree of polymer extension. A large polymer extension occurs at the position between the two vortices that are located above and below the center line on the sheet structure, and a strong viscoelastic force is acting on the fluid. This force acts in opposite directions across the center line, strengthening each of the upper and lower spanwise vortices. We can therefore infer that the viscoelastic force produces the vortices appearing in the small-scale velocity field. Moreover, the viscoelastic force acting continuously around the center of the channel is the cause of the decreasing velocity gradient at the channel center, as shown in Figure 4 and Figure 5.

We additionally confirm the visualization of the large-scale flow in the x–y plane. Figure 14 shows the large-scale flow in the decomposed velocity field $(u^l,v^l,w^l)$. The large-scale flow exhibits a structure that covers the entire channel gap and has an interval longer than that of the small-scale structure visualized in Figure 13. The small-scale structure and large-scale flow are less correlated because their wavelengths are significantly different, and the large-scale flow is slightly disordered with respect to the regularly arranged small-scale structure.

The eigenvalues of the conformation tensor $c_{ij}$ represent the principal stretches of the polymer. The $3\times 3$ tensor of $c_{ij}$ has three eigenvalues. Let us now denote them as ${\sigma }_1$, ${\sigma }_2$, and ${\sigma }_3$. The magnitude of the eigenvalue physically indicates the degree of polymer deformation along the corresponding principal axis. From the viewpoint of thermal equilibrium, the current eigenvalues should be positive, and the symmetric matrix of the tensor $c_{ij}$ is a positive-definite matrix [48,49]. These characters lead to “${\sigma }_1,{\sigma }_2,{\sigma }_3>0$”, which has been satisfied in the whole domain developing in time, according to the present DNS. This fact that our obtained $min\left({\sigma }_i\right)\approx 0$, but was never $\le 0$, clearly demonstrates that the current simulated $c_{ij}$ are strictly positive-definite matrices. Note further that $0<{\sigma }_i<1$ represents a compression of the polymer, whereas $1<{\sigma }_i$ indicates a stretch. Some typical wall-parallel distributions of the three eigenvalues for $Re=10$ and $Wi=20$ at the channel center are shown in Figure 15. Although the current set of eigenvalues is not sorted with respect to their magnitudes, on average ${\sigma }_3>{\sigma }_2>{\sigma }_1$. For instance, the time average of the lowest value of each eigenvalue reveals a relation of $min\left({\sigma }_3\right)\gg min\left({\sigma }_2\right)\approx 1>min\left({\sigma }_1\right)$. In Figure 15, ${\sigma }_3$ clearly exhibits a distribution similar to the small-scale structure, observed via the flow visualization of Figure 11. It can be inferred that the main stretch of the polymer is caused by the small scale. Other-directional stretching is extremely low compared to that in the principal axis relevant to ${\sigma }_3$. As for ${\sigma }_2$, its spatial distribution seems rather large compared to the small-scale structures and this may be related to the large-scale flow. It is interesting to note that two different principal axes of stretching are associated in each case with different well-defined structures, but further discussions, determining these eigenvectors, are required. Figure 16 shows the distributions of the eigenvalues for higher $Wi$ and $Re$ and they again demonstrate clearly the $Wi$-dependence of polymer stretching. The degree of the polymer stretching depends on $Wi$, and it can be confirmed from the figure that local ${\sigma }_i$ values increase as $Wi$ increases. However, it should be noted that the trends observed in the low $Wi$ and $Re$ case in Figure 11 can be seen here as well. The small-scale spanwise-oriented structures are relevant to the largest eigenvalue (${\sigma }_3$), whereas the distribution of intense ${\sigma }_2$ resemblesthe large-scale flow. The minimum eigenvalue of ${\sigma }_1$ is approximately globally distributed with zero values, although there are locally large values, perhaps due to slightly-twisted spanwise vortices.

As discussed in relation to the visualizations and correlation coefficients, it was found that the EIT of plane Couette flow is a flow in which large- and small-scale flows coexist. The energy spectrum of the velocity fluctuation reveals the level of energy with each scale and the behavior of its decay. For a dilute polymer solution, the polymer has no effect on the dissipation of energy in the fluid, and the role of the polymer itself is to transfer the stored kinetic energy to a small-scale flow of the solvent [50]. Therefore, it is known that the EIT is different from the energy cascade of Newton turbulence. Dubief et al. [12] showed a $k^{-14/3}$ steep decay law, where k represents the wavenumber, in the energy spectrum of EIT, but there is no theoretical significance of the $-14/3$ exponent as yet. Additionally, Pereira et al. [51] observed a similar $k^{-14/3}$ decay law in high-drag-reduced turbulence. It is characteristic that decay scaling that is steeper than Newtonian turbulence appears in the relatively-high-wavenumber region. Furthermore, Pereira et al. [51,52] indicated that the spectrum has a steep decay law of $k^{-8}$, which is steeper than during the hibernation events appearing in high-drag-reduced turbulence.Although there is no theoretical support for this decay law, it is recognized as a common property of EIT that the energy spectrum has a steep decay law in the specific wavenumber region at the limit of inertia with high elasticity, that is, the MDR asymptote.

Figure 17 shows the time-averaged longitudinal one-dimensional energy spectra $E_{u_i{u}_j}$ as a function of the streamwise wavenumber $k_x$ for $Re=10$ and $Re=100$ with a varying Weissenberg number at the channel’s center plane. The energy spectra of velocity components in each direction and for each typical case of varying Weissenberg number are plotted. For $Re=10$, as confirmed in the discussion on correlation coefficients, a peak of the energy spectrum that appears in the low wavenumber region, i.e., a large-scale flow, is not formed in $E_{vv}$. The other components, $E_{uu}$ and $E_{ww}$, have similar profiles. The energy spectra of the streamwise and spanwise directions have the first peak in the low-wavenumber region, which corresponds to the large-scale flow, and have the second peak in the region where the wavenumber $k_x\delta \approx 10$, which corresponds to the spanwise-oriented structure. The decay of the energy spectrum in the streamwise direction $E_{uu}$ can be scaled by the $k_x^{-7}$ power law in the moderate-wavenumber domain $10<k_x\delta <30$. This suggests that there is an energy cascade from the scale of spanwise-oriented structures to the smaller scale. It can be seen that the decay of the energy spectrum follows a power law, scaled with the same exponent as for other components in this wavenumber region. As the peak of $E_{vv}$ is in the low-wavenumber region, a large-scale flow is formed for the higher Reynolds number $Re=100$. In well-developed EIT at a high Weissenberg number ($Wi=100$), the energy spectrum of the streamwise direction $E_{uu}$ for $Re=100$ can also be scaled with the same exponent as $E_{uu}$ of $Re=10$ in the same medium-wavenumber region. However, scaling with the same exponent deviates from the decay of the energy spectrum of the other components. Although there is no quantitative agreement between the decay law of the EIT energy spectrum confirmed in plane Couette flow and the decay law obtained in previous studies, based on this result, it can be inferred that EIT has a steeper decay law than that observed in Newtonian turbulence.

4. Conclusions

We performed a direct numerical simulation of viscoelastic-plane Couette flow at very low Reynolds numbers and investigated the characteristics of viscoelasticity-induced instability. The model used for simulating the viscoelastic fluid of dilute solutions was the Giesekus model and the specific parameters of the model that determined the characteristics of the fluid rheology were $\alpha =0.001$ and $\beta =0.9$. The Reynolds number, which is a flow parameter, was fixed at one of two very low values, $Re=10$ and $Re=100$. The Weissenberg number, which is a rheological parameter representing elastic strength, was gradually increased for each Reynolds number, focusing on the Weissenberg-number-dependence of viscoelasticity-induced instability. Instabilities were observed for both low and high Reynolds numbers, and we concluded that this instability represented EIT, based on its characteristic structural features. Although the present simulation for Re was limited to $Re=10$ and 100, the critical $Wi$ that triggered the transition to EIT was approximately $Wi\approx 20$. The structure of EIT had a cylindrical shape and was elongated in the spanwise direction and aligned in the streamwise direction. Moreover, the structure was concentrated in the region of the channel center. The topological structure identified by the second invariant of the velocity gradient tensor had alternating rotational and extensional/compression regions. The mean velocity profiles of the flow field in which this structure was formed had inflection points around the channel center. Except for the inflection points, the mean velocity profiles were almost unchanged, based on the solution of laminar flow, even if EIT occurred. The two-point correlation coefficients revealed that there were structures of different scales in the EIT. The small wavelengths with strong correlations corresponded to spanwise-oriented structures that characterized EIT, whereas the large wavelengths represented a large-scale flow elongated in the streamwise direction. The wavelength of the spanwise-oriented structures was approximately ${\lambda }_x\approx 0.5\delta$, and the wavelength of the large-scale flow was ${\lambda }_x=3$–$4\delta$ in the streamwise direction. The wavelength of spanwise-oriented structures depended on the Reynolds number, and smaller structures were observed at the higher Reynolds number. The spanwise-oriented structure lost its regularity as the Weissenberg number increased and became a disordered structure, whereas the large-scale flow developed and its intensity increased. The distributions of eigenvalues of $c_{ij}$, which represent the degree of polymer deformation, were similar to those of the structures that appeared in the flow field. One of the eigenvalues with a distribution similar to the small-scale structure had a much larger value than the other eigenvalues. The stretching of the polymer was inferred to be mainly due to spanwise-oriented structures. The magnitude of the eigenvalues increased as $Wi$ increased, and thus the degree of polymer elongation depended on $Wi$. The energy spectra of the velocity fluctuations indicated that there could be an energy cascade from the scale of spanwise-oriented structures to the smaller scale. The exponent of the decay law in the energy spectrum observed was larger than $-5/3$, which is well-known as the exponent of the energy cascade for Newtonian turbulence. We obtained a $k_x^{-7}$ power law for the streamwise energy spectrum $E_{uu}\left(k_x\right)$. This steep decay law may be one of the ubiquitous features of EIT.