The research on elasto-inertial turbulence (EIT), a new type of turbulent flow, has reached the stage of identifying the minimal flow unit (MFU). On this issue, direct numerical simulations of FENE-P fluid flow in two-dimensional channels with variable sizes are conducted in this study. We demonstrate with the increase of channel length that the simulated flow experiences several different flow patterns, and there exists an MFU for EIT to be self-sustained. At Weissenberg number ($Wi$) higher than the one required to excite EIT, when the channel length is relatively small, a steady arrowhead regime (SAR) flow structure and a laminar-like friction coefficient is achieved. However, as the channel length increases, the flow can fully develop into EIT characterized with high flow drag. Close to the size of the MFU, the simulated flow behaves intermittently between the SAR state with low drag and EIT state with high drag. The flow falling back to ‘laminar flow’ is caused by the insufficient channel size below the MFU. Furthermore, we give the relationship between the value of the MFU and the effective $Wi$, and explain its physical reasons. Moreover, the intermittent flow regime obtained based on the MFU gives us an opportunity to look into the origin and exciting process of EIT. Through capturing the onset process of EIT, we observed that EIT originates from the sheet-like extension structure located near the wall, which is maybe related to the wall mode rather than the centre mode. The fracture and regeneration of this sheet-like structure is the key mechanism for the self-sustaining of EIT.

Microorganism motility often takes place within complex, viscoelastic fluid environments, e.g. sperm in cervicovaginal mucus and bacteria in biofilms. In such complex fluids, strains and stresses generated by the microorganism are stored and relax across a spectrum of length and time scales and the complex fluid can be driven out of its linear response regime. Phenomena not possible in viscous media thereby arise from feedback between the swimmer and the complex fluid, making swimming efficiency co-dependent on the propulsion mechanism and fluid properties. Here, we parameterize a flagellar motor and filament properties together with elastic relaxation and nonlinear shear-thinning properties of the fluid in a computational immersed boundary model. We then explore swimming efficiency, defined as a particular flow rate divided by the torque required to spin the motor, over this parameter space. Our findings indicate that motor efficiency (measured by the volumetric flow rate) can be boosted or degraded by relatively moderate or strong shear thinning of the viscoelastic environment.

We investigate theoretically the steady incompressible viscoelastic flow in a rigid axisymmetric cylindrical pipe with varying cross-section. We use the Oldroyd-B viscoelastic constitutive equation to model the fluid viscoelasticity. First, we derive exact general formulae: for the total average pressure-drop as a function of the wall shear rate and the viscoelastic axial normal extra-stress; for the viscoelastic extra-stress tensor and the Trouton ratio as functions of the fluid velocity on the axis of symmetry; and for the viscoelastic extra-stress tensor along the wall in terms of the shear rate at the wall. Then we exploit the classic lubrication approximation, valid for small values of the square of the aspect ratio of the pipe, to simplify the original governing equations. The final equations are solved analytically using a regular perturbation scheme in terms of the Deborah number, De, up to eighth order in De. For a hyperbolically shaped pipe, we reveal that the reduced pressure-drop and the Trouton ratio can be recast in terms of a modified Deborah number, Dem, and the polymer viscosity ratio, η, only. Furthermore, we enhance the convergence and accuracy of the eighth-order solutions by deriving transformed analytical formulae using Padé diagonal approximants. The results show the decrease of the pressure drop and the enhancement of the Trouton ratio with increasing Dem and/or increasing η. Comparison of the transformed solutions with numerical simulations of the lubrication equations using pseudospectral methods shows excellent agreement between the results, even for high values of Dem and all values of η, revealing the robustness, validity and efficiency of the theoretical methods and techniques developed in this work. Last, it is shown that the exact solution for the Trouton ratio gives a well-defined and finite solution for any value of Dem and reveals the reason for the failure of the corresponding high-order perturbation series for Dem > 1/2.

We perform simulations of an impulsively started, axisymmetric viscoelastic jet exiting a nozzle and entering a stagnant gas phase using the open-source code Basilisk. This code allows for efficient computations through an adaptively refined volume-of-fluid technique that can accurately capture the deformation of the liquid–gas interface. We use the FENE-P constitutive equation to describe the viscoelasticity of the liquid, and employ the log-conformation transformation, which provides stable solutions for the evolution of the conformation tensor as the jet thins down under the action of interfacial tension. For the first time, the entire jetting and breakup process of a viscoelastic fluid is simulated, including the pre-shearing flow through the nozzle, which results in an inhomogeneous initial radial stress distribution in the fluid thread that affects the subsequent breakup dynamics. The evolution of the velocity field and the elastic stresses in the nozzle are validated against analytical solutions where possible, and the early-stage dynamics of the jet evolution are compared favourably to the predictions of linear stability theory. We study the effect of the flow inside the nozzle on the thinning dynamics of the viscoelastic jet (which develops distinctive ‘beads-on-a-string’ structures) and on the spatio-temporal evolution of the polymeric stresses in order to systematically explore the dependence of the filament thinning and breakup characteristics on the initial axial momentum of the jet and the extensibility of the dissolved polymer chains.

The peculiar migration and rotational dynamics of non-spherical particles in non-Newtonian flows stem from the interplay between fluid rheology and fluid inertial effects. In this paper, the cross-flow migration of a neutrally buoyant oblate spheroid (aspect ratio $AR = 0.5$) immersed in the elasto-inertial duct flow is investigated by particle-resolved simulations with the immersed boundary method. Different from spherical particles, due to the orientation-dependent lift force, the oblate spheroid migrates in an oscillating manner in the duct. The travelling period for particles reaching the duct centreline undergoes a non-monotonic change with elastic number, revealing the existence of a critical elastic number governing the migrating efficiency of oblate particles within the present flow system. For the particle rotation and orientation, the present results indicate that the particle rotation rate and orbit drifting rate are both hindered by the fluid elasticity. With increasing the fluid elasticity, three different orientation modes – log-rolling mode, kayaking-like mode and longside-flow alignment mode – are observed successively during the elasto-inertial migration of the oblate spheroid. Potentially, the present results could be used to design the rheology-based controlling strategy for guiding particles to achieve optimal focusing and orientation in microfluidic applications without the need for external forcing fields.

The elasto-inertial focusing and rotating characteristics of spheroids in a square channel flow of Oldroyd-B viscoelastic fluids are studied by the direct forcing/fictitious domain method. The rotational behaviours, changes in the equilibrium positions and travel distances are explored to analyse the mechanisms of spheroid migration in viscoelastic fluids. Within the present simulated parameters (1 ≤ Re ≤ 100, 0 ≤ Wi ≤ 2, 0.4 ≤ α ≤3), the results show that there are four kinds of equilibrium positions and six (five) kinds of rotational behaviours for the elasto-inertial migration of prolate (oblate) spheroids. We are the first to identify a new rotational mode for the migration of prolate spheroids. Only when the particles are initially located at a corner and wall bisector, some special initial orientations of the spheroids have an impact on the final equilibrium position and rotational mode. In other general initial positions, the initial orientation of the spheroid has a negligible effect. A higher Weissenberg number means the faster the particles migrate to the equilibrium position. The spheroid gradually changes from the corner (CO), channel centreline (CC), diagonal line (DL) and cross-section midline (CSM) equilibrium positions as the elastic number decreases, depending on the aspect ratio, initial orientation and rotational behaviour of the particles and the elastic number of the fluid. When the elastic number is less than the critical value, the types of rotational modes of the spheroids are reduced. By controlling the elastic number near the critical value, spheroids with different aspect ratios can be efficiently separated.

Many natural and industrial processes involve the flow of fluids made of solid particles suspended in non-Newtonian liquid matrices, which are challenging to control due to the fluid's nonlinear rheology. In the present work, a Taylor–Couette canonical system is used to investigate the flow of dilute to semi-dilute suspensions of neutrally buoyant spherical particles in highly elastic base polymer solutions. Friction measurement synchronized with direct flow visualization are combined to characterize the critical conditions for the onset of elasto-inertial instabilities, expected here as a direct transition to elasto-inertial turbulence (EIT). Adding a low particle volume fraction (${\leq }2\,\%$, dilute regime) does not affect the nature of the primary transition and reduces the critical Weissenberg number for the onset of EIT, despite a significant decrease in the apparent fluid elasticity. However, for particle volume fractions ${\geq }6\,\%$ (semi-dilute regime), EIT is no longer observed in the explored Reynolds range, suggesting an apparent relaminarization with yet not further decrease in fluid elasticity. Instead, a new regime, termed here elasto-inertial dissipative (EID), was uncovered. It originates from particle–particle interactions altering particle–polymer interactions and occurring under elasto-inertial conditions comparable to those of EIT. Increasing particle volume fraction in the semi-dilute regime and, in so, the particle contribution to the overall viscosity, delays the onset of EID similarly to what was observed previously for EIT in lower elasticity fluids. After this onset, a decrease in the pseudo-Nusselt number observed with increasing inertia and particle-to-polymer concentration ratio confirms a particle-induced alteration of energy transfer in the flow.

Cross-stream migration of a deformable fluid particle is investigated computationally in a pressure-driven channel flow of a viscoelastic fluid via interface-resolved simulations. Flow equations are solved fully coupled with the Giesekus model equations using an Eulerian–Lagrangian method and extensive simulations are performed for a wide range of flow parameters to reveal the effects of particle deformability, fluid elasticity, shear thinning and fluid inertia on the particle migration dynamics. Migration rate of a deformable particle is found to be much higher than that of a solid particle under similar flow conditions mainly due to the free-slip condition on its surface. It is observed that the direction of particle migration can be altered by varying shear thinning of the ambient fluid. With a strong shear thinning, the particle migrates towards the wall while it migrates towards the channel centre in a purely elastic fluid without shear thinning. An onset of elastic flow instability is observed beyond a critical Weissenberg number, which in turn causes a path instability even for a nearly spherical particle. An inertial path instability is also observed once particle deformation exceeds a critical value. Shear thinning is found to be suppressing the path instability in a viscoelastic fluid with a high polymer concentration whereas it reverses its role and promotes path instability in a dilute polymer solution. It is found that migration of a deformable particle towards the wall induces a secondary flow with a velocity that is approximately an order of magnitude higher than the one induced by a solid particle under similar flow conditions.

Elastoinertial turbulence (EIT) is a chaotic flow resulting from the interplay between inertia and viscoelasticity in wall-bounded shear flows. Understanding EIT is important because it is thought to set a limit on the effectiveness of turbulent drag reduction in polymer solutions. Here, we analyse simulations of two-dimensional EIT in channel flow using spectral proper orthogonal decomposition (SPOD), discovering a family of travelling wave structures that capture the sheetlike stress fluctuations that characterise EIT. The frequency-dependence of the leading SPOD mode contains distinct peaks and the mode structures corresponding to these peaks exhibit well-defined travelling structures. The structure of the dominant travelling mode exhibits shift–reflect symmetry similar to the viscoelasticity-modified Tollmien–Schlichting (TS) wave, where the velocity fluctuation in the travelling mode is characterised by large-scale regular structures spanning the channel and the polymer stress field is characterised by thin, inclined sheets of high polymer stress localised at the critical layers near the channel walls. The travelling structures corresponding to the higher-frequency modes have a very similar structure, but are nested in a region roughly bounded by the critical layer positions of the next-lower-frequency mode. A simple theory based on the idea that the critical layers of mode $\kappa$ form the ‘walls’ for the structure of mode $\kappa +1$ yields quantitative agreement with the observed wave speeds and critical layer positions, indicating self-similarity between the structures. The physical idea behind this theory is that the sheetlike localised stress fluctuations in the critical layer prevent velocity fluctuations from penetrating them.

Motivated by the recent numerical results of Khalid et al. (Phys. Rev. Lett., vol. 127, 2021, 134502), we consider the large-Weissenberg-number ($W$) asymptotics of the centre mode instability in inertialess viscoelastic channel flow. The instability is of the critical layer type in the distinguished ultra-dilute limit where $W(1-\beta )=O(1)$ as $W \rightarrow \infty$ ($\beta$ is the ratio of solvent-to-total viscosity). In contrast to centre modes in the Orr–Sommerfeld equation, $1-c=O(1)$ as $W \rightarrow \infty$, where $c$ is the phase speed normalised by the centreline speed as a central ‘outer’ region is always needed to adjust the non-zero cross-stream velocity at the critical layer down to zero at the centreline. The critical layer acts as a pair of intense ‘bellows’ which blows the flow streamlines apart locally and then sucks them back together again. This compression/rarefaction amplifies the streamwise-normal polymer stress which in turn drives the streamwise flow through local polymer stresses at the critical layer. The streamwise flow energises the cross-stream flow via continuity which in turn intensifies the critical layer to close the cycle. We also treat the large-Reynolds-number ($Re$) asymptotic structure of the upper (where $1-c=O(Re^{-2/3})$) and lower branches of the $Re$–$W$ neutral curve, confirming the inferred scalings from previous numerical computations. Finally, we remark that the viscoelastic centre-mode instability was actually first observed in viscoelastic Kolmogorov flow by Boffetta et al. (J. Fluid Mech., vol. 523, 2005, pp. 161–170).

We study the spreading of Newtonian viscous (aqueous glycerin solution) and viscoelastic (aqueous polymer solution) drops on solid substrates with different wettabilities. For drops of the same zero-shear viscosity, we find in the early stages of spreading that viscoelastic drops (i) spread faster and (ii) their contact radius shows a different power law vs time than Newtonian drops. We argue that the effect of viscoelasticity is only observable for experimental time scales of the order of or larger than the internal relaxation time of the viscoelastic polymer solution. We attribute this behaviour to the shear thinning of the viscoelastic polymer solution. When approaching the contact line, the shear rate increases and the steady-state viscosity of the viscoelastic drop is lower than that of the Newtonian drop. We support our experimental findings with a simple (first-order) perturbation model that qualitatively agrees with our findings.

Pressure-driven flows of viscoelastic fluids in narrow non-uniform geometries are common in physiological flows and various industrial applications. For such flows, one of the main interests is understanding the relationship between the flow rate $q$ and the pressure drop $\Delta p$, which, to date, is studied primarily using numerical simulations. We analyse the flow of the Oldroyd-B fluid in slowly varying arbitrarily shaped, contracting channels and present a theoretical framework for calculating the $q-\Delta p$ relation. We apply lubrication theory and consider the ultra-dilute limit, in which the velocity profile remains parabolic and Newtonian, resulting in a one-way coupling between the velocity and polymer conformation tensor. This one-way coupling enables us to derive closed-form expressions for the conformation tensor and the flow rate–pressure drop relation for arbitrary values of the Deborah number ($De$). Furthermore, we provide analytical expressions for the conformation tensor and the $q-\Delta p$ relation in the high-Deborah-number limit, complementing our previous low-Deborah-number lubrication analysis. We reveal that the pressure drop in the contraction monotonically decreases with $De$, having linear scaling at high Deborah numbers, and identify the physical mechanisms governing the pressure drop reduction. We further elucidate the spatial relaxation of elastic stresses and pressure gradient in the exit channel following the contraction and show that the downstream distance required for such relaxation scales linearly with $De$.

The long- and finite-wavelength instabilities of weakly viscoelastic film on an oscillating inclined plane are investigated. By using the Chebyshev series solution with the Floquet theory, the combined effects of viscoelasticity and forcing amplitude on instability are described when the inclined plane oscillates in streamwise and wall-normal directions. For long-wavelength instability, the solution to the eigenvalue problem is obtained analytically by the asymptotic expansion method. Results show that the Floquet exponent is independent of the wall-normal oscillation amplitude. The effects of inclined angle, gravity and surface tension on the stability of viscoelastic liquid film are also discussed. For finite-wavelength instability, numerical results corresponding to the wall-normal oscillation disclose that with the increase of viscoelasticity, the unstable gravitational boundary moves to a higher wavenumber, and the critical amplitudes of subharmonic and harmonic instabilities are reduced. The neutral curve of gravity instability for streamwise oscillatory flow is divided into two parts, and a stable bandwidth is formed for a large value of the viscoelastic parameter. Besides, a new oscillatory mode is identified at small angles of inclination, which will be enhanced if the effect of viscoelasticity is incorporated.

The two-dimensional, steady, homogeneous flow field proposed by Astarita (J. Rheol., vol. 35, 1991, pp. 687–689) is studied for a range of viscoelastic constitutive equations of differential form including the models due to Oldroyd (the upper and lower convected Maxwell; UCM/LCM), Phan-Thien and Tanner (simplified, linear form; sPTT) and Giesekus. As the flow is steady and homogeneous, the sPTT model results also give the FENE-P model solutions via a simple transformation of parameters. The flow field has the interesting feature that a scalar parameter may be used to vary the flow ‘type’ continuously from solid-body rotation to simple shearing to planar extension whilst the rate of deformation tensor $\boldsymbol{\mathsf{D}}$ remains constant (i.e. independent of flow type). The response of the models is probed in order to determine how a scalar ‘viscosity’ function may be rigorously constructed which includes flow-type dependence. We show that for most of these models – the Giesekus being the exception – the first and second invariants of the resulting extra stress tensor are linearly related, and for models based on the upper convected derivative, this link is simply via a material property, i.e. half the modulus. By defining a frame-invariant coordinate system with respect to the eigenvectors of $\boldsymbol{\mathsf{D}}$, we associate a ‘viscosity’ for each of the flows to a deviatoric stress component and show how this quantity varies with the flow-type parameter. For elliptical motions, rate thinning is always observed and all models give essentially the UCM response. For strong flows, i.e. flow types containing at least some extension, thickening occurs and only a small element of extension is required to remove any shear thinning inherent in the model (e.g. as occurs in steady simple shearing for the sPTT/Giesekus models). Finally, a functional form of a viscosity equation which could incorporate flow type, but be otherwise inelastic, the so-called GNFFTy (generalised Newtonian fluid model incorporating flow type, pronounced ‘nifty’), is proposed. In the frame-invariant coordinate system proposed, this model is also capable of capturing normal-stress differences.

We combine flow visualisation techniques and particle image velocimetry to experimentally investigate the higher-order transition to elastoinertial turbulence of Boger fluids ($El = 0.11\unicode{x2013}0.34$) in Taylor–Couette flows. The observed route to turbulence is associated with the appearance of chaotic inflow jets, termed flame patterns, for increasing inertia, and stable structures of solitons, known as diwhirls, for decreasing inertia. We also report for the first time spatially and temporally resolved flow fields in the meridional plane for the three characteristic viscoelastic flow regimes (diwhirls, flame patterns and elastoinertial turbulence). We observe in all cases coherent structures of dynamically independent solitary vortex pairs. The stability of these coherent structures is jet-dominated and can be mainly ascribed to the high extension of the polymer chains in the inflow boundaries in the $r$–$z$ plane. Solitary pairs are self-sustained when created through random events and do not split; instead, they merge when moving sufficiently close and annihilate when hoop stresses are not sufficient to sustain them. The highly localised and random events result in highly fluctuating, chaotic flow states. We estimate the decay exponent of spatial power spectral density, illustrating a universal scaling of $-2.5$ for elastoinertial turbulence. Based on our observations and in an effort to unify and combine precedent theories with our results, we suggest a mechanism for the origins of elastoinertial instabilities, accounting for both the effect of elasticity on the vortex formation and the effect of increasing/decreasing inertia on the flow dynamics.

It is well known that inertia-free shearing flows of a viscoelastic fluid with curved streamlines, such as the torsional flow between a rotating cone and plate or the flow in a Taylor–Couette geometry, can become unstable to a three-dimensional time-dependent instability at conditions exceeding a critical Weissenberg ($Wi$) number. However, the combined effects of fluid elasticity, shear thinning and finite inertia (as quantified by the Reynolds number $Re$) on the onset of elasto-inertial instabilities are not fully understood. Using a set of cone–plate geometries, we experimentally explore the entire $Wi$–$Re$ phase space for a series of nonlinear viscoelastic fluids (with the dependence on shear rate $\dot{\gamma}$ quantified using a shear-thinning parameter $\beta _P(\dot {\gamma })$). We tune $\beta _P(\dot {\gamma })$ by varying the dissolved polymer concentration in solution. This progressively reduces shear thinning but leads to finite inertial effects before the onset of elastic instability, and thus naturally results in elasto-inertial coupling. Time-resolved rheometric measurements and flow visualization experiments allow us to investigate the effects of flow geometry, and document the combined effects of varying $Wi, Re$ and $\beta _P(\dot {\gamma })$ on the emergence of secondary motions at the onset of instability. The resulting critical state diagram quantitatively depicts the competition between the stabilizing effects of shear thinning and the destabilizing effects of inertia. We extend the curved streamline instability criterion of Pakdel & McKinley (Phys. Rev. Lett., vol. 77, no. 12, 1996, p. 2459) for the onset of purely elastic instability in curvilinear geometries by using scaling arguments to incorporate shear thinning and finite inertial effects. The augmented condition facilitates predictions of the onset of instability over a broader range of flow conditions, and thus bridges the gap between purely elastic and elasto-inertial curved streamline instabilities.

We study the motion and deformation of a single bubble rising inside a cylindrical container filled with a viscoelastic material. We solve numerically the mass and momentum balances along with the constitutive equation for the viscoelastic stresses, without assuming axial symmetry to allow the growth of three-dimensional disturbances. Hence, we may predict the emergence of the notorious knife-edge shape of the bubble, which is a result of a purely elastic instability triggered in the locality of the trailing edge. Our results compare well with existing experiments. We visualize, for the first time to the best of our knowledge, the flow kinematics and dynamics that arise downstream of the bubble. We propose two quantities, one kinematical and one geometrical, for the determination of the onset of the instability. We demonstrate that extension-rate thinning in the constitutive law is necessary for the emergence of the instability. Moreover, our results indicate that increasing (a) the deformability of the bubble and (b) the initial extension rate hardening of the viscoelastic material, prior to thinning, triggers the instability earlier. These novel findings help us formulate and propose a mechanism that controls the onset of the instability and explain why the knife-edge shape is not encountered as frequently.

The classical Cox–Voinov theory of contact line motion provides a relation between the macroscopically observable contact angle, and the microscopic wetting angle as a function of contact-line velocity. Here, we investigate how viscoelasticity, specifically the normal stress effect, modifies the wetting dynamics. Using the thin film equation for the second-order fluid, it is found that the normal stress effect is dominant at small scales yet can significantly affect macroscopic motion. We show that the effect can be incorporated in the Cox–Voinov theory through an apparent microscopic angle, which differs from the true microscopic angle. The theory is applied to the classical problems of drop spreading and dip coating, which shows how normal stress facilitates (inhibits) the motion of advancing (receding) contact lines. For rapid advancing motion, the apparent microscopic angle can tend to zero, in which case the dynamics is described by a regime that was already anticipated in Boudaoud (Eur. Phys. J. E, vol. 22, 2007, pp. 107–109).

This study explores the dynamics of finite-size fibres suspended freely in a viscoelastic turbulent flow. For a fibre suspended in Newtonian flows, two different flapping regimes were identified by Rosti et al. (Phys. Rev. Lett., vol. 121, issue 4, 2018, 044501): one dominated by time scales from the flow, and another dominated by time scales associated with its natural frequency. We explore in this work how the fibre dynamics is modified by the elasticity of the carrier fluid. For this, we perform direct numerical simulations of a two-way coupled fibre–fluid system in a parametric space spanning different Deborah numbers, fibre bending stiffness (flexible to rigid) and linear density difference between the fibre and the flow (neutrally buoyant to denser-than-fluid fibres). We examine how these parameters influence various fibre characteristics such as the frequency of flapping, curvature, and alignment with the fluid strain and polymer stretching directions. Results reveal that the neutrally buoyant fibres, depending on their flexibility, oscillate with large and small time scales transpiring from the flow, but the smaller time scales are suppressed as the polymer elasticity increases. Polymer stretching is uncommunicative to denser-than-fluid fibres, which flap with large time scales from the flow when flexible, and with their natural frequency when rigid. Thus the characteristic elastic time scale has a subdominant effect when the fibres are neutrally buoyant, while its effect is absent when the fibres become more inertial. In addition, we also explore the fibre's bending curvature and its preferential alignment with the flow to identify the other roles of viscoelasticity in modifying the coupled fluid–structure dynamics. Inertial fibres have larger curvatures and are less responsive to the polymer presence, whereas the neutrally buoyant fibres show quantitative changes. The perceptible passivity of the denser fibres is again reflected in the way they align preferentially with the polymeric stretching directions: the neutrally buoyant fibres show a higher alignment with the polymer stretching directions compared to the denser ones. In a nutshell, the polymers exert a larger influence on neutrally buoyant fibres, which are more reflective of the polymeric influence in the flow. The study addresses comprehensively the interplay between polymer elasticity and the fibre structural properties in determining its response behaviour in an elasto-inertial turbulent flow.

The viscoelasticity of a dilute bubble suspension is theoretically derived from the constitutive equation originally for a dilute emulsion proposed by Frankel & Acrivos (J. Fluid Mech., vol. 44, issue 1, 1970, pp. 65–78). Non-dimensionalization of the original tensor equation indicates that the viscoelasticity is systematized for a given void fraction by the capillary number $Ca$ and dynamic capillary number $Cd$, representing the bubble deformability and unsteadiness of bubble deformation. Comprehensive evaluation of the viscoelasticity according to the volume fraction, $Ca$ and $Cd$ reveals that whether the viscosity increases or decreases depends on whether $Ca$ or $Cd$ exceeds a common critical value. In addition, it is indicated that the bubble suspension has the most prominent viscoelasticity when the time scale of the shear deformation is the same as the relaxation time of the suspended bubble and when the bubbles keep a spherical shape, that is, $Ca \ll 1$ and $Cd = 1$. The applicability of this theory in flow prediction was examined in a Taylor–Couette system, and experimentally good agreement was confirmed.