1. Introduction
Active matter represents a range of non-equilibrium systems comprising self-propelled units such as active particles (Datt & Elfring Reference Datt and Elfring2019; Maity & Burada Reference Maity and Burada2022), bacteria, cells (Dombrowski et al. Reference Dombrowski, Cisneros, Chatkaew, Goldstein and Kessler2004; Sokolov et al. Reference Sokolov, Aranson, Kessler and Goldstein2007; Dunkel et al. Reference Dunkel, Heidenreich, Drescher, Wensink, Bär and Goldstein2013) or biopolymers (Fürthauer et al. Reference Fürthauer, Lemma, Foster, Ems-McClung, Yu, Walczak, Dogic, Needleman and Shelley2019; Lemma et al. Reference Lemma, Norton, Tayar, DeCamp, Aghvami, Fraden, Hagan and Dogic2021; Chandrakar et al. Reference Chandrakar2022; Zarei et al. Reference Zarei, Berezney, Hensley, Lemma, Senbil, Dogic and Fraden2023). Active fluids are a typical class of active matter which can form a gas or liquid phase or be suspended in another liquid. The constituents in an active fluid can convert other forms of energy into mechanical work by inducing stresses and interacting with the environment, such as driving flows in the background fluid, moving against neighbour units and undergoing phase separation. Active fluids are promising for emergent applications in biomedicine (Ghosh et al. Reference Ghosh, Xu, Gupta and Gracias2020), robotics (Ceron et al. Reference Ceron, Kimmel, Nilles and Petersen2021) and materials science (Needleman & Dogic Reference Needleman and Dogic2017), and have the potential to revolutionise our understanding of how living systems operate and to provide new design principles for biomimetic, functional and autonomous materials (Zhang et al. Reference Zhang, Mozaffari and de Pablo2021). However, our current understanding of their dynamic behaviours is still overwhelmed by the microscopic details of the specific system.
A paradigmatic class of active fluids is active liquid crystals (LCs) (Zhang et al. Reference Zhang, Mozaffari and de Pablo2021), which comprise locally aligned dense units with anisotropic shapes, forming different LC phases in non-equilibrium conditions. One of such active LC phases is namely active nematic (Doostmohammadi et al. Reference Doostmohammadi, Ignés-Mullol, Yeomans and Sagués2018), which widely exists in many biological and synthetic systems (Narayan, Ramaswamy & Menon Reference Narayan, Ramaswamy and Menon2007; Sanchez et al. Reference Sanchez, Chen, DeCamp, Heymann and Dogic2012; Zhou et al. Reference Zhou, Sokolov, Lavrentovich and Igor2014; Kawaguchi, Kageyama & Sano Reference Kawaguchi, Kageyama and Sano2017; Saw et al. Reference Saw, Doostmohammadi, Nier, Kocgozlu, Thampi, Toyama, Marcq, Lim, Yeomans and Ladoux2017; Kumar et al. Reference Kumar, Zhang, De, Juan and Gardel2018; Li et al. Reference Li, Shi, Huang, Chen, Xiao, Liu, Chaté and Zhang2019). Though active LC systems are diverse, their spontaneous flows can be well understood by a hydrodynamic model (Denniston, Orlandini & Yeomans Reference Denniston, Orlandini and Yeomans2001; Marenduzzo et al. Reference Marenduzzo, Orlandini, Cates and Yeomans2007a ), in which active stresses drive positive topological defects into self-propulsion (Aditi Simha & Ramaswamy Reference Aditi Simha and Ramaswamy2002; Shankar et al. Reference Shankar, Ramaswamy, Marchetti and Bowick2018). Nevertheless, the success of the hydrodynamic model of active LCs cannot be extended directly to active fluids without LC order (Marchetti et al. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Simha2013), in which a nematic director field does not exist.
In the present work, we propose an active viscoelastic model in an attempt to provide a unified framework to describe distinct isotropic active fluids. Here we focus on two experimental systems. One is microtubule-based active fluids, in which the active motion of the kinesin motor clusters can lead to the relative slidings of microtubule bundles, driving different spontaneous flow patterns tuneable by confinement geometry (Wu et al. Reference Wu, Hishamunda, Chen, DeCamp, Chang, Fernández-Nieves, Fraden and Dogic2017). The other active fluid is bacterial suspensions (Liu et al. Reference Liu, Shankar, Marchetti and Wu2021; Nishiguchi et al. Reference Nishiguchi, Shiratani, Takeuchi and Aranson2024). In addition to geometric confinements (Wioland, Lushi & Goldstein Reference Wioland, Lushi and Goldstein2016; Bhattacharjee & Datta Reference Bhattacharjee and Datta2019), the inclusion of different amounts of DNA molecules in the bacteria-based active fluid can also be used to control its flow state via tuning its viscoelasticity (Liu et al. Reference Liu, Shankar, Marchetti and Wu2021). As DNA concentration increases, the spontaneous flow in a disk can transition from a turbulent-like state to a unidirectional giant vortex, and to an oscillatory mode in which a single vortex flow can periodically switch between clockwise (CW) and counter-clockwise (CCW) directions (Liu et al. Reference Liu, Shankar, Marchetti and Wu2021). Despite the richness of the flow patterns observed in these different isotropic active fluidic systems, little is known about the universality of these flow patterns and how they depend on the system specifics.
Figure 1. Flow state transitions of an active viscoelastic fluid in a 2-D channel. (a) Flow state diagram in terms of the activity
$\zeta$
and Weissenberg number
$\textrm {Wi}$
. The background colour represents the magnitude of the time-averaged flow rate
$Q$
. The solid magenta line marks
$\zeta _c$
obtained by linear stability analysis. (b) Schematic of the network structure for the PTT model. (c) Unidirectional flow state (
$\textrm {Wi}=10^{-0.5}$
,
$\zeta =3.35$
). The left half-part of the channel shows the velocity field
$\textbf {U}$
with background colour indicating the magnitude
$|\textbf {U}|$
. The other half-part shows the orientation
$\textbf {r}$
of the polymer molecules and their order magnitude
$S$
. (d) Travelling-wave state (
$\textrm {Wi}=10^{-0.5}$
,
$\zeta =3.4$
). (e) Vortex-roll state (
$\textrm {Wi}=10^{-0.5}$
,
$\zeta =3.5$
).
$\pm 1/2$
defects are present. (f) Dancing state (
$\textrm {Wi}=10^{-0.5}$
,
$\zeta =3.8$
). (g) Turbulent-like state (
$\textrm {Wi}=10^{-0.5}$
,
$\zeta =4.2$
).
The hydrodynamic model of active LCs has been used to understand isotropic active fluids in both two-dimensional (2-D) (Fielding, Marenduzzo & Cates Reference Fielding, Marenduzzo and Cates2011; Caballero, You & Marchetti Reference Caballero, You and Marchetti2023) and three-dimensional (3-D) geometries (Chandragiri et al. Reference Chandragiri, Doostmohammadi, Yeomans and Thampi2020; Varghese et al. Reference Varghese, Baskaran, Hagan and Baskaran2020). In these models, a nematic order parameter is required. For many polymer-based viscoelastic liquids, however, there is no LC ordering but instead require a stress tensor.
In this study, we try to ask whether a hydrodynamic model of active fluids without LC order can be used to describe various isotropic active fluidic phenomena. To this end, we modify a viscoelastic liquid model by including an active stress term (Aditi Simha & Ramaswamy Reference Aditi Simha and Ramaswamy2002). We then use the model to study different 2-D confinement geometries. Our model produces rich spontaneous flow states (see figure 1 a) and can reproduce distinct dynamical flow patterns observed in the various experiments, which we elaborate in the following.
2. Problem formulation
2.1. Governing equations
We simulate a 2-D active viscoelastic fluid using the Phan-Thien–Tanner (PTT) model (Thien & Tanner Reference Thien and Tanner1977; Phan-Thien Reference Phan-Thien1978) in which the constitutive equation is derived from a network theory for polymeric liquids (see figure 1 b). The network junctions in this model are allowed to exhibit a certain degree of effective slip relative to the background continuum, and the relaxation time is assumed to be dependent on the local stress. The dimensionless governing equations of the active viscoelastic system (derivation presented in the supplementary material available at https://doi.org/10.1017/jfm.2025.177) consist of the incompressible Navier–Stokes equation and the constitutive equation:
\begin{gather} \nabla \cdot \mathbf {U}=0, \end{gather}
\begin{gather} \frac {\partial \mathbf {U}}{\partial t}=-\nabla P +\nabla ^2\mathbf {U}+\frac {1-\beta }{\beta }(1-\zeta ) \nabla \cdot \boldsymbol {\tau } - \nu \mathbf {U}, \end{gather}
\begin{gather} \frac {\partial \boldsymbol {\tau }}{\partial t} + \mathbf {U}\cdot \nabla \boldsymbol {\tau } + \xi \left (\boldsymbol {\tau }\cdot \mathbf {D} +\mathbf {D}\cdot \boldsymbol {\tau }\right ) - \boldsymbol {\tau }\cdot \nabla \mathbf {U} - \left [\boldsymbol {\tau }\cdot \left (\nabla \mathbf {U}\right ) \right ]^T =\frac {2}{\mathrm {\mathrm {Wi}}}\mathbf {D} - \frac {f(\tau _{kk})}{\mathrm {Wi}}\boldsymbol {\tau } + \kappa \nabla ^2\boldsymbol {\tau }, \end{gather}
where
$P$
denotes pressure, and
$\mathbf {U}$
and
$\boldsymbol {\tau }$
are the velocity vector and polymer-induced symmetric stress tensor, respectively. Here
$\mathbf {D}=[\nabla \mathbf {U}+(\nabla \mathbf {U})^T]/2$
is the rate-of-strain tensor. The two dimensionless control parameters are the Weissenberg number
$\mathrm {Wi}=U_c\lambda /L_c$
with
$U_c=\eta _s/(\rho _0L_c)$
, where
$\nu$
is the friction coefficient and
$\rho_0$
denotes the liquid density, and the viscosity ratio
$\beta =\eta _s/(\eta _s+\eta _p)$
, where
$\lambda$
is the polymer relaxation time for
$\tau _{kk}\to 0$
,
$L_c$
represents a characteristic length of the system,
$\tau_{kk}$
refers to the trace of the stress tensor, and
$\eta _s$
and
$\eta _p$
are the contributions to viscosity from the solvent and the polymer, respectively. To account for activity of the self-propelled units, we adopt an active stress term
$-\zeta \boldsymbol {\tau }$
that is linear with respect to the stress tensor, which implies that any anisotropy in the conformation of the polymer can induce an anisotropic stress that enhances such anisotropy; this term is similar to the swim pressure proposed by Omar, Wang & Brady (Reference Omar, Wang and Brady2020) and the active stress term in the active LC theory (Marenduzzo et al. Reference Marenduzzo, Orlandini and Yeomans2007b
). The effect of substrate friction from the third dimension is described by
$-\nu \mathbf {U}$
in equation (2.1
b). We consider Stokes limit by dropping the non-linear inertial term in equation (2.1
b). An artificial diffusion term
$\kappa \nabla ^2\boldsymbol {\tau }$
is introduced in equation (2.1
c) to enhance numerical stability (here
$\kappa =10^{-3}$
). In the PTT model, the dimensionless parameter
$\xi$
(
$0\lt \xi \lt 1$
) accounts for the slip between the network and the background continuous medium. Additionally, the linear relaxation function
$f(\tau _{kk})=1+\epsilon \mathrm {Wi} \tau _{kk}$
(Thien & Tanner Reference Thien and Tanner1977) is used, where the dimensionless parameter
$\epsilon$
controls the extensibility of the polymer (Thien & Tanner Reference Thien and Tanner1977) which is different from the constant molecular length in the active LC model. In addition, compared with the pure active LC model (Marenduzzo et al. Reference Marenduzzo, Orlandini and Yeomans2007b
) and the combined model for active LCs and viscoelastic polymers (Hemingway et al. Reference Hemingway, Maitra, Banerjee, Marchetti, Ramaswamy, Fielding and Cates2015), the governing equations (2.1) have different forms, particularly in the momentum equation.
2.2. Simulation details
The governing equations (2.1) are numerically solved using the open-source pseudo-spectral package Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). The Fourier basis is applied for the periodic directions and the Chebyshev polynomials for the wall-normal direction. The time integration is performed using the RK222 scheme which is a second-order, two-step, implicit/explicit Runge–Kutta method. On the wall, a no-slip boundary condition is used for the velocity field, and a Neumann condition
$\partial \boldsymbol {\tau }/\partial w_{\perp }=0$
is used for the polymer stress, with
$w_{\perp }$
denoting the wall-normal direction. In contrast to LC models, here no preferred molecular alignment direction is imposed on the boundary. All the simulations reported here are started with zero velocity and a randomised polymer stress of small magnitude. The flow rate
$Q$
is defined by
$Q=\iint U_{\parallel } \textrm {d}\textbf {x}/\mathcal {V}$
, where
$U_{\parallel }$
denotes the velocity component parallel to the wall, and
$\mathcal {V}$
denotes the volume of the system. Here
$Q$
is used as an order parameter to distinguish the coherent and turbulent flow states. When
$Q\neq 0$
, a coherent flow emerges. In the following, the viscosity ratio is set to
$\beta =0.7$
which is close to the experimental observations (Wu et al. Reference Wu, Hishamunda, Chen, DeCamp, Chang, Fernández-Nieves, Fraden and Dogic2017; Gagnon et al. Reference Gagnon, Dessi, Berezney, Boros, Chen, Dogic and Blair2020). The active viscoelastic fluid is assumed to exhibit a strong slip effect and weak extensibility, with parameters
$\xi =0.7$
and
$\epsilon =10^{-2}$
. The friction force is not considered (i.e.
$\nu =0$
) unless specified otherwise.
3. Results
3.1. Channel geometry
We start with a flat channel with lateral length
$L=4$
and grid resolution
$N_x\times N_y= 144\times 36$
. Above a critical activity
$\zeta _c$
, there is an onset of a spontaneous flow (i.e. unidirectional flow) as shown in figure 1(a). Weissenberg number is defined as
$\textrm {Wi}=\eta _s \lambda / \{{\rho _0 L_c^2}\}$
(with
$L_c=H^*$
in figure 1(b), the dimensionless channel height would be
$1$
). As
$\textrm {Wi}$
increases,
$\zeta _c$
becomes higher, which is consistent with our linear stability analysis (the solid magenta line in figure 1
a) as well as the mixed model of LC and viscoelastic liquids (Hemingway et al. Reference Hemingway, Maitra, Banerjee, Marchetti, Ramaswamy, Fielding and Cates2015). Note that
$\zeta _c$
approaches a constant value at a vanishing
$\textrm {Wi}$
. This limit can be understood by considering an unconfined 2-D domain. Indeed, our stability analysis of an infinite 2-D domain shows that
$\zeta _c$
is only determined by the viscosity ratio, namely
$\zeta _c=1/(1-\beta )$
(see supplementary material). For
$\beta =0.7$
used in the simulation, the critical activity
$\zeta _c=1/(1-\beta )\approx 3.3$
agrees quantitatively well with the vanishing
$\textrm {Wi}$
result in figure 1(a).
As activity further increases, the steady unidirectional flow sequentially transitions into travelling-wave, vortex-roll, dancing and finally turbulent state at the highest
$\zeta$
as shown in figure 1(c–g) (see supplementary movie 1). Many of these 2-D channel flow states except the vortex-roll state are also found in active nematic models (Shendruk et al. Reference Shendruk, Doostmohammadi, Thijssen and Yeomans2017; Chandragiri et al. Reference Chandragiri, Doostmohammadi, Yeomans and Thampi2019). In the unidirectional flow state, the flow rate
$Q$
increases as activity
$\zeta$
increases. Upon further increase of
$\zeta$
,
$Q$
drops abruptly as the system transitions into the travelling-wave state, in which the spontaneous flow starts to develop in the vertical direction. Owing to the random initial conditions of the stress tensor field, rightward and leftward flows in the unidirectional and travelling-wave states can be observed in our simulation with equal probability (see supplementary material). Based on the linear stability analysis (see supplementary material), the wavelength
$l_{tw}$
of the travelling-wave state increases with
$\textrm {Wi}$
and decreases with
$\zeta$
as shown in figures 2(a)and 2(b). In the vortex-roll state, the stretching of the polymers gives rise to a birefringent field, in which
$\pm 1/2$
defects can be identified (see figure 1
e), and this resembles active nematic systems (Zhang et al. Reference Zhang, Mozaffari and de Pablo2021). The
$\pm 1/2$
topological defects refer to regions where the principal axis corresponding to the largest eigenvalue of the local stress tensor changes abruptly, which is similar to that seen in nematic LCs. The development of the local nematic order has also been observed in 3-D microtubule suspensions (Wu et al. Reference Wu, Hishamunda, Chen, DeCamp, Chang, Fernández-Nieves, Fraden and Dogic2017).
Figure 2. Effect of (a) Weissenberg number
$\textrm {Wi}$
and (b) activity
$\zeta$
on the wavelength
$l_{tw}$
in the travelling-wave state in the 2-D channel by linear stability analysis (see supplementary material). Here
$\zeta =3.4$
in (a) and
$\textrm {Wi}=10^{0.5}$
in (b). (c) Normalised velocity-velocity correlation function
$C(r)$
at different activity
$\zeta$
. The inset shows the velocity correlation length
$l_v$
determined by
$C=0.8$
. Region I denotes the vortex-roll state, and region II represents the dancing and turbulent states. (d) Effect of
$\zeta$
on the r.m.s. velocity,
$V_{rms}$
. The error bars denote the standard deviation for time series of r.m.s. velocity. In (c,d), the Weissenberg number is
$\textrm {Wi}= 10^{-0.5}$
.
Figure 2(c) shows the effect of activity
$\zeta$
on the normalised velocity-velocity correlation function
$C(r)=\langle \textbf {U}(r)\cdot \textbf {U}(0) \rangle /\langle \textbf {U}(0)^2 \rangle$
, where
$r$
denotes the dimensionless horizontal distance and
$\langle \cdot \rangle$
denotes the ensemble (spatio–temporal) average. The velocity correlation length
$l_v$
is found to be independent of
$\zeta$
except at the transition from the vortex-roll into the turbulent state, at which
$l_v$
undergoes an abrupt decrease. Such a feature is also observed in the active microtubule experiment (Sanchez et al. Reference Sanchez, Chen, DeCamp, Heymann and Dogic2012) and in nematic LC models (Thampi, Golestanian & Yeomans Reference Thampi, Golestanian and Yeomans2013). Additionally, the measured root-mean-square (r.m.s.) velocity
$V_{rms}$
saturates at high
$\zeta$
as presented in figure 2(d) which is consistent with the experiments (Sanchez et al. Reference Sanchez, Chen, DeCamp, Heymann and Dogic2012; Henkin et al. Reference Henkin, DeCamp, Chen, Sanchez and Dogic2014; Tan et al. Reference Tan, Roberts, Smith, Olvera, Arteaga, Fortini, Mitchell and Hirst2019), while the nematic LC model did not exhibit such saturation trend (Thampi et al. Reference Thampi, Golestanian and Yeomans2013).
3.2. Annulus geometry
In the annulus geometry, the aspect ratio is defined by
$\Gamma =R_i/R_o$
, where
$R_i$
and
$R_o$
denote the inner and outer radii, respectively. The grid resolution is set to
$N_{\theta }\times N_r=216\times 36$
. Here the outer radius is chosen as the characteristic length, namely
$L_c=R_o$
. This geometry parameter
$\Gamma \in (0,1)$
can be used to investigate the effect of geometric confinement on the flow states. Different from the channel geometry, the annulus geometry exhibits stationary, unidirectional, travelling-wave and turbulent states only as shown in figure 3(a), consistent with the previous numerical results which were obtained for 2-D apolar active suspension confined in an annulus by using a coarse-grained LC model (Chen, Gao & Gao Reference Chen, Gao and Gao2018). As
$\Gamma$
increases, the confinement becomes stronger, and the system undergoes transitions from the turbulent flow into the travelling-wave state and into the unidirectional flow state (see figure 3
a,b). Figure 3(b) shows that for
$\Gamma \to 1$
(thin annulus) or
$\Gamma \to 0$
(wide annulus), the system is in the stationary or the turbulent state, in which a net flow vanishes. Note that in a 3-D experiment of a microtubule-based active fluid, the flow rate also shows a non-monotonic dependence on the aspect ratio of the confinement (Wu et al. Reference Wu, Hishamunda, Chen, DeCamp, Chang, Fernández-Nieves, Fraden and Dogic2017). Because the azimuthal direction of the annulus requires periodicity, the wavelength of the travelling-wave state should be a fraction of the lateral length. As shown in figure 3(c), it is observed that such wavelength
$\lambda _a$
increases with the increasing width
$1-\Gamma$
, consistent with experimental results (Chandrakar et al. Reference Chandrakar, Varghese, Aghvami, Baskaran, Dogic and Duclos2020).
Figure 3. Flow state transitions of an active viscoelastic fluid in a 2-D annulus. (a) Flow state diagram in terms of
$\zeta$
and aspect ratio
$\Gamma$
at
$\textrm {Wi}=1$
. (b) Effect of aspect ratio
$\Gamma$
on the average flow rate at
$\zeta =5.2$
. The error bars denote the standard deviation of time series for flow rate. (c) Variation of the arc wavelength
$\lambda _a$
with the thickness
$(1-\Gamma )$
of the annulus at
$\zeta =4.4$
. The inset contour denotes the velocity magnitude under different aspect ratio
$\Gamma$
, which helps to identify the number of repeated arc segments,
$N_a$
.
Figure 4. Chirality-switching vortex flows of the active viscoelastic fluid in a 2-D disk. (a) Chirality-switching frequency
$f$
as a function of
$\textrm {Wi}$
at
$\zeta =4.5$
. The friction coefficient is set to
$\nu =1.1$
. The bottom-left inset shows the time series of
$Q$
at
$\textrm {Wi}=2.4$
, and the top-right inset shows the transient two-vortex structure in our model and in former experiment (Liu et al. Reference Liu, Shankar, Marchetti and Wu2021). (b) Temporal evolution of viscous dissipation rate and frictional dissipation rate. (c) Sequential snapshots show the chirality switching from CW to CCW. The colour indicates velocity magnitude.
3.3. Disk geometry
Next we study disk confinement where its radius is taken as the characteristic length
$L_c$
, and the grid resolution is
$N_{\theta }\times N_r=216\times 36$
. Different from the previous two geometries, a disk only has one single wall which may induce more interesting flow phenomena. Similar to that in a channel, the active flow in a disk also tends to be turbulent at high
$\zeta$
and small
$\textrm {Wi}$
(see supplementary material for the flow state diagram of the disk geometry). This effect of
$\textrm {Wi}$
agrees well with the experimental result for low DNA concentration (and hence small
$\textrm {Wi}$
) reported in Liu et al. (Reference Liu, Shankar, Marchetti and Wu2021). The effective friction due to the confining substrates can tune the flow structure in active fluids (Liu et al. Reference Liu, Shankar, Marchetti and Wu2021; Caballero et al. Reference Caballero, You and Marchetti2023). Here we incorporate it into our model by setting
$\nu =1.1$
. We observe the flow which shows periodic chirality switching, with one full-scale vortex followed by the next one with opposite chirality (see figure 4(c) and supplementary movie 2). Moreover, the chirality switching shows a frequency that decreases as
$\textrm {Wi}$
increases, which is proportional to the polymer relaxation time as shown in figure 4(a). The detailed process of switching from CW to CCW vortex is presented in figure 4(c). The system is started with a full-scale CW vortex at
$t_1$
, which gradually decays (at
$t_2$
), while a secondary CCW vortex is formed close to the wall (at
$t_3$
). The decay of the CW vortex progresses. Concomitantly, the CCW vortex continues to grow (at
$t_4$
and
$t_5$
). Finally, the CCW vortex reaches its full scale (at
$t_6$
). It is interesting to note that the centres of the vortices also rotate, first in the CW sense (before
$t_4$
) and then in the CCW sense (after
$t_4$
). This chirality switching repeats periodically. These numerical observations agree well with the experiment of bacterial active fluids, where the switching can be made slower by adding more DNAs into the fluid (Liu et al. Reference Liu, Shankar, Marchetti and Wu2021). A similar chirality-reversal vortex flow was observed in a microtubule-based active nematic confined to a disk region; however, the reversal does not appear to persist periodically (Opathalage et al. Reference Opathalage, Norton, Juniper, Langeslay, Aghvami, Fraden and Dogic2019). Figure 4(b) shows the frictional dissipation rate
$D_f=\nu \int |\textbf {U}|^2\textrm {d}\textbf {x}$
and the viscous dissipation rate
$D_v=\frac {1}{2}\int [ \nabla \textbf {U}+(\nabla \textbf {U})^T ]: [ \nabla \textbf {U}+(\nabla \textbf {U})^T ]\textrm {d}\textbf {x}$
, and their time variations are consistent with a fast decrease of the absolute flow rate
$|Q|$
from a peak (from
$t_1$
to
$t_3$
) followed by a relatively slow increase of
$|Q|$
toward the next peak (from
$t_4$
to
$t_6$
). A proper orthogonal decomposition (POD) analysis (Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993) of the vorticity field reveals that the dominant modes are 1-, 2- and 4-vortex flows as shown in figure 5 which is consistent with a recent bacterial active fluid experiment (Nishiguchi et al. Reference Nishiguchi, Shiratani, Takeuchi and Aranson2024). Using different
$\zeta$
and
$\textrm {Wi}$
, our system can generate chirality-switching flows with more vortices (see supplementary material). As
$\textrm {Wi}$
decreases, a transition from 2-vortex to 4-vortex state can take place, again in good agreement with the bacteria experiment (Nishiguchi et al. Reference Nishiguchi, Shiratani, Takeuchi and Aranson2024). A similar periodic flow reversal phenomenon was predicted in a channel of an active nematic sandwiched between two viscoelastic layers (Mori et al. Reference Mori, Bhattacharyya, Yeomans and Thampi2023).
Figure 5. Energy of different modes for chirality-switching vortex flows of the active viscoelastic fluid in a 2-D disk with
$\nu =1.1$
,
$\zeta =4.5$
and
$\textrm {Wi}=2.4$
. Inset shows the first three modes obtained by the POD analysis.
4. Conclusion
We have developed an active viscoelastic liquid model to describe active liquids without LC order. We studied its spontaneous flow patterns in 2-D channel, annulus and disk geometries. In addition to the geometric confinement (i.e. aspect ratio of annulus here), the current model is also capable of controlling the flow state by Weissenberg number which is one macroscopic parameter measuring the viscoelasticity of the liquid. Our model features a different mathematical form compared with existing models, and it can reproduce active-flow patterns and characteristics typical of distinct experimental systems within a single framework. Therefore, we argue that flow patterns observed in bacteria-based systems, e.g. periodic chirality-switching flow, should also occur in microtubule-based systems, and vice versa: the non-monotonic dependence of channel flow on aspect ratio found in microtubule-based systems should also be found in bacteria-based systems (Wioland et al. Reference Wioland, Lushi and Goldstein2016). Note that for gel-like active LCs exhibiting a nematic or polar order (Hemingway et al. Reference Hemingway, Maitra, Banerjee, Marchetti, Ramaswamy, Fielding and Cates2015, Reference Hemingway, Cates and Fielding2016), the well-developed active LC model will be needed. In addition, the present 2-D simulations are unable to explicitly capture the inherently 3-D phenomena observed in isotropic active fluids except for quasi-2-D phenomena which can be reproduced through the addition of a frictional term to account for the viscous damping by the confining walls in the third direction. In future work, we plan to extend our research to 3-D systems, leveraging the inherently 3-D nature of the PTT model used in current work. Furthermore, exploring the influence of pre-defined polymer orientations at the boundaries also offers a compelling direction for further investigation.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.177.
Funding
R.Z. acknowledges the Hong Kong RGC General Research Fund (grant no. 16300221). T.Q. acknowledges the Hong Kong RGC General Research Fund (grant no. 16306121) and the Key Project of the National Natural Science Foundation of China (no. 12131010).
Declaration of interests
The authors report no conflict of interest.
Open access
