3. Results

When dealing with EVP fluids in practice, a macroscopic point of view of the yield transition is usually adopted, where it is assumed that an applied shear above the yield criteria $\sigma _y$ will sufficiently shear-thin the material, allowing it to flow (Dennin Reference Dennin2008; Bonn et al. Reference Bonn, Denn, Berthier, Divoux and Manneville2017). Here, we will model directly how much of the domain remains solid (unyielded) under the imposed flow field. The solid volume fraction is calculated by $\phi (t)=\mathcal{V}_{\mathcal{F}=0}(t)/\mathcal{V} = \lvert \{k:\mathcal{F}(x_k,t)=0\}\rvert / (n_x n_y N^2)$ where ‘:’ denotes the set formed by the unyielded regions ( $\mathcal{F}=0$ ) of $k$ , and ‘ $\lvert \boldsymbol{\cdot }\rvert$ ’ denotes cardinality (Dzanic et al. Reference Dzanic, From and Sauret2024). The polymer stretching $\operatorname {tr}{\unicode{x1D63E}}$ field at various $\textit{Bi}$ is shown in figure 1. (See supplementary movies 1–4 for $\textit{Bi}=0,\ 0.5,\ 2$ and $3$ , respectively.) The volume fraction $\phi$ time series and statistics are shown in figure 2(a) and (b,c), respectively. Here, $\phi$ increases with $\textit{Bi}$ non-trivially, with fluctuations of $\phi (t)$ (yielding and unyielding) in the statistically homogenous region $\bar {t}$ varying non-monotonically with $\textit{Bi}$ (figure 2 b,c). At $\textit{Bi}=3$ , the flow slowly transitions to a nearly completely jammed state, arresting the flow (as shown figure 2 a and supplementary movie 4).

Figure 1. The EVP Kolmogorov flow in the ET regime. Representative snapshots of the polymer stretching $\operatorname {tr}{\unicode{x1D63E}}$ field $x=[0,6\times 2\pi )$ and $y=[0,4\times 2\pi )$ at various Bingham numbers, from $\textit{Bi} = 0$ to $ 3$ (columns), at different instances in time (rows): from the initial steady state $\sim\!100T$ (top) to the transition $\sim\!300T$ (middle), and the statistically homogenous regime $\overline {t}\,{\gtrsim}\, 600T$ (bottom). The grey areas represent the instantaneous unyielded $\mathcal{F}=0$ regions. The first column $\textit{Bi} = 0$ corresponds to a purely viscoelastic ET. Note, the colour bar is truncated below the true maximum at $\textrm{tr}(\boldsymbol{C}) = 1500$ for visual clarity.

Figure 2. The solid (unyielded) volume fraction ( ${\textit{a}}$ ) time series $\phi (t)$ and its fluctuation statistics in the statistically homogeneous regime $\overline {t}\,{\gtrsim}\, 600T$ , including ( ${\textit{b}}$ ) violin distribution density superimposed with box-whisker statistics relative to the median $\widetilde {\phi }$ (white line) and ( ${\textit{c}}$ ) the root-mean-squared fluctuations, where $\overline {\phi }$ is the temporal mean. The colour scheme in ( ${\textit{a}}$ ) refers to $\textit{Bi}\gt 0$ as in ( ${\textit{b}}$ ) and ( ${\textit{c}}$ ). In ( ${\textit{b}}$ ), the box plots summarise the lower and upper quartile range of $\phi (t\geqslant \overline {t})$ with violins visualising the density and shape of the distribution. Notably, the strongest and broadest distribution of fluctuations in $\phi (t\geqslant \overline {t})$ is at $\textit{Bi}=2$ . The system energy balance of the spatiotemporal mean for ( ${\textit{d}}$ ) instantaneous kinetic energy (3.1), ( ${\textit{e}}$ ) the viscous dissipation $\epsilon _\nu$ and elastic dissipation $\epsilon _p$ , and ( ${\textit{f}}$ ) their corresponding fluctuations (3.2). The directory including the data and the notebook that generated this figure can be accessed at https://www.cambridge.org/S0022112025104588/JFM-Notebooks/files/Figure_2/Fig2.ipynb.

A notable feature of viscoelastic Kolmogorov flow is the formation of coherent structures (CS) of the stress field, the aforementioned arrowhead structures (Page et al. Reference Page, Dubief and Kerswell2020), which are well-known for the purely viscoelastic case ( $\textit{Bi}=0$ ) at $Wi\gt 10$ (Boffetta et al. Reference Boffetta, Celani, Mazzino, Puliafito and Vergassola2005b ; Berti et al. Reference Berti, Bistagnino, Boffetta, Celani and Musacchio2008; Berti & Boffetta Reference Berti and Boffetta2010; Lewy & Kerswell Reference Lewy and Kerswell2025) and manifest as travelling elastic waves in the streamwise direction $t\,{\gtrsim}\,\, 300T$ (see $\textit{Bi}=0$ case in figure 1 and supplementary movie 1). For EVP cases, similar CS appear due to viscoelastic instabilities of high elasticity, where the addition of spatiotemporal interplay between its viscoelastic-solid and-fluid behaviour leads to further dynamic structural stress deformations of the CS arrowheads (see supplementary movies). These modifications to the CS, which manifest even through the introduction of minimal plasticity at $\textit{Bi}=0.25$ (figure 1), immediately imply modifications to the modes of elastic waves, altering the transition route to the chaotic self-sustaining ET state (Kerswell & Page Reference Kerswell and Page2024; Lewy & Kerswell Reference Lewy and Kerswell2025). The spectral scaling exponent $E\propto k^{-\alpha }$ at $\textit{Bi}=0$ follows the exponent $\alpha =4$ in agreement with Lellep, Linkmann & Morozov (Reference Lellep, Linkmann and Morozov2024), Lewy & Kerswell (Reference Lewy and Kerswell2025) and all cases $\textit{Bi}\leqslant 2.5$ are within the ET regime, as shown in supplemetary material, § S2, figure S3. Notably, increasing $\textit{Bi}$ progressively flattens the scaling exponent $4\gt \alpha \gt 3$ and reduces the length of the inertial subrange. The unyielded (viscoelastic-solid) regions initially ( $t\sim 100T$ , top row in figure 1) form between shear layers (i.e. low-shear-rate regions) before manifesting behind the CS front $t\,{\gtrsim}\, 300T$ . These low-shear regions effectively act as low-energy wells (Donley et al. Reference Donley, Narayanan, Wade, Park, Leheny, Harden and Rogers2023) that facilitate the formation of unyielded regions where, due to mass conservation, the stress is redistributed between shear layers (see the middle row $\sim\!300T$ in figure 1). During the initial transition $\sim\!300T$ , the length of the unyielded regions in the CS increases with $\textit{Bi}$ due to a combination of greater $\phi$ (a consequence of increased yield criteria) and increased elastic effects with increasing $\textit{Bi}$ (at a given constant $\textit{Wi}$ ), causing longer streaks. This behaviour is also apparent in the shape of unyielded regions as they become increasingly deformed and elongated. In the statistically homogeneous regime $\overline {t}\,{\gtrsim}\, 600T$ for $\textit{Bi}=0.25$ to $2.5$ , the unyielded regions manifest as local rearrangements with a broad distribution of sizes and shapes leading to further deformations of the CS (see figure 1 and figure 2 b–c) – concomitant with the view that plastic events lead to a redistribution of elastic stresses in the system (Nicolas et al. Reference Nicolas, Ferrero, Martens and Barrat2018; N’Gouamba et al. Reference N’Gouamba, Goyon and Coussot2019). For $\textit{Bi}\geqslant 1$ , unyielded regions interact across shear layers, merging or splitting each other, where increasing $\textit{Bi}$ increasingly disorganises the base flow until at the extreme $\textit{Bi}=3$ where the polymers are stretched in thin and highly localised regions. Notably, $\textit{Bi}=3$ is the only case with two distinct transitions at $\sim\!300T$ and $\sim\!800T$ , reaching a statistically homogeneous regime $\overline {t}\,{\gtrsim}\, 1000T$ (figure 2 a).

Figure 3. Features of jamming. Spatiotemporal mean of the velocity components: (a) the mean streamise flow profile $\langle \overline {U_x} (y^*)\rangle = (K n_y)^{-1}\sum _{k=0}^{K n_y-1}\lvert \langle \overline {u_x}(y^* + k\pi )\rangle _{{x}} \rvert$ along $y^* \in [0,\pi )$ , (b) $\langle \overline {u_y} (y) \rangle _{{x}}$ and (c) the spatiotemporal mean of the polymer force $\langle \boldsymbol{\overline {\boldsymbol{\nabla }\boldsymbol{\cdot }\sigma }} \rangle _{\boldsymbol{x}} = \langle \overline {\boldsymbol{F}_{p}} \rangle _{\boldsymbol{x}}$ . In (a), velocity profiles are superimposed (grey dashed lines) with base sinusoidal profile scaled by the amplitude of each $\textit{Bi}$ case, i.e. $\max _y(\langle \overline {U_x} (y)\rangle )\operatorname {cos}(Ky)$ . (d) Influence of plasticity on the energy injection rate per unit area, measured by flow resistance $\mathcal{R}$ , the ratio between the power injected in the statistically homogeneous regime to the base laminar fixed point. (e) Volume fraction $\phi (t)$ as a function of Bingham number $\textit{Bi}$ , comparing (diamonds) the temporal mean $\bar {\phi }=\langle \phi \rangle _{\overline {t}}$ numerical simulations and (circles) the theoretically approximated $\phi$ (3.3). For $\phi \lt \phi _J$ , we find that $\phi \sim \sqrt {\textit{Bi}}$ (blue line) with the linear fit $\phi =0.387 \sqrt {\textit{Bi}} -9.7\times 10^{-2}$ with a squared correlation coefficient $R^2=0.991$ . Beyond the jamming transition $\phi _J\simeq 0.54$ at $\textit{Bi}=2$ (red dash-dot line) $\phi \sim \exp {\textit{Bi}}$ (green line) with the linear fit $\phi =3.2\times 10^{-2} \exp {\textit{Bi}} - 0.24$ with ${R}^2=0.999$ . The directory including the data and the notebook that generated this figure can be accessed at https://www.cambridge.org/S0022112025104588/JFM-Notebooks/files/Figure_3/Fig3.ipynb.

Increasing $\textit{Bi}$ decreases the overall energy (figure 2 d), inferring a shift in the instantaneous kinetic energy balance $\partial E_k/\partial t\approx 0$ in the statistically homogeneous regime,

(3.1) \begin{equation} \frac {\partial E_k}{\partial t} = \langle \varepsilon _i\rangle _{\bar {t}}-\langle \varepsilon _{\nu }\rangle _{\bar {t}} - \langle \varepsilon _p\rangle _{\bar {t}} = 0, \quad t\geqslant \bar {t}, \end{equation}

where $\varepsilon _i = \langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{F}_0\rangle _{\boldsymbol{x}}$ is the energy input from the Kolmogorov forcing ( $\boldsymbol{F}_0 = (F_0 \cos (y/\ell _y), 0)$ ), which is dissipated due to contributions from both the viscous Newtonian component $\varepsilon _{\nu }=2\mu _s\langle {\unicode{x1D63F}}:{\unicode{x1D63F}}\rangle _{\boldsymbol{x}}$ and the non-Newtonian polymer component $\varepsilon _p=\langle \boldsymbol{u}\boldsymbol{\cdot }(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma })\rangle _{\boldsymbol{x}}$ . Here, ${\unicode{x1D63F}} = ({1}/{2})(\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^{\mathsf{T}})$ is the rate-of-strain tensor. The contributions of $\varepsilon _p$ (elastic and plastic behaviour) increase with $\textit{Bi}$ and eventually dominate for $\textit{Bi}\geqslant 2.5$ , absorbing more energy than the dissipation of viscous kinetic energy (figure 2 e). Fluctuations of the kinetic energy,

(3.2) \begin{equation} \frac {\partial {E_k}^{\prime }}{\partial t} = \langle {\varepsilon _i}^{\prime }\rangle _{\bar {t}}-\langle {\varepsilon _{\nu }}^{\prime }\rangle _{\bar {t}} - \langle {\varepsilon _p}^{\prime }\rangle _{\bar {t}}= 0, \quad t\geqslant \bar {t}, \end{equation}

in which the production term ${\varepsilon _i}^{\prime } =\langle {u_i}^{\prime }{u_j}^{\prime }({(\partial \overline {u_i})}/{(\partial x_j)})\rangle _{\boldsymbol{x}} \approx 0$ (figure 2 f) due to negligible inertial ( $Re=1$ ) contributions to the velocity fluctuations ${\boldsymbol{u}^{\prime }}=\boldsymbol{u}-\overline {\boldsymbol{u}}$ (Gotoh & Yamada Reference Gotoh and Yamada1984; Boffetta et al. Reference Boffetta, Celani and Mazzino2005a ), where $\overline {\boldsymbol{u}}=\langle {\boldsymbol{u}}\rangle _{\bar {t}}$ . Instead, the fluctuating non-Newtonian polymer term ${\varepsilon _p}^{\prime }=2\mu _s\langle {\unicode{x1D63F}}^{\,\prime }:{\unicode{x1D63F}}^{\,\prime }\rangle _{\boldsymbol{x}}$ is the positive source term to sustain and counteract the fluctuating viscous dissipation ${\varepsilon _{\nu }}^{\prime }=\langle \boldsymbol{u}^{\prime }\boldsymbol{\cdot }(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }^{\prime })\rangle _{\boldsymbol{x}}$ , which vary non-monotonically with $\textit{Bi}$ (figure 2 f). Interestingly, the viscous ${\varepsilon _{\nu }}^{\prime }$ and elastic dissipation ${\varepsilon _p}^{\prime }$ contributions to the energy budget clearly indicate a transition at $\textit{Bi}=2$ , which aligns directly with the strongest and broadest distribution of fluctuations in $\phi (t)$ (figure 2 b and c). Fluctuations in the dissipation then decrease for $\textit{Bi}\gt 2$ (figure 2 f) due to the polymer dissipation exceeding the viscous dissipation (figure 2 e), i.e. the system becomes unable to sustain the flow due to strong fluctuations. The non-monotonic relationship between $\textit{Bi}$ and the fluctuations $\phi (t)$ along with the impaired flow energy for $\textit{Bi}\geqslant 2$ are clear features of a typical phase transition – such a phase transition from a fluid-like to solid-like state is specifically known as jamming (Bonn et al. Reference Bonn, Denn, Berthier, Divoux and Manneville2017).

A remarkable feature of the Kolmogorov flow is that even in the chaotic ET regime, the mean velocity (figure 3 a) and conformation tensor are accurately described by sinusoidal profiles of the base flow with smaller amplitudes with respect to the laminar fixed point (Boffetta et al. Reference Boffetta, Celani and Mazzino2005a , Reference Boffetta, Celani, Mazzino, Puliafito and Vergassolab ; Berti & Boffetta Reference Berti and Boffetta2010). We observe the $\langle \overline {U_x} \rangle$ profile for all $\textit{Bi}$ retain this feature (see figure S4 in the supplementary material) with the velocity front decreasing as $\textit{Bi}$ increases (figure 3 a). Cases for $\textit{Bi}\leqslant 0.5$ show traces of CS (figure 1) resembling purely elastic flow behaviour. The decreasing front velocity amplitude $\langle \overline {U_x} \rangle$ is minor for $\textit{Bi}\leqslant 1$ but abruptly shifts for $\textit{Bi}\geqslant 2$ decreasing dramatically, approaching static rest $\langle \overline {U_x} \rangle \rightarrow 0$ at $\textit{Bi}=3$ in figure 3 a (see supplemetary movie 4). This infers a shift in the energy balance, in particular energy dissipation (figure 2 d–f), further evident from the increase in the transverse velocity component $\langle \overline {u_y} (y) \rangle _x$ and the polymer force $\boldsymbol{F}_{p} = \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }$ with increasing $\textit{Bi}$ in figures 3(b) and 3(c), respectively. The deformation and locality of plastic events redistribute stresses anisotropically (Nicolas et al. Reference Nicolas, Ferrero, Martens and Barrat2018) as observed in figures 1 and 2, where the growth of these unyielded regions is strongly correlated with secondary flow effects (figures 3 b and 3 c), leading to a departure from the Kolmogorov base sinusoidal flow profile as $\textit{Bi}$ increases (figure 3 a). Notably, the small transverse loads $F_{p,y}$ for $\textit{Bi}\gt 1$ in figure 3(c) arise due to finite deformability (figure 1), which is a well-known feature of traditional jamming in dense suspension flows (see e.g. Cates et al. Reference Cates, Wittmer, Bouchaud and Claudin1998).

We estimate the jamming transition for the present Kolmogorov flow configuration to be at the transition $\textit{Bi}=2$ , $\phi _J\simeq 0.54$ (figure 3 e), a value common in many traditional jammed systems, such as dense suspensions (Peters, Majumdar & Jaeger Reference Peters, Majumdar and Jaeger2016; Bonn et al. Reference Bonn, Denn, Berthier, Divoux and Manneville2017). Approaching the jamming transition from below $\phi \rightarrow \phi _J$ ( $\textit{Bi}\lt 2$ ), we find that $\phi \sim \sqrt {\textit{Bi}}$ in figure 3(e). Above jamming $\phi \geqslant \phi _J$ is a clear deviation from $\phi \sim \sqrt {\textit{Bi}}$ , scaling as $\phi \sim \exp {\textit{Bi}}$ ; such a dramatic change in behaviour further supports the features of phase transition observed in figure 2. To theoretically approximate $\phi$ in figure 3(e), we derive the laminar solution (see supplementary material, § S3),

(3.3) \begin{equation} {\unicode{x1D63E}}_{\textit{lam}}(y) = \begin{pmatrix} f^{-1}\left (1 + \dfrac {2 K^2Wi^2}{\mathcal{F}^2 f} \operatorname {sin}^2(Ky)\right ) & K\dfrac {Wi}{\mathcal{F} f}\operatorname {sin}(Ky) \\[12pt] K\dfrac {Wi}{\mathcal{F} f}\operatorname {sin}(Ky) & 1 \end{pmatrix}\!. \end{equation}

Figure 3(e) shows the scaling behaviour, prior to and above $\phi _J$ , is in direct agreement with $\phi$ approximated by the laminar solution (3.3). The close agreement at intermediate $\textit{Bi}\leqslant 2.5$ is surprising given the drastic flow deformations (figure 1). Interestingly, for $\textit{Bi}=3$ , while (3.3) does not predict $\langle \phi \rangle _{\bar {t}}$ in $\overline {t}\,{\gtrsim}\, 1000T$ (figure 3 e), it is consistent with $\langle \phi \rangle _{t}$ in the initial steady-state region $500T\lesssim t\lesssim 800T$ (figure 2 a). Analogue scaling behaviour $\phi \sim \sqrt {\textit{Bi}}$ for $\phi \leqslant \phi _J$ is commonly observed in traditional jamming transitions of dense suspensions (Peters et al. Reference Peters, Majumdar and Jaeger2016; Bonn et al. Reference Bonn, Denn, Berthier, Divoux and Manneville2017). Recent findings by Abdelgawad et al. (Reference Abdelgawad, Cannon and Rosti2023) showed that EVP fluids in inertial turbulence ( $Re\gg 10^3$ ), with subdominant elastic effects $Wi\lesssim 1$ , plasticity $\textit{Bi}\gt 1$ increases intermittency. We find their results for $\textit{Bi}\gt 1$ agree with our scaling $\phi \sim \sqrt {\textit{Bi}}$ approaching the jamming transition (as shown in supplementary material, figure S5). Whether an analogue exponential scaling beyond the jamming transition $\phi \gt \phi _J$ is relevant in such inertia-dominated turbulent regimes remains to be observed.

The potential implications in practical applications of this elastic-plastic-induced jamming will be inferred by measuring the flow resistance to the power injected (Berti et al. Reference Berti, Bistagnino, Boffetta, Celani and Musacchio2008). The energy injection rate per unit area (power injected), $ \mathcal{P}_{\textit{inj}} = \langle \overline {\boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{F}} \rangle _{\boldsymbol{x}}$ , with $\boldsymbol{F}=\boldsymbol{F}_0+\boldsymbol{F}_{\!p}$ . In the laminar steady state, denoted by $\widehat {\boldsymbol{\cdot }}$ , $\partial _x \widehat {P}$ is constant and $\langle u_x \rangle _{\boldsymbol{x}} \approx \widehat {u}$ , where $\widehat {\mathcal{P}_{\textit{inj}}} = ({1}/{2}) V F_0$ (see supplemetary material, § S3.1). Due to the body force and the energy dissipated by the elastic-plastic-viscous effects, the flow resistance, $\mathcal{R}=\mathcal{P}_{\textit{inj}}/\widehat {\mathcal{P}_{\textit{inj}}} = \langle \overline {\boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{F}} \rangle _{\boldsymbol{x}}/({1}/{2} V F_0)$ , gradually increases for $\textit{Bi}\leqslant 1$ before shifting as a consequence of the jamming phase change (figure 3 d). Moreover, our jamming transition scaling behaviour $\phi \sim \sqrt {\textit{Bi}}$ has a potentially direct connection with other known scaling, namely, the Fanning friction factor in channel flows (Rosti et al. Reference Rosti, Izbassarov, Tammisola, Hormozi and Brandt2018) and the pressure drop in porous media flows (De Vita et al. Reference De Vita, Rosti, Izbassarov, Duffo, Tammisola, Hormozi and Brandt2018).

The macroscopic point of view of jamming in EVP materials is typically considered as the yield criteria, an immediate shift to a flow state once a certain stress threshold is reached, e.g. if the flow applied induces sufficiently high shear (Dennin Reference Dennin2008; Bonn et al. Reference Bonn, Denn, Berthier, Divoux and Manneville2017). We demonstrate the statistically homogenous dynamics of EVP Kolmogorov flows in the ET regime (dominant elastic instabilities) across a range of $\textit{Bi}$ ( $\phi$ ) align with the laminar scaling prediction below and above the transition $\phi _J$ , i.e. $\phi \sim \sqrt {\textit{Bi}}$ and $\phi \sim \exp {\textit{Bi}}$ , respectively (figure 3 e). The connection to jamming in this work has profound implications on EVP fluid flows; it describes plastic events, such as plastic ‘plugs’ in channel flows (Rosti et al. Reference Rosti, Izbassarov, Tammisola, Hormozi and Brandt2018; Izbassarov et al. Reference Izbassarov, Rosti, Brandt and Tammisola2021; Villalba et al. Reference Villalba, Daneshi, Chaparian and Martinez2023), as a phase-change phenomenon, whose consequences in dynamic behaviour are analogous to traditional jamming but different in the sense that it is induced by dynamic structural stress deformations arising from viscoelastic flow instabilities. This dependency on the locality of these plastic-induced stress deformations is all important as these depend on the local geometry and flow field, implying that the nature of this jamming transition $\phi \rightarrow \phi _J$ varies from one configuration to another. For example, for the square-root scaling in the form $\phi \sim c_1 \sqrt {\textit{Bi}} + c_2$ (see e.g. figure 3), the constants $c_1$ and $c_2$ depend on the flow configuration, $Wi$ , and $Re$ . Similarly, the exponential scaling $\phi \sim \exp \textit{Bi}$ may also be influenced by such factors, and its precise form could vary across flow configurations. In support of this, our previous work (Dzanic et al. Reference Dzanic, From and Sauret2024) found that two extensional flow benchmark problems, with the same dimensionless variables and similar flow-type distribution, result in very different dynamic responses to plasticity with different $\phi$ across a range of $\textit{Bi}$ . The consequence of this dependency, for example, is that the bulk shear rheology characterisation of the material (Cheddadi, Saramito & Graner Reference Cheddadi, Saramito and Graner2012; Varchanis et al. Reference Varchanis, Haward, Hopkins, Syrakos, Shen, Dimakopoulos and Tsamopoulos2020) – commonly performed in Couette-type geometries – will be subject to jamming dynamics different from those experienced in the actual flow configuration of the application, making their performance in practice unpredictable. Such issues in translating rheological characterisation to quantify flow performance in actual flow configurations have recently been reported as major challenges in predicting the flow performance of functional EVP materials (Bertsch et al. Reference Bertsch, Diba, Mooney and Leeuwenburgh2022). Predicting and controlling jamming is crucial in, e.g. 3-D extrusion biomaterial printing, where the jamming of soft materials during extrusion has a negative impact on both the print quality and cell viability (Xin et al. Reference Xin, Deo, Dai, Pandian, Chimene, Moebius, Jain, Han, Gaharwar and Alge2021).