L4 · Chapter 7

How a four-roll flow loses symmetry at its stagnation point, and how much the numerics change it

Sandbox Physics · 2026-10-29

Push an inertia-free polymer solution with a steady, symmetric force. Can the flow break the symmetry by itself and start to oscillate? The four-roll model of Thomases and Shelley gives a numerical yes: above a Weissenberg number of about 5 the stress strand at the stagnation point tilts to one side, and further up the flow oscillates persistently and mixes. A later study argues that the symmetry breaking comes, at least in part, from the stress diffusion added to the computation. The course sets out both sides and gives a way to judge whether a numerical result can be trusted.

Four vortices and one stagnation point

The domain is a periodic square box of side $2\pi$ . A steady force $\mathbf f$ drives four vortices around a central stagnation point, where fluid arrives along one axis and leaves along the other, as in the planar extension of Chapter 3. Without polymers the velocity is exactly half the force.

The time unit $\tau_f$ is the inverse strain rate at the stagnation point of the Newtonian flow, and the Weissenberg number is the relaxation time divided by it. The polymer–solvent coupling $\beta$ is adjusted with the Weissenberg number so that $\beta\mathrm{Wi}=1/2$ stays fixed, which fixes the ratio of polymer to solvent viscosity at 0.5 to compare with the cross-slot experiments of Arratia and co-workers. The stress equation also carries a diffusion term $\nu_p\Delta\mathbf S$ , discussed below.

-\nabla p+\Delta\mathbf u=-\beta\nabla\cdot\mathbf S+\mathbf f,\quad \nabla\cdot\mathbf u=0,\qquad \overset{\nabla}{\mathbf S}+\mathrm{Wi}^{-1}\left(\mathbf S-\mathbf I\right)=\nu_p\Delta\mathbf S,\qquad \mathbf f=\left(2\sin x\cos y,\,-2\cos x\sin y\right)

Eqs. (1), (2) and (6) of Thomases, Shelley and Thiffeault; the conformation tensor is written S there.

Source: B. Thomases · B. Thomases and M. Shelley

The stress strand at the stagnation point

Near the stagnation point the velocity is close to $\alpha\,(x,-y)$ , so the polymers feel an effective Weissenberg number $\alpha\mathrm{Wi}$ . The polymer stress in turn slows the flow at the stagnation point: as the Weissenberg number goes to zero $\alpha$ tends to 2/3, and it falls further as the Weissenberg number grows.

Thomases and Shelley (2007) computed this model without stress diffusion on grids of up to 4096 points per side. When the effective value reaches 1/3 the stress profile along the outgoing direction develops a cusp, at $\mathrm{Wi}\approx0.5$ ; when it reaches 1/2 the stress grows without bound, at $\mathrm{Wi}\approx0.9$ . The latter is exactly the planar-extension threshold of Chapter 3.

The region of stretched polymers therefore contracts into a strand along the outgoing direction, ever narrower and stronger. With FENE-P the stress stays bounded, but the cusp-like structures remain.

\mathbf u\approx\alpha\,(x,\,-y),\qquad \mathrm{Wi}_{\mathrm{eff}}=\alpha\,\mathrm{Wi}

The local linear flow near the stagnation point, as written by Thomases and Shelley (2007).

Source: B. Thomases and M. Shelley

Symmetry breaking and oscillation

With stress diffusion $\nu_p=10^{-3}$ , Thomases, Shelley and Thiffeault first found the symmetric steady solution and then added a small perturbation: at $\mathrm{Wi}=4.8$ it decays, at $4.9$ it grows (Fig. 2). They used a $256\times256$ grid and report that doubling the resolution leaves the dynamics unchanged.

Between about 4.9 and 9 the asymmetric flow oscillates for a while and then settles into an asymmetric steady state. Above about 9 one vortex centre comes unpinned, the flow oscillates persistently, and fluid from the four quadrants begins to mix. At $\mathrm{Wi}=10$ the velocity spectrum has two periods of about 330 and 55 $\tau_f$ (Fig. 6(d)); the 2009 letter reports 317 and 59 and a symmetry breaking near $t\approx200$ . At Weissenberg numbers of 12 and 15 the motion is periodic, and above about 20 it is aperiodic.

Both papers check their computations with the balance of polymer strain energy: the energy $E$ is the integral of $\tfrac12\operatorname{tr}(\mathbf S-\mathbf I)$ , and its change is set by the work of the force, viscous dissipation and relaxation.

\dot E+\mathrm{Wi}^{-1}E=-\mathcal D+\mathcal W,\qquad E=\frac12\int\operatorname{tr}\left(\mathbf S-\mathbf I\right)\,d\mathbf x

Eqs. (4)–(5) of Thomases, Shelley and Thiffeault; see the paper for the exact dissipation and work terms.

Source: B. Thomases · B. Thomases and M. Shelley

The question raised by stress diffusion

Stress diffusion keeps the equations smooth at the strand, at the cost of changing the physics. Gupta and Vincenzi point out that $\nu_p=10^{-3}$ corresponds to a Schmidt number of about $10^3$ , three orders of magnitude below that of real polymer solutions. They ran a comparison with a cellular force similar to the four-roll mill at a Deborah number of 10. With that much diffusion the flow stays almost frozen for a long time and then breaks the cellular symmetry; with a positivity-preserving scheme that needs no added diffusion, the flow becomes chaotic quickly but keeps the symmetry of the force.

They conclude that in their setup the symmetry breaking is a spurious effect of the diffusion, and that forcings mixing strain and vorticity may be affected less. Thomases, Shelley and Thiffeault calibrate the diffusion coefficient against FENE-P in their Appendix B but do not report how the threshold changes with it.

The two sets of computations differ in the form of the force and in the numerical scheme, so it cannot yet be said how much of the four-roll symmetry breaking is physical. The question is open. The review by Datta and co-authors adds that Oldroyd-B is almost never quantitative at high Weissenberg number.

Source: A. Gupta and D. Vincenzi · B. Thomases · S. S. Datta et al.

How to check a numerical result

A numerical threshold or period should pass several independent checks before it counts as a property of the model. Refine the grid by a factor of two and the result should not change. Halve the time step and it should not change. Change the stress diffusion coefficient, at least to one smaller value, and if possible switch to a scheme that needs no artificial diffusion. Monitor the smallest eigenvalue of the conformation tensor throughout; it must stay positive, which log-conformation schemes guarantee automatically. Finally, switch to a bounded model such as FENE-P and see whether the conclusion survives.

One more check is cheap: accumulate the mismatch between the two sides of the energy balance. If it does not grow in time, the time integration is not quietly adding or removing energy. The four-roll experiment that accompanies this course shows these quantities on the page and records grid size, time step, diffusion coefficient and perturbation seed with every run. Its discretisation is coarser than the papers’, which is exactly why these checks are there.

Source: S. S. Datta et al. · Sandbox Physics FR26

Worked example: the papers’ numbers in relaxation times

Time is measured in $\tau_f$ , and the relaxation time is $\mathrm{Wi}\,\tau_f$ . At $\mathrm{Wi}=10$ the periods of 330 and 55 are 33 and 5.5 relaxation times, and the symmetry-breaking time of 200 reported in 2009 is about 20 relaxation times. The slow period is far longer than the memory of a single polymer, so it is a collective behaviour of the whole flow.

In one relaxation time, diffusion smooths over a distance of about $\sqrt{\nu_p\,\mathrm{Wi}\,\tau_f}$ . At $\nu_p=10^{-3}$ and $\mathrm{Wi}=10$ that is $\sqrt{10^{-2}}=0.1$ , against a box side of $2\pi\approx6.3$ . Above a Weissenberg number of 10 the paper uses $\nu_p=0.01/\mathrm{Wi}$ , which fixes the product $\nu_p\mathrm{Wi}$ at 0.01, so the distance stays 0.1.

The precision of the threshold: decay at 4.8 and growth at 4.9 place it at $4.85\pm0.05$ , about 1%. That holds for one grid and one diffusion coefficient; the checks of the previous section are what it takes to call it the threshold of the model.

Source: B. Thomases · B. Thomases and M. Shelley

Try it in the Lab

Find the symmetry breaking, then test it

Choose the Symmetric stress strand preset (Weissenberg number 4), then Symmetry breaks, then settles (Weissenberg number 6). Watch the asymmetry chart: does one decay towards zero while the other grows along a straight line on its log axis?
Switch to a Weissenberg number of 10, record the oscillation period of the probe velocity in units of $\tau_f$ , and compare it with the papers’ values of about 330 and 55. Write down the grid size and the stress diffusion coefficient.
Keep the Weissenberg number at 10 and change one thing at a time: a finer grid, half the time step, more stress diffusion. Which change alters the result? Then watch the smallest stretch eigenvalue: started from rest on the default 64 by 64 grid it dips below zero while the strand forms, a sign that the grid does not resolve the strand; on the 128 by 128 grid it stays positive.

Compare with the threshold of uniform planar extension

Run Oldroyd-B at Weissenberg numbers of 0.45 and 0.55 and confirm that the stretch in extension grows without bound only in the second case. That is what happens at the four-roll stagnation point once the effective Weissenberg number reaches $1/2$ .

Exercises with explained answers

If $\alpha$ at the four-roll stagnation point stayed at 2/3, at what Weissenberg number would the effective value reach 1/2? Thomases and Shelley found about 0.9. Why is it larger?
Show the answer
$\mathrm{Wi}=(1/2)/(2/3)=0.75$ . The actual value is larger because the polymer stress slows the stagnation-point flow further as the Weissenberg number grows, so $\alpha$ drops below 2/3 and a larger Weissenberg number is needed for the product to reach 1/2.
With $\nu_p=0.01/\mathrm{Wi}$ at a Weissenberg number of 20, how far does stress diffuse in one relaxation time? And over a slow period of 330, taking the period found at 10?
Show the answer
In one relaxation time $\sqrt{(0.01/20)\times20}=0.1$ . Over a period of 330, $\sqrt{5\times10^{-4}\times330}=0.41$ , already about 7% of the box side.
What evidence would you need to convince someone that the symmetry breaking belongs to the model rather than to the stress diffusion?
Show the answer
The same conclusion after refining the grid and halving the time step; the same behaviour at a smaller diffusion coefficient, or in a scheme that needs no artificial diffusion; a conformation tensor that stays positive definite throughout; and ideally the same result with FENE-P. A threshold found at one diffusion coefficient is a result of that discretisation only.
Why do Thomases and co-workers hold the product of the coupling and the Weissenberg number at one half instead of fixing the coupling?
Show the answer
That keeps the ratio of polymer to solvent viscosity at 0.5, so changing the Weissenberg number changes only the elasticity, not the viscous make-up of the solution. It matches the cross-slot experiments of Arratia and co-workers. Thomases and Shelley explain that for a given working fluid this ratio does not change with the experimental conditions.

Sources and what we recomputed

S. S. Datta et al., Perspectives on viscoelastic flow instabilities and elastic turbulence, arXiv:2108.09841v1 (2021); Physical Review Fluids 7, 080701 (2022)
arXiv version 1, 22 August 2021; the journal version was not checked. Section II, Eq. (1): the Pakdel–McKinley criterion, which combines streamline curvature with the ratio of normal stress to shear stress. No critical value is given for Taylor–Couette or cone-and-plate flow. Section III: simulations of cross-slot flows lose symmetry or become time dependent when the largest local value of that criterion exceeds about 4 to 5; Eq. (2) and Table II write Oldroyd-B, FENE-P and other models with a conformation tensor, which allows a positive-definiteness check and log-conformation schemes; a two-dimensional Taylor–Couette simulation became chaotic near a Weissenberg number of 10. The workshop summary states that Oldroyd-B is almost never quantitative at high Weissenberg number.
B. Thomases and M. Shelley, Emergence of singular structures in Oldroyd-B fluids, Physics of Fluids 19, 103103 (2007)
Stokes–Oldroyd-B flow driven by the four-roll force without stress diffusion. Near the central stagnation point, the local strain rate times the Weissenberg number acts as an effective Weissenberg number; that strain rate tends to 2/3 as the Weissenberg number goes to zero. Critical effective values 1/3 and 1/2 correspond to Weissenberg numbers of about 0.5 (cusp in the stress) and 0.9 (unbounded stress). Grids up to 4096 points on a side. With a FENE-P penalization the stress stays bounded but cusp-like structures remain.
B. Thomases and M. Shelley, Transition to mixing and oscillations in a Stokesian viscoelastic flow, Physical Review Letters 103, 094501 (2009)
Eqs. (1)–(2): Stokes flow with Oldroyd-B stress in a periodic box, driven by a steady four-roll force; the polymer-to-solvent viscosity ratio fixed at 0.5 to match the cross-slot experiments of Arratia et al.; stress diffusion 1e-3. At a Weissenberg number of 10 the flow breaks symmetry near time 200 in the paper’s units and later develops oscillations with periods of about 317 and 59; below about 5 the flow stays slaved to the forcing.
B. Thomases, M. Shelley and J.-L. Thiffeault, A Stokesian viscoelastic flow: transition to oscillations and mixing, Physica D 240, 1602 (2011)
Eqs. (1)–(2) and (6): the same model with stress diffusion 1e-3 for Weissenberg numbers up to 10. Eqs. (4)–(5): the polymer strain energy balance. Section 3.1, Fig. 2: a small perturbation of the symmetric steady state decays at 4.8 and grows at 4.9; 256 by 256 grid, and doubling the resolution does not change the dynamics. Fig. 6(d): at 10 the velocity spectrum has periods of about 330 and 55. Section 5: between about 4.9 and 9 the asymmetric flow settles down again; above about 9 a vortex unpins and persistent oscillations mix the fluid; periodic at 12 and 15, aperiodic above about 20. Appendix B: stress diffusion is calibrated against FENE-P; its effect on the threshold is not reported.
A. Gupta and D. Vincenzi, Effect of polymer-stress diffusion in the numerical simulation of elastic turbulence, arXiv:1809.09648v2 (2019)
Version 2, 7 March 2019; the journal version was not checked. Two-dimensional Oldroyd-B flow driven by a cellular force at Deborah number 10, solved with a Cholesky decomposition and the Kurganov–Tadmor scheme. With the stress diffusion used by Thomases and co-workers (Schmidt number 1e3) the flow stays almost frozen for a long time and then breaks the cellular symmetry (Figs. 2–4); without diffusion it becomes chaotic quickly and keeps the symmetry of the forcing (Fig. 5). Section 3: in their setup the symmetry breaking is a spurious effect of the diffusion. Section 4: the same holds for FENE-P; for forcings that mix strain and vorticity the effect may be weaker.
Sandbox Physics FR26, four-roll elastic instability: numerical model (docs/frontier/fr26/four-roll-elastic-instability-model.md, with long-run records in docs/frontier/fr26/records/)
Model version fr26-four-roll-elastic-instability-v1. The model note states that the doubly periodic four-roll cell replaces the curved channel of the 2001 experiment: it has curved streamlines and a stretching stagnation point but no walls, no through-flow, no channel width and no radius, so channel thresholds and spectra must not be compared with it. It records the equations, the pseudo-spectral discretisation and the diagnostics shown in the Lab. Each run is stored with its grid size, time step, stress diffusion and perturbation seed. These are numerical results of this route, not experimental data.