When a solution of polymers flows past obstacles or through a porous medium, the flexible chains that are transported undergo strong velocity gradients that change their conformation and generate elastic stress (Dauben & Menzie Reference Dauben and Menzie1967; Seright et al. Reference Seright, Fan, Wavrik and de Carvalho Balaban2011; Kawale et al. Reference Kawale, Bouwman, Sachdev, Zitha, Kreutzer, Rossen and Boukany2017a; Browne, Shih & Datta Reference Browne, Shih and Datta2020b). In turn, this stress modifies the velocity field, leading to a two-way coupling that may induce strong nonlinear effects, such as elastic turbulence (Groisman & Steinberg Reference Groisman and Steinberg2000; Steinberg Reference Steinberg2021). Contrary to inertial turbulence, nonlinear viscoelastic effects are predominant for the smallest pores (Poole Reference Poole2019; Browne et al. Reference Browne, Shih and Datta2020b), not the largest. Modern experimental approaches have thus exploited microfluidic technologies with pore sizes of several tens to hundreds of micrometres to explore elastic instabilities in complex structures (Kawale et al. Reference Kawale, Marques, Zitha, Kreutzer, Rossen and Boukany2017b).
The most recent microfluidic experiments have used slender cylinders to generate quasi-two-dimensional (quasi-2-D) flows that can be visualized using standard optical microscopy (Haward, Toda-Peters & Shen Reference Haward, Toda-Peters and Shen2018; Hopkins, Haward & Shen Reference Hopkins, Haward and Shen2020; Haward et al. Reference Haward, Hopkins, Varchanis and Shen2021a). This technology has revealed fascinating physics. Spontaneous symmetry breaking and bistability have been shown to occur in the flow of wormlike micellar or hydrolysed polyacrylamide solutions past a single cylinder confined in a channel (Haward et al. Reference Haward, Kitajima, Toda-Peters, Takahashi and Shen2019; Haward, Hopkins & Shen Reference Haward, Hopkins and Shen2020; Hopkins et al. Reference Hopkins, Haward and Shen2020). In the case of two cylinders aligned with the flow, a stagnation zone with counter-rotating recirculation vortices appears between the cylinders (Varshney & Steinberg Reference Varshney and Steinberg2017, Reference Varshney and Steinberg2019; Kumar & Ardekani Reference Kumar and Ardekani2021). Flows past two side-by-side cylinders have revealed a series of bifurcations depending on both the gap between cylinders and the Weissenberg number, with symmetric and asymmetric, converging and diverging, bi- and tri-stable states (Hopkins, Haward & Shen Reference Hopkins, Haward and Shen2021; Khan & Sasmal Reference Khan and Sasmal2021).
Further increasing in complexity, arrays of slender cylinders (Walkama, Waisbord & Guasto Reference Walkama, Waisbord and Guasto2020; Haward, Hopkins & Shen Reference Haward, Hopkins and Shen2021b) have been used to explore elastic instabilities and transition to chaos in porous structures. Walkama et al. (Reference Walkama, Waisbord and Guasto2020) showed that the spatio-temporal velocity fluctuations that develop in a staggered hexagonal lattice of cylinders are suppressed by randomly displacing each cylinder from the crystalline position. They used this observation to argue that disorder suppresses chaos. Haward et al. (Reference Haward, Hopkins and Shen2021b) showed that introducing disorder into an aligned hexagonal lattice – the set of points that defines the aligned structure being congruent with the staggered one through a rotation of 
${{\rm \pi} }/{6}$ modulo 
${{\rm \pi} }/{3}$ – yields the opposite effect, with disorder promoting instability. They hypothesized that the important parameter that controls the flow is not disorder, but rather the number of stagnation points accessible to the polymeric chains, with more stagnation points leading to less stability.
Stagnation points induce large molecular strain and represent starting points for the formation of localized regions of large polymer stress known as birefringent strands (Becherer, van Saarloos & Morozov Reference Becherer, van Saarloos and Morozov2009). Strands further induce a feedback on the flow field, essentially acting as a line distribution of forces within an otherwise Newtonian fluid (Harlen, Rallison & Chilcott Reference Harlen, Rallison and Chilcott1990). We have recently shown in Mokhtari et al. (Reference Mokhtari, Latché, Quintard and Davit2022) that these structures play a fundamental role in steady viscoelastic flows past obstacles. In arrays of cylinders, they modify the global flow distribution with an increase of stagnation zones and strong channelization. This yields an increase of viscous dissipation in the flow channels and of the rate of entropy generation in the strands. Both these mechanisms generate a surge in resistance to flow through the structure – a phenomenon observed experimentally but not fully understood yet (Galindo-Rosales et al. Reference Galindo-Rosales, Campo-Deano, Pinho, Van Bokhorst, Hamersma, Oliveira and Alves2012; Clarke et al. Reference Clarke, Howe, Mitchell, Staniland, Hawkes and Leeper2015, Reference Clarke, Howe, Mitchell, Staniland and Hawkes2016; Mitchell et al. Reference Mitchell, Lyons, Howe and Clarke2016; Qin et al. Reference Qin, Salipante, Hudson and Arratia2019; Browne & Datta Reference Browne and Datta2021). This evidence that strands are a fundamental component of steady viscoelastic flows past obstacles raises important new questions: What are the feedback mechanisms between flow in these complex geometries, stagnation points on solid obstacle and birefringent strands? Could localized structures of stress, such as strands, be a key to a better understanding of viscoelastic instabilities and chaos in complex structures?
In this article, we use high-performance computing to simulate the flow of viscoelastic fluids through a 2-D hexagonal lattice of obstacles at zero Reynolds number. We solve incompressible Oldroyd-B, FENE-CR (Finite Extendable Non-linear Elastic – Chilcott and Rallison) and FENE-P (Finite Extendable Non-linear Elastic – Peterlin) models with an imposed force density in the momentum transport equation and biperiodic boundary conditions. By varying the orientation of the forcing term, the Weissenberg number 
$Wi$ and the ratio of viscosities 
${\beta }$, we explore the physics of viscoelastic flows in a variety of configurations, including both the aligned and staggered cases.
1. Methods
1.1. Domain and boundary conditions
We consider a simulation domain that is a rectangle with biperiodic left/right and top/bottom conditions and a no-slip/no-penetration condition on the boundaries of the obstacles. We have two geometries with a hexagonal lattice of cylindrical obstacles (figure 1). The first one, which is used for studying multistability and steady flow, consists of 
$30$ cylinders in the aligned configuration. The second one, which is used for unsteady flow, is larger and consists of 
$120$ cylinders in the staggered configuration.
Figure 1. Geometry of the hexagonal lattice of obstacles with aligned and staggered configurations, coordinate systems and angle definitions. (a) Aligned configuration and corresponding coordinate system 
$(\boldsymbol {x},\boldsymbol {y})$. The radius of each circle is 
$R$ and the shortest distance between centres of obstacles is 
$S$. The unit vector 
$\boldsymbol {f}$ is the force density in the non-dimensionalized momentum transport equation and forms an angle 
$\theta$ with 
$\boldsymbol {x}$. Here, 
$\langle \boldsymbol {u}\rangle$ is the intrinsic average of the velocity field 
$\langle \boldsymbol {u}\rangle ={1}/{|\varOmega |}\int _{\varOmega }\boldsymbol {u}\,{\rm d}\varOmega$ and forms an angle 
$\alpha$ with 
$\boldsymbol {x}$. (b) Staggered configuration and corresponding coordinate system 
$(\boldsymbol {x}^{\star },\boldsymbol {y}^{\star })$. The staggered periodic pattern can be obtained by rotating the aligned one by an angle of 
$30^\circ$.
1.2. Flow equations
We consider an incompressible flow with 
$\text {div}\,\boldsymbol {\mathfrak {u}}=0$ and steady momentum transport given by
\begin{equation} 0=-\boldsymbol{\nabla}\mathfrak{p}+\text{div}\left(\tau_{s}+\tau_{p}\right)+\boldsymbol{{\mathfrak{F}}}, \end{equation}
with 
$\boldsymbol {\mathfrak {u}}$ the velocity, 
$\mathfrak {p}$ the pressure, 
$\tau _{s}$ the solvent stress tensor, 
$\tau _{p}$ the polymer stress tensor and 
$\boldsymbol {\mathfrak {F}}$ the imposed force density. This choice of imposing a pressure gradient through 
$\boldsymbol {\mathfrak {F}}$, rather than through the boundary conditions, is natural when working with periodic boundary conditions and eliminates boundary effects or instabilities specific to Dirichlet and Neumann conditions – the idea is that of an ‘infinite’ porous medium.
The constitutive stress relation for the solvent is
\begin{equation} \tau_{s}=\eta_{s}\left(\boldsymbol{\nabla}\boldsymbol{\mathfrak{u}}+\left(\boldsymbol{\nabla} \boldsymbol{\mathfrak{u}}\right)^{\mathsf{T}}\right) \end{equation}
with 
$\eta _{s}$ the viscosity. The general form of the polymer stress tensor is
\begin{equation} \tau_{p}=\frac{\eta_{p}}{\lambda}\left(g\left(\boldsymbol{c}\right)\boldsymbol{c}- h\left(\boldsymbol{c}\right)\boldsymbol{I}_{d}\right) \end{equation}
with 
$\eta _{p}$ the polymer viscosity, 
$\lambda$ a characteristic time for polymer relaxation and 
$\boldsymbol {c}$ the conformation tensor. We consider the Oldroyd-B, FENE-CR and FENE-P models with
\begin{equation} g(\boldsymbol{c})=h(\boldsymbol{c})=1, \end{equation}for the Oldroyd-B,
\begin{equation} g(\boldsymbol{c})=h(\boldsymbol{c})=\frac{b}{b-\mathrm{tr}(\boldsymbol{c})}, \end{equation}for the FENE-CR and
\begin{equation} g(\boldsymbol{c})=\frac{b}{b-\mathrm{tr}(\boldsymbol{c})},\quad h(\boldsymbol{c})=\frac{b}{b-\mathrm{tr}(\boldsymbol{I}_{d})}, \end{equation}
for the FENE-P. In both FENE models, 
$b$ is the parameter that controls the maximum stretching of the polymers with 
$\mathrm {tr}(\boldsymbol {c})< b$. The Oldroyd-B model is based upon an affine relation between 
$\tau _{p}$ and 
$\boldsymbol {c}$, which can be interpreted as a Hookean entropic spring description in a molecular dumbbell representation (Bird et al. Reference Bird, Curtiss, Armstrong and Hassager1987). This model imposes no limit upon the stretching of the molecules and can lead to stress singularity (Shaqfeh & Khomami Reference Shaqfeh and Khomami2021). FENE-type models limit the trace of the conformation tensor, thus limiting the maximum stretching of the chains (Bird et al. Reference Bird, Curtiss, Armstrong and Hassager1987). FENE-CR and FENE-P models differ in their response to shear flow: the FENE-CR model represents a Boger fluid, whereas the FENE-P model is shear thinning (Herrchen & Öttinger Reference Herrchen and Öttinger1997).
The last constitutive equation is the transport of the conformation tensor, which reads
\begin{equation} \partial_{t}\boldsymbol{c}+(\boldsymbol{\mathfrak{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{c}=\boldsymbol{\nabla}\boldsymbol{\mathfrak{u}}\,\boldsymbol{c}+\boldsymbol{c} (\boldsymbol{\nabla}\boldsymbol{\mathfrak{u}})^{\mathsf{T}}-\frac{1}{\eta_{p}}\,\tau_{p}. \end{equation}1.3. Non-dimensionalization
To non-dimensionalize the problem, we use the minimum distance between two obstacles, 
$H=S-2R$, as a characteristic length. We also use the expression
\begin{equation} U=\frac{H^{2}}{\eta_{s}+\eta_{p}}\Vert\boldsymbol{\mathfrak{F}}\Vert, \end{equation}
for the reference velocity and 
${H}/{U}$ as a characteristic time. This reference velocity can be understood as an estimation of the Darcy velocity in the Newtonian case, with 
$H^{2}$ an estimate of the permeability, 
$\eta _{s}+\eta _{p}$ the viscosity and 
$\Vert \boldsymbol {\mathfrak {F}}\Vert$ a norm of the pressure gradient.
The final dimensionless equations read
$$\begin{gather} 0 =-\boldsymbol{\nabla} p+\nabla^{2}\boldsymbol{u}+\frac{\beta}{Wi}\text{div}\left(g\left(\boldsymbol{c}\right)\boldsymbol{c}-h \left(\boldsymbol{c}\right)\boldsymbol{I}_{d}\right)+(1+\beta)\boldsymbol{f}, \end{gather}$$
$$\begin{gather}0 =\text{div}\,\boldsymbol{u}, \end{gather}$$
$$\begin{gather}\partial_{t}\boldsymbol{c}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{c} = \boldsymbol{\nabla}\boldsymbol{u}\,\boldsymbol{c}+\boldsymbol{c}(\boldsymbol{\nabla}\boldsymbol{u})^{\mathsf{T}}- \frac{1}{Wi}\left(g\left(\boldsymbol{c}\right)\boldsymbol{c}-h(\boldsymbol{c})\boldsymbol{I}_{d}\right), \end{gather}$$with the viscosity ratio defined as
\begin{equation} \beta=\frac{\eta_{p}}{\eta_{s}}, \end{equation}the Weissenberg number as
\begin{equation} Wi=\frac{\lambda H}{\eta_{s}+\eta_{p}}\Vert\boldsymbol{\mathfrak{F}}\Vert,\end{equation}and the pressure, velocity and force density as
\begin{equation} p=\frac{1+\beta}{\Vert\boldsymbol{\mathfrak{F}}\Vert}\mathfrak{p}, \quad \boldsymbol{u}=\frac{\mathfrak{u}}{U}, \quad \boldsymbol{f}=\frac{\boldsymbol{\mathfrak{F}}}{\Vert\boldsymbol{\mathfrak{F}}\Vert}. \end{equation} It is important to note that, with this definition of the reference velocity, the Weissenberg number 
$Wi$ does not depend on the actual average flow velocity, but rather on the volume source term 
$\Vert \boldsymbol {\mathfrak {F}}\Vert$. This is a natural choice in our setting since flow results from this volume source term and not from an imposed flow on the boundary conditions. Working with the average flow velocity is possible but unpractical as it is not known a priori but rather results from each simulation. One caveat of our definition in (1.13) is that comparison with transition values of 
$Wi$ from works using the more standard definition is not straightforward – keeping in mind that no comparison with experiments is ever straightforward as viscoelastic models are widely accepted as not being accurate in estimating these transition values. Our largest 
$Wi$ values, for example, would be much smaller with a definition based on the average flow velocity. For comparison with other works, we provide in the supplementary material available at https://doi.org/10.1017/jfm.2023.916 the equivalent numbers for a more classical definition of the Weissenberg number.
1.4. Numerical scheme and implementation
A detailed presentation of the numerical scheme is available in Mokhtari et al. (Reference Mokhtari, Davit, Latché and Quintard2023). Momentum transport is solved via a staggered grid scheme using the Rannacher and Turek low-order finite elements (Rannacher & Turek Reference Rannacher and Turek1992). The discrete conformation tensor unknowns are located at the centres of the cells and the corresponding transport equation (1.11) is solved using a finite-volume scheme.
The mesh is structured and uniform. It consists of quadrilateral cells with obstacles obtained by making holes into the uniform mesh and circle boundaries approximated as stair steps (see Mokhtari et al. (Reference Mokhtari, Latché, Quintard and Davit2022) and more details in the supplementary material). In studying the transition to unsteady flows in the staggered case, we took care of working with structured meshes that were symmetric about the horizontal configuration of the strands corresponding to base flow at steady state. We used a constant mesh size 
${\simeq }0.04R$, where 
$R$ is the circle radius. This corresponds to approximately 
$3\times 10^{5}$ cells and approximately 
$1.2$ million cells for the two geometries presented in figure 1.
At time 
$t=0$, the velocity, pressure and conformation fields are initialized as 
$\boldsymbol {u}=(0\ 0)^{\mathsf {T}}$, 
$p=0$ and 
$\boldsymbol {c}=\boldsymbol {I}_{d}$. The time discretization consists of a fractional step approach where (1.9) and (1.10) are decoupled from (1.11) using a beginning-of-step approximation of the latter. Given 
$\boldsymbol {c}^{n}$, the value of the conformation tensor at time iteration 
$n$, we calculate the pressure and velocity at time iteration 
$n+1$ by solving
$$\begin{gather} 0 =-\boldsymbol{\nabla} p^{n+1}+\nabla^{2}\boldsymbol{u}^{n+1}+\boldsymbol{S}^{n}, \end{gather}$$
$$\begin{gather}0 =\text{div}\,\boldsymbol{u}^{n+1}, \end{gather}$$
with 
$\boldsymbol {S}^{n}=({\beta }/{Wi})\text {div}(g(\boldsymbol {c}^{n})\boldsymbol {c}^{n}- h(\boldsymbol {c}^{n})\boldsymbol {I}_{d})+(1+\beta )\boldsymbol {f}$ the source term. For each time step, this Stokes problem is solved in an iterative way through a standard projection scheme that decouples the momentum balance equation from the divergence constraint. For each sub-iteration 
$k+1$, the method makes an initial prediction 
$\tilde {\boldsymbol {u}}_{k+1}$ of the velocity using the pressure gradient 
$\boldsymbol {\nabla } p_{k}$, starting from 
$\boldsymbol {u}_{0}={{\boldsymbol {u}^{n}}}$ and 
$p_{0}=p^{n}$. The correction step then consists in determining 
$p_{k+1}$ and correcting 
$\boldsymbol {u}_{k+1}$ in such a way that the incompressibility condition (1.16) is respected. This prediction–correction scheme reads
\begin{equation} \left.\begin{gathered} {Prediction\ step}-\text{Solve for }{\tilde{\boldsymbol{u}}_{k+1}}:\\ \xi\left(\tilde{\boldsymbol{u}}_{k+1}-\boldsymbol{u}_{k}\right)+\boldsymbol{\nabla} p_{k}-\nabla^{2}\tilde{\boldsymbol{u}}_{k+1}+\boldsymbol{S}^{n}=0.\\ {Correction\ step}-\text{Solve for }p_{k+1}\text{ and }{\boldsymbol{u}_{k+1}:}\\ \xi\left(\boldsymbol{u}_{k+1}-\tilde{\boldsymbol{u}}_{k+1}\right)+\boldsymbol{\nabla}\left(\,p_{k+1}-p_{k}\right)=0,\\ \text{div}\,\boldsymbol{u}_{k+1}=0, \end{gathered}\right\} \end{equation}
where 
$\xi >0$ is a constant that can be understood as the inverse of a time step. Convergence is considered to be reached when the relative difference between two consecutive sub-iterations satisfies the following criterion:
\begin{equation} \max\left(\frac{\left\Vert \boldsymbol{u}_{k+1}-\boldsymbol{u}_{k}\right\Vert _{\infty}}{\max\left(1,\left\Vert \boldsymbol{u}_{k}\right\Vert _{\infty}\right)},\frac{\left\Vert \,p_{k+1}-p_{k}\right\Vert _{\infty}}{\max\left(1,\left\Vert \,p_{k}\right\Vert _{\infty}\right)}\right)\leq10^{-8}, \end{equation}which is small enough to guarantee the independence of the global solution and is easily obtained after a few sub-iterations.
The transport equation of the conformation tensor (1.11) is solved using a Strang operator splitting combined, for accuracy, with a 
$\log$ transform of the hyperbolic transport equation
\begin{equation} \left.\begin{gathered} {Advection \; I}-\text{Solve for }{\boldsymbol{c}^{n+({1}/{3})}:}\\ \dfrac{1}{{\rm d}t/2}\left(\log\mathsf{\boldsymbol{c}}^{n+({1}/{3})}-\log\boldsymbol{c}^{n}\right)+ \left(\boldsymbol{u}^{n+1}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\log \boldsymbol{c}^{n+({1}/{3})}=0.\\ {ODE}-\text{Set }{\boldsymbol{c}(t_{n})=\boldsymbol{c}^{n+({1}/{3})}}\text{ and solve for } \boldsymbol{c}^{n+({2}/{3})}=\boldsymbol{c}(t_{n}+{\rm d}t):\\ \boldsymbol{\dot{c}}-\left(\boldsymbol{\nabla}\boldsymbol{u}^{n+1}\right) \boldsymbol{c}-\boldsymbol{c}\left(\boldsymbol{\nabla}\boldsymbol{u}^{n+1}\right)^{\mathsf{T}}+ \frac{1}{Wi}\left(g\left(\boldsymbol{c}\right)\boldsymbol{c}-h \left(\boldsymbol{c}\right)\boldsymbol{I}_{d}\right)=0.\\ {Advection\; II}-\text{Solve for }{\boldsymbol{c}^{n+1}:}\\ \dfrac{1}{{\rm d}t/2}\left(\log\boldsymbol{c}^{n+1}-\log\boldsymbol{c}^{n+({2}/{3})}\right)+ \left(\boldsymbol{u}^{n+1}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\log\boldsymbol{c}^{n+1}=0, \end{gathered}\right\} \end{equation}
with 
$\textrm {d}t$ the time step and ODE the ordinary differential equation. In practice, 
$dt$ is adaptive and monitored in the interval 
$[10^{-6},10^{-3}]$ in order to reduce the total computation time, while always satisfying the following criterion:
\begin{equation} \max\left(\frac{\left\Vert \boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}\right\Vert _{\infty}}{\max\left(1,\left\Vert \boldsymbol{u}^{n}\right\Vert _{\infty}\right)},\frac{\left\Vert \,p^{n+1}-p^{n}\right\Vert _{\infty}}{\max\left(1,\left\Vert \,p^{n}\right\Vert _{\infty}\right)},\frac{\left\Vert \boldsymbol{c}^{n+1}-\boldsymbol{c}^{n}\right\Vert _{\infty}}{\max\left(1,\left\Vert \boldsymbol{c}^{n}\right\Vert _{\infty}\right)}\right)\leq10^{-4}. \end{equation}For each cell, the local ordinary differential equation corresponding to the relaxation is solved using a first-order implicit Euler scheme with local sub-cycling
\begin{align} \frac{1}{\delta t}\left(\boldsymbol{c}_{k+1}-\boldsymbol{c}_{k}\right)=\boldsymbol{\nabla}\boldsymbol{u}{}^{n+1} \boldsymbol{c}_{k+1}+\boldsymbol{c}_{k+1}\left(\boldsymbol{\nabla}\boldsymbol{u}^{n+1}\right)^{\mathsf{T}}- \frac{1}{Wi}\left(g\left(\boldsymbol{c}_{k+1}\right)\boldsymbol{c}_{k+1}-h \left(\boldsymbol{c}_{k+1}\right)\boldsymbol{I}_{d}\right),\end{align}
where 
$\delta t$ is a local time step on each mesh element which is set to 
$\delta t={\textrm {d}t}/{a}$, with 
$a$ the smallest integer number such that 
$\delta t\le 1/(100\Vert \boldsymbol {\nabla }\boldsymbol {u}^{n+1}\Vert _{\infty })$, to preserve the positive definiteness of the conformation tensor; see Mokhtari et al. (Reference Mokhtari, Davit, Latché and Quintard2023).
The solver is implemented as a specific module of CALIF
$^{3}$S (CALIF3S 2021). Calculations were performed on TotalEnergies’ supercomputer Pangea II. Parallel computations were carried out using a domain decomposition approach based on the METIS (4.0.3) graph partitioner (Karypis & Kumar Reference Karypis and Kumar1998) and the library OpenMPI (3.1.5) (Graham, Woodall & Squyres Reference Graham, Woodall and Squyres2005). As an example to illustrate the cost of the calculations, one point in the stationary case required approximately 
$10^{3}\ \textrm {cores}\times \textrm {hours}$ while the case (
$\beta =10, Wi=84$) alone required approximately
$10^{5}\ \textrm {cores}\times \textrm {hours}$.
1.5. Time averaging and root mean square variations
The time-averaged value of a variable 
$\varphi$ over the interval 
$[t_{1},t_{2}]$ is defined as
\begin{equation} \bar{\varphi}=\frac{1}{t_{2}-t_{1}}\int_{t_{1}}^{t_{2}}\varphi(\tau)\,{\rm d}\tau, \end{equation}and the corresponding root mean square (r.m.s.) is defined as
\begin{equation} \varphi_{rms}=\sqrt{\overline{(\varphi-\bar{\varphi})^{2}}}. \end{equation}1.6. Kymographs
For the kymographs, the variable 
$\psi$ is averaged along the vertical direction using the operator 
$\langle \psi \rangle ^{y}=({1}/{L_{y}})\int _{0}^{L_{y}}\psi (x,y)\,{\textrm {d}y}$, where 
$L_{y}$ corresponds to the domain height. The spatial average is then normalized by the time average as 
$({\langle \psi \rangle ^{y}-\langle \bar {\psi }\rangle ^{y}})/{\langle \bar {\psi }\rangle ^{y}}$. Kymographs display this field – averaged in the 
$y$ direction – along the horizontal direction over time.
1.7. Power spectral density
The power spectral density 
$S(\,f)$ of a signal 
$\varphi$ discrete in time is given by
\begin{equation} S(\,f)=\sum_{k=1}^{K}A(k)\exp(-2{\rm \pi} {\rm i}fk), \end{equation}
with 
$K$ the number of discrete values. This corresponds to the discrete Fourier transform of the auto-correlation function 
$A(k)$ defined by
\begin{equation} A(k)=\frac{1}{K}\sum_{n=1}^{K}\varphi(n)\varphi(n-k). \end{equation}2. Results
2.1. Steady base flows
Our first objective is to study how steady-state strands generated in the wake of an obstacle (Haward et al. Reference Haward, Toda-Peters and Shen2018) interact with other obstacles to guide the flow through the lattice. Since strands direct the flow through the porous structures (Mokhtari et al. Reference Mokhtari, Latché, Quintard and Davit2022) and thus define the orientation of the flow, our main observable is the angle 
$\alpha$ between the average velocity field 
$\langle \boldsymbol {u}\rangle$ and the 
$\boldsymbol {x}$ axis, which is a proxy for the orientation of the strands. Our idea is to study how the system responds to changes in the forcing term, first starting in figure 2(a) with variations of the Weissenberg number 
$Wi$ for fixed values of the angle 
$\theta$ of the force density 
$\boldsymbol {f}$ (as defined in figure 1). In our simulations, we proceeded by progressively increasing the Weissenberg number, starting from 
$Wi=0$, for each value of 
$\theta$. The case 
$Wi=0$ is Newtonian, so that the observed flow angle is the same as the angle of the imposed force density, 
$\theta =\alpha (Wi=0)$. Upon increasing 
$Wi$, results in figure 2(a) indicate the clear formation of two preferential flow directions at 
$\alpha =0^{\circ }$ and 
$\alpha =30^{\circ }$ for three different constitutive relations (Oldroyd-B, FENE-P and FENE-CR). Figure 3 shows that the preferential directions correspond to strands that stick to other obstacles (e.g. for 
$\theta =15^{\circ }$), either to the closest neighbours at 
$0^{\circ }$ and 
$60^{\circ }$ or to the second closest at 
$30^{\circ }$. The fact that three models describing slightly different physics all yield similar results demonstrates robustness across constitutive laws. An example corresponding to the case 
$\theta =25^{\circ }$ is given in supplementary movie 1 for the Oldroyd-B model.
Figure 2. Steady viscoelastic flow through a hexagonal lattice. (a) Plots of the angle 
$\alpha$ between the direction of the average flow velocity and 
$\boldsymbol {x}$ as a function of the Weissenberg number, 
$Wi$, for the viscosity ratio 
$\beta =1$, different values of 
$\theta$ and the FENE-P, FENE-CR and Oldroyd-B models. Each curve is obtained by fixing the angle 
$\theta$ and progressively increasing 
$Wi$, starting from 
$Wi=0$. At 
$Wi=0$, the flow is Newtonian so that the average flow velocity is along the direction of the forcing term and the value of 
$\theta$ for each curve can be determined as 
$\theta =\alpha (Wi=0)$. (b) Plots of the angle 
$\alpha$ as a function of 
$Wi$ for the Oldroyd-B model and different values of 
$\theta$ and 
$\beta$. In (a,b) points are actual computations and lines are just guides for the eyes; graphs are limited to the range 
$\theta \in [0,30^{\circ }]$ because of the 6-fold rotational symmetry combined with the reflection symmetry about the 
$\boldsymbol {x}$ axis; and the Weissenberg number is defined with a reference velocity different from the averaged flow velocity, so that 
$Wi$ values may seem larger than in other studies. For comparison purposes, the same figure is provided in figure 3 of the supplementary material for a Weissenberg number defined with the average flow velocity.
Figure 3. Steady viscoelastic flow through a hexagonal lattice of cylindrical obstacles. Fields of polymer stretching – the trace of the conformation tensor, 
$\textrm {tr}(\boldsymbol {c})$ – for the Oldroyd-B model with a ratio of polymer to solvent viscosities 
$\beta =1$ and different values of the Weissenberg number 
$Wi$ and angles 
$\theta$ from the aligned configuration.
We hypothesized that the mechanism that makes strands appear sticky is a feedback between the transport of the conformation tensor and the flow, which is reminiscent of the amplification mechanism of preferential flow paths in the case of side-by-side obstacles (Hopkins et al. Reference Hopkins, Haward and Shen2021) (see also Mokhtari et al. Reference Mokhtari, Latché, Quintard and Davit2022). To illustrate this effect, let us consider again the case 
$\theta =25^{\circ }$ in supplementary movie 1. Initially, strands tend to form an angle of 
$25^{\circ }$ in the wake of each circle, such that the distance with the circle immediately below is smaller than the distance with the circle above. The flow path of least resistance is above the strands with a positive feedback loop that tends to curve the streamlines towards the path of lowest flow rate and thus displace the strands further down, until eventually the path below is completely closed.
To test this hypothesis, we studied the behaviour of the strands for different values of the parameter 
$\beta$. As detailed in Mokhtari et al. (Reference Mokhtari, Latché, Quintard and Davit2022), 
$\beta$ modulates the feedback of the strands upon the flow field – it is a prefactor of the divergence of the conformation tensor and thus acts as a weight for the polymer-associated force in the momentum transport equation. For instance, the limit 
$\beta \rightarrow 0$, which may be thought of as the limit of vanishingly small polymer concentration, eliminates completely the feedback of the polymers upon the flow with a velocity field given by Stokes flow and the conformation tensor by (1.11) with a fixed Stokes velocity. The graph in figure 2(b) shows that, for 
$\beta =0.1$, the preferential flow directions almost completely disappear and strands do not stick to other obstacles. When 
$\beta =10$, however, strands are sticky for even smaller Weissenberg numbers, thus suggesting that stickiness does indeed result from a feedback between the transport of the conformation tensor and flow.
To further understand the sticky properties of the strands, we next study the variations of the flow angle 
$\alpha$ with the force angle 
$\theta$ for fixed values of 
$Wi$. Figure 4 shows in blue how 
$\alpha$ evolves when progressively increasing 
$\theta$ from 
$0$ to 
$60^{\circ }$. In the limit of small Weissenberg numbers, we observe a linear behaviour that is characteristic of the Newtonian case. For a moderate value of the Weissenberg number (
$Wi=8.4$), the evolution of the blue curve resembles a sigmoid, which is a signature of strands interacting with the closest obstacles (distance 
$S$ at 
$0^{\circ }$ and 
$60^{\circ }$). For 
$Wi=14$, the shape of the blue curve changes for 
$\theta \simeq 30^{\circ }$ with strands that are long enough to reach the second closest obstacle (distance 
$\sqrt {3}S$ at 
$30^{\circ }$ in the staggered direction), thus creating a preferential flow direction at 
$\alpha \simeq 30^{\circ }$ and breaking the linearity of the blue curve in the corresponding neighbourhood. This effect evolves non-monotonically with the Weissenberg number, with only a small perturbation in the blue curve for 
$Wi=8.4$, a clear preferential direction at 
$Wi=14$ and 
$Wi=19.6$ and a complete disappearance of this effect for 
$Wi=28$. For the largest 
$Wi$ values, this is because the interactions of the horizontal strands with the closest obstacles at 
$\alpha \simeq 0^{\circ }$ get stronger, so that ultimately detachment from the horizontal position occurs for 
$\theta >30^{\circ }$ and the strands switch directly from 
$\alpha =0^{\circ }$ to 
$\alpha =60^{\circ }$. At obstacle scale, this increase in the strength of the interaction and in the stability of the flow at 
$\alpha \simeq 0^{\circ }$ corresponds to a doubling of the strand with the formation of an envelope of stress wrapped around the obstacles – a minimal example showing that these envelopes result from the interaction of strands with downstream obstacles is provided in figure 5 and additional details can be found in Mokhtari et al. (Reference Mokhtari, Latché, Quintard and Davit2022).
Figure 4. Bifurcation diagrams with the angle 
$\theta$ of the force density showing multistability and hysteresis for the Oldroyd-B model with 
$\beta =1$. In blue, we present the results obtained when increasing 
$\theta$ from 
$0$ to 
$60^{\circ }$ and in red those obtained when decreasing 
$\theta$ from 
$60^{\circ }$ to 
$0$. To obtain the green curves, we proceeded by either using the results from a point 
$\alpha \simeq 45^{\circ }$ on the blue curve as the initial condition and progressively decreasing 
$\theta$, or from a point 
$\alpha \simeq 15^{\circ }$ on the red curve and increasing 
$\theta$. Points are actual computations and lines are just guides for the eyes: (a) 
$Wi=1.4$; (b) 
$Wi=8.4$; (c) 
$Wi=14.0$; (d) 
$Wi=19.6$; (e) 
$Wi=25.2$; and (f) 
$Wi=28.0$.
Figure 5. Illustration of the concept of the envelope. (a) Strand in the wake of a single cylinder (b) envelope of stress wrapping two cylinders aligned with the flow (
$\beta =1$, 
$Wi=200$ and the distance between the centres of the two cylinders is equal to seven times the radius of the cylinders). Colours indicate the value of 
$\textrm {tr}(\boldsymbol {c})$. The boundary conditions are left–right and top–bottom periodicity with a channel of sufficient length to avoid boundary effects – only a small part of this channel is shown here. The characteristic distance is the channel height 
$H=20$. Details of these calculations and a more thorough discussion of this problem can be found in Mokhtari et al. (Reference Mokhtari, Latché, Quintard and Davit2022).
This threshold effect in 
$\theta$ with strands strongly sticking to closest neighbours is associated with subcritical bifurcations and hysteresis. Figure 4 also shows in red the path obtained when decreasing 
$\theta$ starting from 
$\theta =60^{\circ }$. Symmetrically, the strands are now stuck to the closest neighbours at 
$60^{\circ }$ and the transition occurs for lower values of 
$\theta$ compared with the blue curve. We can also recover an intermediate path at 
$30^{\circ }$ – in green in figure 4 – by either decreasing 
$\theta$ starting from an initial condition corresponding to a point at 
$\alpha \simeq 45^{\circ }$ on the blue curve or increasing 
$\theta$ from a point at 
$\alpha \simeq 15^{\circ }$ on the red curve. This green path corresponds to strands sticking to the second closest neighbours but can only be obtained when the strands are not initially stuck to the closest neighbours. We further observe two hysteresis regimes with the Weissenberg number. For 
$Wi=19.6$, there are two separate hysteretic loops featuring bistability around 
$\alpha =0^{\circ }$ and 
$30^{\circ }$ or around 
$\alpha =30^{\circ }$ and 
$60^{\circ }$. These loops then meet to yield tristability around 
$\alpha =0^{\circ }$, 
$30^{\circ }$ and 
$60^{\circ }$ (see example for 
$Wi=24$ in supplementary movie 2 for the Oldroyd-B model). The 
$(\theta,Wi)$ phase diagram in figure 6 summarizes the different regimes and shows insets of the corresponding stretching field.
Figure 6. The 
$(\theta,Wi)$ phase diagram with the different stability domains for the Oldroyd-B model with 
$\beta =1$. Insets show the corresponding polymer stretching fields and the position of the strands; (
$\circ$) indicates a single stable state, (
$\blacktriangle$, blue) bistability and (
$\blacklozenge$, red) tristability. Stripes indicate a transition towards unsteady flows.
2.2. Destabilization and unsteady flows
Upon further increasing the Weissenberg number, we observe a spontaneous transition to unsteady flows with strands that initially form along the direction of the force and then destabilize. To study this transition, we focus on the FENE-CR model to avoid complications associated with shear thinning (vs FENE-P) (Herrchen & Öttinger Reference Herrchen and Öttinger1997) and to maximize accuracy (vs Oldroyd-B) (Kim et al. Reference Kim, Kim, Chung, Ahn and Lee2005; Mokhtari et al. Reference Mokhtari, Davit, Latché and Quintard2023). Since multistability occurs for the staggered direction 
$\theta ={30}^{\circ }$, we expect the most interesting dynamics to develop for this case and therefore consider it first – as we will see, the multistability previously described in the steady state is important to understand the dynamics. We proceed by placing ourselves in the frame of reference 
$(\boldsymbol {x}^{\star },\boldsymbol {y}^{\star })$ (see definition in figure 1), rather than changing the angle of the forcing term in the aligned direction 
$\theta ={0}^{\circ }$. Since we use a structured mesh, this approach makes it easier to match the symmetries of the mesh with that of the problem – closest neighbours are at 
${-30}^{\circ }$ and 
${+30}^{\circ }$ angles in 
$(\boldsymbol {x}^{\star },\boldsymbol {y}^{\star })$, rather than 
${0}^{\circ }$ and 
${60}^{\circ }$ in 
$(\boldsymbol {x},\boldsymbol {y})$.
Based upon our simulations, we propose a classification of the different regimes as
(I) Self-sustained flapping oscillations about base flow.
(II) Travelling waves with fronts separating zones with different configurations of the strands.
(III) An intermediate regime with remaining flapping oscillations combined with large amplitude/low frequency fluctuations.
(IV) A quasi-steady regime with no oscillations and large amplitude/low frequency fluctuations leading to a steady state.
(V) Strands pulsations and fluctuations over a range of frequencies.
Regime (I) appears for relatively low values of 
$\beta$ and 
$Wi$, at the onset of the transition to unsteady flow. Oscillations develop about the direction of the force density with what appears to be a Hopf bifurcation (see supplementary movie 3). Figure 7(a) shows that the average position of the strands, as visible on the field of the time average of the trace of the conformation tensor, 
$\overline {\textrm {tr}(\boldsymbol {c})}$, is horizontal with r.m.s. variations about this position. The time-averaged flow in the wake of each cylinder is similar to the steady case, with strands essentially acting as a tangential force opposed to the flow (Mokhtari et al. Reference Mokhtari, Latché, Quintard and Davit2022). Figures 7(b)–7(e) further show out-of-phase flapping oscillations from one row to the next.
Figure 7. Regime (I): self-sustained flapping oscillations, obtained for the FENE-CR model with 
$b=1000$ and 
$(\beta =0.5, Wi=47.6)$. (a) The first and second columns show the time-averaged field of the trace of the conformation tensor and the corresponding r.m.s. variations, respectively, for the time interval 
$(7500\leqslant t\leqslant 9500)$. The third and fourth columns show the normalized time-averaged field of the velocity and the corresponding r.m.s. variations, respectively. (b) Definition of the points of interest 
$P_{1}$ (blue) and 
$P_{2}$ (red) in the staggered configuration. Here, 
$\boldsymbol {f}$ shows the direction of the volume force. (c) Time variations of the local velocity angle 
$\alpha _{P}$ at 
$P_{1}$ and 
$P_{2}$. (d) Phase diagrams with 
$\alpha _{P}$ as a function of the local trace of the conformation tensor at 
$P_{1}$ and 
$P_{2}$. (e) Kymographs of 
$\textrm {tr}(\boldsymbol {c})$ averaged over the vertical direction and normalized by its time average.
In regime (II), strands are still relatively short but start interacting more strongly with closest neighbours at 
$\pm 30{}^{\circ }$. Time-averaged fields for the trace of the conformation tensor, 
$\overline {\textrm {tr}(\boldsymbol {c})}$, in figure 8(a) now show that strands tend to stick to closest neighbours and periodically switch positions. Orbits in figure 8(c,d) confirm that strands do not oscillate about the base flow configuration but rather that they stick to one side or the other, while generating a travelling wave solution in the direction of the flow (kymograph in figure 8(e) and supplementary movie 4). To understand this phenomenon, consider that each strand may be roughly thought of as a barrier allowing or preventing flow through a given pore throat. These barriers have a limited number of accessible positions. They are either stuck to one or the other neighbour and, in so doing, allow or prevent flow. Globally, strands also arrange in a way that maintains flow through the structure, this leading to a set of possible configurations. Figure 8(f) presents such solutions with strands successively sticking up and down so as to form tortuous ‘zigzag’ channels in the 
$x$ direction. Figure 8(f) further shows that two topologically equivalent configurations are possible for the same type of zigzag channels. Some of the observed travelling wave solutions consist of vertically invariant zones of one of these two configurations separated by a front where strands switch side (figure 8(g) and supplementary movie 4).
Figure 8. Regime (II): travelling waves for the FENE-CR model with 
$b=1000$ and 
$(\beta =2, Wi=47.6)$. (a) The first and second columns show the time-averaged field of the trace of the conformation tensor and the corresponding r.m.s. variations, respectively, for the time interval 
$(7500\leqslant t\leqslant 9500)$. The third and fourth columns show the normalized time-averaged field of the velocity and the corresponding r.m.s. variations, respectively. (b) Definition of the points of interest 
$P_{1}$ (blue) and 
$P_{2}$ (red) in the staggered configuration. Here, 
$\boldsymbol {f}$ shows the direction of the volume force. (c) Time variations of the local velocity angle 
$\alpha _{P}$ at 
$P_{1}$ and 
$P_{2}$. (d) Phase diagrams with 
$\alpha _{P}$ as a function of the local trace of the conformation tensor at 
$P_{1}$ and 
$P_{2}$. (e) Kymographs of 
$\textrm {tr}(\boldsymbol {c})$ averaged over the vertical direction and normalized by its time average. (f) Illustration of the two equivalent strand configurations with a travelling front (dashed line) corresponding to strands switching side. (g) Results of the simulations for 
$\textrm {tr}(\boldsymbol {c})$ for the case 
$(\beta =2, Wi=47.6)$ showing the position of the front aligned with the dashed line in (f).
In regimes (III), (IV) and (V), strands get longer and interact more strongly with other obstacles, so that initial flapping oscillations tend to disappear. We observe a stabilization of the zigzag channels (figure 9b,e) and, for regimes (IV) and (V), the appearance of envelopes of localized stress along the 
$\pm 30^{\circ }$ angle, whereby the cylinders are encased between two strands of high polymeric stress instead of each one generating a strand in their wake. An example from regime (IV) in figure 9(c,f) for the case 
$(\beta =5, Wi=19.6)$ shows a steady-state mix of zigzag channels and envelopes. Figure 10 shows that envelopes tend to be more prominent as we increase 
$Wi$ and 
$\beta$ in regime (V), ranging from only a few instances for 
$(\beta =5, Wi=42)$ to covering the entire lattice for 
$(\beta =10, Wi=58.8)$ with a pattern of time-averaged localized stress that is similar to that of the aligned case at steady state in figure 3. Contrary to the aligned case, however, this pattern now implies that there is a misalignment between the flow and the imposed force, with strands directing the flow in the (
$\pm 30^{\circ }$) directions, as visible from the velocity fields in figure 10. We hypothesized that this phenomenon is the consequence of the stability of the envelope pattern at 
$\pm 30{}^{\circ }$. To test this hypothesis, we performed calculations in the aligned case and show in figure 4 of the supplementary material that this pattern is associated with a much simpler dynamics, where envelopes remain stable over a wide range of values of 
$\beta$ and 
$Wi$.
Figure 9. Examples of steady-state fields for 
$\mathrm {tr}(\boldsymbol {c})$ and the corresponding normalized amplitude of the velocity; (a) (
$\beta =0.5, Wi=42$) corresponds to the trivial steady solution, (b) (
$\beta =1, Wi=53.2$) to the regime (III) and (c) (
$\beta =5, Wi=19.6$) to regime (IV).
Figure 10. Regime (V): strand pulsations for the FENE-CR model with 
$b=1000$. Each line corresponds to a set of dimensionless numbers with results presented for different time intervals. The first and second columns show the time-averaged fields of the trace of the conformation tensor and the corresponding r.m.s. variations, respectively. The third and fourth columns show the normalized time-averaged fields of the velocity and the corresponding r.m.s. variations, respectively.
In these regimes, we also observe defects at the junction of zones with different stress configurations, for example between envelopes and zigzags as visible on the top right of figure 10(b) for 
$\textrm {tr}(\boldsymbol {c})$. These defects are sometimes displaced through the lattice of obstacles, as shown in figure 10(a), or are stable in a given position, as shown in figure 10(b). The propagation of a defect through the structure occurs in a variety of ways, ranging from slow propagation from neighbour to neighbour to large puffs of rapid strand movements, as illustrated in figure 10(a). In many instances, the defects are only ephemerous and ultimately disappear (supplementary movie 6a and b), leading to a stable pattern of strands. In regime (V), we observe beatings/pulsations that develop along the strands (supplementary movies 5 and 6a and b) with a stable pattern of stress and no change in the network of flow channels. Examples of such pulsations are presented in fields of 
$\textrm {tr}(\boldsymbol {c})_{rms}$ in figure 10 at both 
$\beta =5$ and 
$\beta =10$.
Figure 11 summarizes the different regimes in a 
$(Wi, \beta )$ diagram and illustrates the corresponding variations in time of the global flow angle 
$\alpha$. Regime (I) at 
$(\beta =0.5, Wi=47.6)$ features distinct oscillations of 
$\alpha$ that correspond to the periodic flapping motion of the birefringent strands – the power spectrum has a clear dominant frequency (see figure 12). In regime (II), we also have clear oscillations with strands that successively stick to closest neighbours and a distinct frequency in the power spectrum (see figure 12). For regimes (III), (IV) and (V), strands stick strongly to obstacles and the initial flapping oscillations tend to disappear, with few remaining oscillations in regime (III) (figure 11) and only low frequency/large amplitude temporal fluctuations in regimes (IV) and (V) (figure 11). For regime (V), we observe pulsations at later times that are very different in nature from the self-sustained oscillations in regime (I). They feature a more complex temporal response and a power spectrum (see figure 12) that has no distinct frequency, but rather shows some of the hallmarks of chaos and elastic turbulence (Groisman & Steinberg Reference Groisman and Steinberg2000; Steinberg Reference Steinberg2021). Calculations in the aligned case, as presented in the supplementary material figure 4, also show pulsations with a similar spectral response (see figure 12). This regime is the only unsteady bifurcation that we were able to observe in the aligned case, thus suggesting that strand beating pulsations are integral to the transition towards elastic turbulence.
Figure 11. Diagram of the different regimes obtained for the FENE-CR with 
$b=1000$ in the staggered configuration. The main graph shows the different regimes as a function of 
$\beta$ and 
$Wi$. The other plots are examples of time evolutions of 
$\alpha$ for each regime. The colour code indicates whether the strands remain short (blue), without envelope formation or become long and form envelopes (red). Here, (
${\times }$, blue) corresponds to the trivial steady solution (figure 9(a,d), (
${\triangle }$, blue) to oscillations in regime (I), (
${\triangledown }$, blue) to travelling fronts in regime (II), (
${\blacktriangle }$, blue) to the intermediate in regime (III), (
${\square }$, red) to the quasi-steady regime (IV) and (
${\blacksquare }$, red) to the pulsations in regime (V).
Figure 12. Power spectral density associated with the temporal fluctuation of 
$\alpha$ for different values of the parameters 
$\beta$ and 
$Wi$. The figure on the left-hand side corresponds to regimes (I) and (II) in the staggered configuration, the one at the centre to regime (V) in the staggered configuration and the one on the right-hand side to regime (V) in the aligned configuration. The dimensionless frequency 
$f$ corresponds to the non-dimensionalization presented in the article (Methods). We also show the same figure in the supplementary material figure 6 for a frequency corresponding to time non-dimensionalized with the characteristic time for polymer relaxation.
3. Discussion
Our simulations do not represent a perfect description of experiments, such as those performed in Walkama et al. (Reference Walkama, Waisbord and Guasto2020) or Haward et al. (Reference Haward, Hopkins and Shen2021b). Due to computational limitations, we have considered a 2-D geometry based on a relatively small representative portion of the porous medium and may have thus neglected important 3-D effects. Further, we have used periodic conditions with a constant forcing term, whereas unstable flow may generate spatio-temporal variations of the locally averaged pressure gradient in real large pore structures. Oldroyd and FENE models are also imperfect; for instance, by assuming some form of uniformity of the concentration of polymers in the porous medium. The question that we ask in this discussion section is whether, despite these limitations, we can extract useful information from our simulations to understand viscoelastic instability mechanisms in porous media. To this end, we first discuss links between our simulations and experimental results from the literature.
Sequences of transitions from self-sustained oscillations to aperiodic fluctuations have been observed in different systems involving viscoelastic fluids, such as axisymmetric contractions (McKinley et al. Reference McKinley, Raiford, Brown and Armstrong1991) and microfluidic cylinders (Haward et al. Reference Haward, Kitajima, Toda-Peters, Takahashi and Shen2019). Hopf bifurcations have also been found to be linked to birefringent strands in cross-slot flows (Xi & Graham Reference Xi and Graham2009). Recent experiments have specifically studied the dynamics of viscoelastic flows through a staggered hexagonal lattice and point towards an excellent agreement with our simulations:
(i) Self-sustained oscillations are evident in supplementary movie S1 in Haward et al. (Reference Haward, Hopkins and Shen2021b), clearly showing strands flapping from one side to the other in a way that is directly comparable to our simulations in supplementary movie 3. Similarly, the retardation field in figure 3 in Haward et al. (Reference Haward, Hopkins and Shen2021b) for the lowest Weissenberg numbers shares similarities with
$\overline {\textrm {tr}(\boldsymbol {c})}$ and $\textrm {tr}(\boldsymbol {c})_{rms}$ in figures 7–8.(ii) A transition from short oscillating strands to zigzag channels is visible in figure 12D in Datta et al. (Reference Datta, Ardekani, Arratia, Beris, Bischofberger, McKinley, Eggers, López-Aguilar, Fielding and Frishman2022), when increasing the Weissenberg number in the staggered configuration. Supplementary movie S1 in Haward et al. (Reference Haward, Hopkins and Shen2021b) also shows the temporary formation of structures that are similar to our zigzag channels.
(iii) Supplementary movies S2 and S3 in Walkama et al. (Reference Walkama, Waisbord and Guasto2020) show travelling waves with fronts that delineate zones with different strand configurations.
(iv) Figure 3 and supplementary movie S3 in Haward et al. (Reference Haward, Hopkins and Shen2021b) show two strands per cylinder in the aligned case and the formation of an envelope of polymer stress.
(v) Supplementary movies S2 and S4 in Haward et al. (Reference Haward, Hopkins and Shen2021b) show stable patterns of retardation that seem to exhibit some of the pulsations observed in our simulations (supplementary movie 5).
de Blois, Haward & Shen (Reference de Blois, Haward and Shen2023) and Ström, Beech & Tegenfeldt (Reference Ström, Beech and Tegenfeldt2023) have also demonstrated the existence of waves in polymer flows through a square lattice of obstacles. de Blois et al. (Reference de Blois, Haward and Shen2023) show the emergence of elastic waves in a microfluidic canopy of either rigid or flexible slender pillars. Although the geometry and the macroscopic boundary conditions are different from our simulations, the development of large-scale coherent structures with a deflection of the flow direction could indicate that sticky birefringent strands are involved. Ström et al. (Reference Ström, Beech and Tegenfeldt2023) also show waves in a microfluidic square lattice, but with important differences in the experimental approach. Ström et al. (Reference Ström, Beech and Tegenfeldt2023) use DNA with a size of the polymer that is similar to that of the gap between obstacles, whereas de Blois et al. (Reference de Blois, Haward and Shen2023) use hyaluronic acid with pores that are much larger than the polymers. The observed waves in Ström et al. (Reference Ström, Beech and Tegenfeldt2023) correspond to polymer concentration and conformation, whereas de Blois et al. (Reference de Blois, Haward and Shen2023) show waves in the velocity vector field and the deflection of flexible pillars. The geometry is also slightly different, in particular the positions of the inlet/outlets and the gap between the top of the pillars and the device cover in de Blois et al. (Reference de Blois, Haward and Shen2023). Notwithstanding these technical differences, Ström et al. (Reference Ström, Beech and Tegenfeldt2023) also show that there is a transition from steady depletion to cyclic fluctuations and then waves; that the waves correspond to areas of stretched polymers; that the concentration of DNA – similar to changes in 
$\beta$ in our model – modulates the instability with no waves for sufficiently low concentrations; that there is a depletion zone between pillars and the formation of vortices; that a slight randomization of pillar position eliminates waves. These observations suggest that the mechanisms at play could in fact be similar to that described in Walkama et al. (Reference Walkama, Waisbord and Guasto2020), Haward et al. (Reference Haward, Hopkins and Shen2021b), Mokhtari et al. (Reference Mokhtari, Latché, Quintard and Davit2022) and the present work. Future simulations in square lattices should clarify similarities and differences.
All of our simulations are two-dimensional. One can argue that the nature of viscoelastic instabilities arising in microfluidic systems is in essence three-dimensional. In studying the elastic instability upstream of a single cylinder – height 
$60\ \mathrm {\mu } \textrm {m}$, diameter 
$50\ \mathrm {\mu } \textrm {m}$ and channel width 
$100\ \mathrm {\mu } \textrm {m}$ – using 3-D holographic particle velocimetry, Qin et al. (Reference Qin, Salipante, Hudson and Arratia2019) have shown that upstream vortices initiate at the corners between the cylinder and the wall, so that these corners actually play a central role in the instability (Poole Reference Poole2019). Three-dimensional flows may also develop independently from the wall effect, directly due to instabilities (Gutierrez-Castillo, Kagel & Thomases Reference Gutierrez-Castillo, Kagel and Thomases2020). And it is difficult to think that, at very large Weissenberg numbers, fully developed elastic turbulence in porous media can remain two-dimensional, even in microfluidic systems with slender obstacles. Three-dimensional unstable simulations would therefore be extremely valuable. On the other hand, even if we could overcome the tremendous computational cost of such 3-D computations, the substantial agreement between our simulations and microfluidic experiments in Walkama et al. (Reference Walkama, Waisbord and Guasto2020) and Haward et al. (Reference Haward, Hopkins and Shen2021b) suggests that 2-D simulations capture some important physics. It also suggests that 2-D instabilities are a fundamental part of the problem for the slender microfluidic geometry and the transitions observed in Walkama et al. (Reference Walkama, Waisbord and Guasto2020) and Haward et al. (Reference Haward, Hopkins and Shen2021b).
Details in the geometry of the system and in the rheology of the fluid may further modulate the relative weight of 3-D effects and modify regimes and transitions. For instance, an obvious difference between experiments in Qin et al. (Reference Qin, Salipante, Hudson and Arratia2019) and Haward et al. (Reference Haward, Hopkins and Shen2021b) is the ratio of height to diameter, which is close to 1 for Qin et al. (Reference Qin, Salipante, Hudson and Arratia2019) and 10 in Haward et al. (Reference Haward, Hopkins and Shen2021b). It is possible that vortices contribute to the dynamics in both cases, but that the relative influence of the 2-D geometry and of the corners in controlling the stability strongly depends upon the ratio of height to diameter. Another difference between these experiments is that Qin et al. (Reference Qin, Salipante, Hudson and Arratia2019) considers a single cylinder whereas Haward et al. (Reference Haward, Hopkins and Shen2021b) deals with an array of cylinders. The instability mechanism in the case with one cylinder involves the growth of the vortex length, which can reach several times the cylinder diameter. In inertial turbulence through porous media, the size of turbulent structures is restricted by the size of the pores (Jin et al. Reference Jin, Uth, Kuznetsov and Herwig2015). It is possible that a similar constraint exists for viscoelastic flows with upstream obstacles limiting the size of the vortex and thus the impact of the walls upon the dynamics. Ribeiro et al. (Reference Ribeiro, Coelho, Pinho and Alves2014) have shown that, for the case of a single confined cylinder, shear-thinning effects have an impact on the development of elastic instabilities for a variety of aspect ratios.
What then did we learn from 2-D simulations? First, our interpretation of spatio-temporal fluctuations in terms of the stability of a web of sticky strands provides an intuitive basis for interpreting experimental results. It clarifies the link between the ‘screening’ of the stagnation points described in Haward et al. (Reference Haward, Hopkins and Shen2021b) and the development of instabilities. Haward et al. (Reference Haward, Hopkins and Shen2021b) hypothesized that stagnation points are screened from the flow in the aligned geometry – resting their argument upon example calculations of creeping flows through the different structures. They suggest that screened stagnation points lose their ‘strength’ as a source of polymer stress, inducing a delay in the occurrence of the instability. We have previously shown in Mokhtari et al. (Reference Mokhtari, Latché, Quintard and Davit2022) that this phenomenon is not only due to the geometry in Stokes flow but also stems from a feedback between flow and birefringent strands for viscoelastic flows. When strands form an envelope around cylinders, they generate a stagnant zone between cylinders that essentially behaves like a two-sided lid-driven cavity. In the present work, we further provide a connection between the patterns of localized stress and the stability: (i) instabilities develop more easily in the absence of envelope patterns; (ii) the initial instability in the staggered case corresponds to a destabilization of individual strands in the wake of each cylinder with an initially tristable configuration at steady state and what looks like a Hopf bifurcation; (iii) the transition towards elastic turbulence corresponds to a destabilization of the envelope pattern in all our simulations – pulsations regime (V). We therefore argue that what controls the stability of the flow is the pattern of localized stress, in particular the very stable envelope pattern.
Our results also show that, in the staggered configuration, there is first a transition to unsteady flow, then a restabilization and finally a transition towards elastic turbulence. We initially observe the formation of a single strand in the wake of each cylinder, so that the first transition is likely to fit within the classical framework of viscoelastic flow past a single cylinder. The elastic instability in flows with curved streamlines is widely considered as being driven by hoop stress (Shaqfeh Reference Shaqfeh1996). The Pakdel–McKinley criterion (McKinley, Pakdel & Öztekin Reference McKinley, Pakdel and Öztekin1996; Pakdel & McKinley Reference Pakdel and McKinley1996) quantifies this using a dimensionless parameter, 
$M$, that takes into account the relative importance of the tensile stress and the ratio of the relaxation length of a perturbation to the radius of curvature of streamlines. We thus hypothesize that the onset of the instability (regime (I)) may be treated in a way similar to § 4.2 in McKinley et al. (Reference McKinley, Pakdel and Öztekin1996). The restabilization observed in regime (IV) is more puzzling. Our hypothesis is that, as the strands grow longer, they modify the stress pattern – a mixture of envelopes and zigzag channels – leading to lower values of the tensile stress and reduced curvature of the streamlines – a phenomenon that is linked to the screening of stagnation points, as discussed in Haward et al. (Reference Haward, Hopkins and Shen2021b). Thus, it may be possible to define a second critical value of the dimensionless number 
$M$ for the onset of the regime with pulsations of the strands forming envelopes.
This picture is consistent with the non-monotonic variation of the spatially averaged retardation fluctuations with the Weissenberg number observed in Haward et al. (Reference Haward, Hopkins and Shen2021b). Our simulations suggest that the first increase is due to self-sustained oscillations and waves (regimes I and II), that the following decrease results from the stabilization of long strands (regimes III and IV) and that the second increase is caused by the destabilization of the envelope patterns and the appearance of strand pulsations (regime V). This is also consistent with the observation in Haward et al. (Reference Haward, Hopkins and Shen2021b) that the staggered and aligned configurations behave similarly at sufficiently high Weissenberg numbers: the instability in both cases corresponds to a destabilization of the envelopes. In fact, this may also suggest that the transition to elastic turbulence and chaos – not just unsteady flow – actually develops for similar Weissenberg numbers in both configurations (see figure 5 in the supplementary material), thus shining a different light on transitions in Walkama et al. (Reference Walkama, Waisbord and Guasto2020).
Another interesting observation in our simulations is what happens as we get away from transitions. We show that the problem rapidly becomes strongly nonlinear and non-local. The mechanism that drives self-sustained oscillations and waves (regimes I and II) is not only a combination of streamline curvature and large tensile stress, but also stems from interactions between the strands and obstacles and from a coupling between flow channels through the entire structure. Similarly, once pulsations develop, fluctuations in one pore modify the behaviour of neighbouring pores – an idea found in recent studies (Browne, Shih & Datta Reference Browne, Shih and Datta2020a; Kumar et al. Reference Kumar, Aramideh, Browne, Datta and Ardekani2021) – emphasizing the importance of spatial correlation, non-locality and perturbation propagation in viscoelastic flows through porous media.
4. Conclusions
At steady state, we have shown that interactions between birefringent strands and obstacles guide the flow through the 2-D structure and control the direction of the average flow. Strands tend to stick to obstacles and, in so doing, create preferential flow directions through the lattice. This phenomenon leads to a succession of subcritical bifurcations with the angle of the forcing term, which yields multistability and hysteresis. Upon further increasing the Weissenberg number 
$Wi$ and the ratio of viscosities 
$\beta$, we discovered a rich zoology of spatio-temporal fluctuations, including self-sustained flapping oscillations, travelling wave solutions, defect propagation and strand pulsation fluctuations – this last regime being integral to the transition towards elastic turbulence. Our work demonstrates that strands and their interactions with flow and obstacles control the steady base flow, the spontaneous transition to unsteady flow and the route to chaos in 2-D lattices of obstacles. More generally, this suggests that the stability of localized stress patterns is key in understanding the transition towards elastic turbulence in complex structures.
Supplementary material and movies
Supplementary material and movies are available at https://doi.org/10.1017/jfm.2023.916.
Funding
The authors thank R. de Loubens, from TotalEnergies, for valuable discussions on this topic and F. Babik, from the 
$\mathrm {CALIF}^{3}\mathrm {S}$ development team at IRSN, for his precious support. The authors would also like to thank TotalEnergies for funding and supporting this work, in particular through access to the PANGEA II supercomputer.
Declaration of interests
The authors report no conflict of interest.
Open access
