1. Introduction
High molecular weight polymers in solution can break in strong extensional flows. Industrial and research-scale fluidic systems are prone to this effect whenever they feature high flow rates of macromolecule solutions through pipes, constrictions or porous media. For example, polymers used for drag reduction purposes or in enhanced oil recovery fluids are mechanically degraded over time (Seright Reference Seright1983; Al-Shakry et al. Reference Al-Shakry, Skauge, Shiran and Skauge2018; Soares Reference Soares2020). Large biopharmaceuticals can break or lose their activity when processed through narrow channels (Lengsfeld & Anchordoquy Reference Lengsfeld and Anchordoquy2002; Rathore & Rajan Reference Rathore and Rajan2008; Cook, Wang & Derby Reference Cook, Wang and Derby2010; Hawe et al. Reference Hawe, Wiggenhorn, van de Weert, Garbe, Mahler and Jiskoot2012). Long polymers incorporated in inkjet inks can be degraded in the jetting flow or recirculation circuit, changing the functional properties or the printability of the fluid (A-Alamry et al. Reference A-Alamry, Nixon, Hindley, Odel and Yeates2010; McIlroy, Harlen & Morrison Reference McIlroy, Harlen and Morrison2013). On the other hand, extensional forces in flows can be used on purpose to fragment DNA in next generation sequencing technologies (Shui et al. Reference Shui, Bomer, Jin, Carlen and van den Berg2011), or more broadly to activate force-sensitive compounds in the context of polymer mechanochemistry (May & Moore Reference May and Moore2013; Willis-Fox et al. Reference Willis-Fox, Rognin, Aljohani and Daly2018, Reference Willis-Fox, Rognin, Baumann, Aljohani, Göstl and Daly2020).
Although extensive research has been done on flow-induced polymer scission, advanced modelling and simulations mainly focus on the molecular scale. Bond rupture can be simulated via quantum chemistry (Stauch & Dreuw Reference Stauch and Dreuw2016), bond angles motion and short chains bending via all-atoms molecular dynamics (Ribas-Arino & Marx Reference Ribas-Arino and Marx2012), and unravelling of long chains of polymers via coarse-grained molecules in implicit solvent, such as bead-spring or bead-rod models (Knudsen, Martínez & Torre Reference Knudsen, Martínez and de la Torre1998; Maroja et al. Reference Maroja, Oliveira, Cieśla and Longa2001; Hsieh, Park & Larson Reference Hsieh, Park and Larson2005; Sim, Khomami & Sureshkumar Reference Sim, Khomami and Sureshkumar2007). Some approaches to model polymer degradation at continuum scale have been focused on predicting the onset of chain scission rather than the complete description of the reaction kinetics. For turbulent flows at high Reynolds numbers, experiments have revealed that the onset of degradation is a function of the Reynolds number and physico-chemical parameters of the polymer–solvent system only, so that the exact geometry of the ducts as well as the detailed patterns of the flow are irrelevant to the problem (Nguyen & Kausch Reference Nguyen and Kausch1991; Vanapalli, Ceccio & Solomon Reference Vanapalli, Ceccio and Solomon2006). In the field of enhanced oil recovery fluids, models have been developed to predict the time evolution of the polymer average molecular weight (Sorbie & Roberts Reference Sorbie and Roberts1984; Brakstad & Rosenkilde Reference Brakstad and Rosenkilde2016; Lohne et al. Reference Lohne, Nødland, Stavland and Hiorth2017), but have remained largely empirical, because the degradation mechanism is modelled by macroscopic averaged quantities, and not at pore scale. In a recent work Garrepally et al. (Reference Garrepally, Jouenne, Olmsted and Lequeux2020) used multiple passes though a microfluidic constriction to study polymer degradation via pressure losses. They found scaling behaviours and could predict polymer degradation for a range of flow rates. However, it is not clear how their findings would translate directly to other constriction geometries or other types of flow fields.
On the other hand, efforts to model polymer degradation in terms of local reaction kinetics and velocity fields have been limited so far. López Cascales & García de la Torre (Reference López Cascales and García de la Torre1992) used coarse-grained molecular models to study the kinetics of rupture of large ensembles of chains in a sudden elongational flow. They found two steps in the degradation process: a first period of time without damage, corresponding to the unravelling of the molecules, and a second step with damage well described by first-order reaction kinetics. Although the scission rate was analysed in terms of the strain rate of the flow and molecular lengths, the model was not generalized to arbitrary flows in a form suitable for computational fluid dynamics (CFD). More recently, Pereira, Mompean & Soares (Reference Pereira, Mompean and Soares2018) presented a series of simulations of polymer degradation in turbulent flows. In their approach, the contour length of the polymer is not a constant parameter but a scalar field convected by the flow. Degradation is simulated by an arbitrary geometric criterion: the local contour length (or molecular weight) is reduced by a small amount as soon as the average conformation reaches an extensibility threshold. The viscoelastic stress is computed based on the local polymer length, and the turbulent flow is solved by a Direct Navier–Stokes (DNS) approach. Although this work is an interesting proof of concept that a local modelling of chain scission is achievable even in the case of turbulent flows, the underpinning scission model was not derived from mechanochemical principles, and the generalization and predictive ability of this approach is still to be confirmed.
However, with the move to integrate mechanically activated functional macromolecules into materials and processes, it is critical to develop a quantitative model that can predict the rate of chain scission while being geometry agnostic and integrated within a CFD analysis. For example, such elaborate models exist for the specific case of worm-like micellar solutions, for which closed-form chemical kinetics for scission and recombination of the micelles have been developed (Carl, Makhloufi & Kröger Reference Carl, Makhloufi and Kröger1997; Vasquez, McKinley & Cook Reference Vasquez, McKinley and Cook2007; Germann, Cook & Beris Reference Germann, Cook and Beris2013; Dutta & Graham Reference Dutta and Graham2018).
In this work we validate through experimental results a continuum mechanochemical model of polymer chain scission containing a minimal number of parameters. We introduce a set of partial differential equations inspired by rheological constitutive equations and molecular simulations previously reported (Rognin et al. Reference Rognin, Willis-Fox, Aljohani and Daly2018). The equations are implemented and solved with the open source CFD library OpenFOAM. The model is validated using a series of scission experiments of high molecular weight polymer solutions flowing through narrow constrictions. Linear poly(ethylene oxide) (PEO) in water solutions are used for this validation. Provided that near-equilibrium properties (zero-shear viscosity, Zimm relaxation time) can be measured beforehand with standard lab tools, the present mechanochemical model needs only two parameter fits for a given polymer–solvent system: the first one is a critical strain rate calibrating the amount of chain scission, and the second one is a maximum extensional viscosity impacting the pressure loss and viscoelastic flow pattern.
The paper is organised as follows: in the next section we present the modelling framework; then the experimental method is reported, followed by an introduction to the simulation work. Experimental results and simulations are then compared and discussed. Final remarks and the potential impact across disciplines and applications are presented in the conclusion.
2. Model
An abundant literature partly cited above has shown that polymer backbone mechanical scission depends on three main components. The first one is the strength of the bond itself, usually expressed as a critical force (in nanoNewtons). The second component is the amount of strain rate in the fluid, as this sets the tension in polymer chains, imparted by the solvent through viscous friction. The third component is the conformation state of the polymer molecules, as significant tension can only build up in long straight segments. Our model, first presented in Rognin et al. (Reference Rognin, Willis-Fox, Aljohani and Daly2018), includes those three aspects in a way that is summarised below.
2.1. Mechanochemical model
We consider an initially monodisperse population of polymer chains in dilute solution. The conformation state of the linear molecules is described by the second tensorial moment 
$\boldsymbol{\mathsf{C}}$ of the end-to-end vectors 
$\boldsymbol {R}$, 
$\boldsymbol{\mathsf{C}}=\left \langle \boldsymbol {RR}\right \rangle$, with the evolution equation
\begin{equation} \frac{\mathrm{D}\boldsymbol{\mathsf{C}}}{\mathrm{D}t}=\left(\boldsymbol{\nabla}\boldsymbol{u}\right)^{\textrm{T}}\boldsymbol{\cdot}\boldsymbol{\mathsf{C}}+\boldsymbol{\mathsf{C}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}-\left(\frac{7}{2}-\frac{3\mathrm{tr}(\boldsymbol{\mathsf{C}})}{2L^{2}}\right)\frac{\left\lfloor \boldsymbol{\nabla}\boldsymbol{u}:\boldsymbol{\mathsf{C}}\right\rfloor }{L^{2}}\boldsymbol{\mathsf{C}}-\frac{\boldsymbol{\mathsf{C}}-\dfrac{R_{EE}^{2}}{3}\boldsymbol{\mathsf{I}}}{\tau_{Z}}, \end{equation}
where 
$\mathrm {D/D}t$ is the Lagrangian time derivative, 
$\boldsymbol {u}$ is the fluid velocity field, 
$L$ is the polymer contour length, 
$R_{EE}$ is the root mean square end-to-end distance of the unperturbed chains, 
$\boldsymbol{\mathsf{I}}$ is the unit tensor and 
$\tau _{Z}$ is the relaxation time for the coil-stretch transition (Zimm time, usually). The brackets 
$\lfloor \boldsymbol {\cdot }\rfloor$ indicate that only positive values of the product 
$\boldsymbol {\nabla }\boldsymbol {u}:\boldsymbol{\mathsf{C}}$ are kept in the evolution equation. This introduces a hysteresis behaviour describing the far-from-equilibrium stretching dynamics of flexible polymers, known as kink dynamics (Larson Reference Larson1990; Hinch Reference Hinch1994). When 
$\boldsymbol {\nabla }\boldsymbol {u}:\boldsymbol{\mathsf{C}}>0$, chains are stretching but in a non-affine manner with respect to the fluid because of their various folding states. In addition, this formulation naturally yields a finite extensibility of the polymer chains, with 
$\mathrm {tr}(\boldsymbol{\mathsf{C}})$ always smaller than 
$L^{2}$. By contrast, when 
$\boldsymbol {\nabla }\boldsymbol {u}:\boldsymbol{\mathsf{C}}<0$, chains are contracting along their principal axis and no significant mechanism can prevent an affine recoiling. Hence, the term involving 
$\boldsymbol {\nabla }\boldsymbol {u}:\boldsymbol{\mathsf{C}}<0$ should vanish.
A normalised and more suitable form of this equation for numerical purposes is obtained by dividing 
$\boldsymbol{\mathsf{C}}$ by 
$L^{2}$,
\begin{equation} \frac{\mathrm{D}\boldsymbol{\mathsf{A}}}{\mathrm{D}t}=\left(\boldsymbol{\nabla}\boldsymbol{u}\right)^{\textrm{T}}\boldsymbol{\cdot}\boldsymbol{\mathsf{A}}+\boldsymbol{\mathsf{A}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}-\frac{7-3\lambda}{2}\left\lfloor \boldsymbol{\nabla}\boldsymbol{u}:\boldsymbol{\mathsf{A}}\right\rfloor \boldsymbol{\mathsf{A}}-\frac{\boldsymbol{\mathsf{A}}-\dfrac{1}{3\xi^{2}}\boldsymbol{\mathsf{I}}}{\tau_{Z}},\end{equation}
where 
$\boldsymbol{\mathsf{A}}=\boldsymbol{\mathsf{C}}/L^{2}$, 
$\lambda =\mathrm {tr}(\boldsymbol{\mathsf{A}})$ is the normalised mean square polymer extension and 
$\xi =L/R_{EE}$ is the polymer extensibility. By monitoring the conformation of intact polymer chains only, 
$\xi$ and 
$\tau _{Z}$ remain constant throughout the flow. How viscoelastic stress is affected by chain scission will be described in § 2.3.
Let 
$c$ be the mass concentration field of intact polymer chains (i.e. chains of initial molecular weight). The flow-induced scission is modelled by a first-order reaction,
\begin{equation} \frac{\mathrm{D}c}{\mathrm{D}t}={-}kc,\end{equation}
where 
$k$ is the reaction rate, and where polymer diffusivity and shear-induced migration are neglected with respect to advection (Graham Reference Graham2011). Neglecting diffusivity is a safe assumption where chains break since the shear rate is already strong enough to overcome fast internal recoiling mechanisms. We define the local degradation field, 
$\varphi$, ranging from 0 to 1, 1 being the case of complete degradation, by
\begin{equation} \varphi=1-\frac{c}{c_{0}},\end{equation}
where 
$c_{0}$ is the initial mass concentration intact polymers.
In our previous study we obtained a closed form for 
$k$ inspired by molecular simulations (Rognin et al. Reference Rognin, Willis-Fox, Aljohani and Daly2018). Although the original expression depends on Lagrangian time derivatives, here, we suggest the following simplified form:
\begin{equation} k=\begin{cases} 0 & \text{if }\left\lfloor \boldsymbol{\nabla}\boldsymbol{u}:\boldsymbol{\mathsf{A}}\right\rfloor /\lambda<\dot{\varepsilon}_{c},\\ \alpha\dfrac{\left\lfloor \boldsymbol{\nabla}\boldsymbol{u}:\boldsymbol{\mathsf{A}}\right\rfloor ^{2}}{\dot{\varepsilon}_{c}\lambda} & \text{otherwise}. \end{cases}\end{equation}
Here 
$\alpha$ is a coefficient of the order of unity, 
$\dot {\varepsilon }_{c}$ is a critical strain rate and 
$\lfloor \boldsymbol {\nabla }\boldsymbol {u}:\boldsymbol{\mathsf{A}}\rfloor /\lambda$ is a measure of the strain rate in the direction of polymer elongation. The first line of (2.5) describes the cut-off case where the strain rate is lower than the critical value. From a theoretical point of view, if the strain rate is exactly 
$\dot {\varepsilon }_{c}$ then the scission rate should tend towards the rate at which chains approach their fully unravelled state, because scission events then occur only in stretched conformation. The second line of (2.5) describes this threshold rate as tending towards 
$\alpha \dot {\varepsilon }_{c}$ when 
$\lambda \rightarrow 1$. Free-draining bead-rod models give 
$\alpha \approx 0.5$, which is the value adopted in this work, but it can be anticipated that, for real chains, 
$\alpha$ would be lower because of the hydrodynamically hindered unravelling dynamics (Hsieh & Larson Reference Hsieh and Larson2004). If the strain rate is higher than 
$\dot {\varepsilon }_{c}$ then scission can occur even in non-fully stretched chains, and at a rate which depends quadratically on the strain rate.
Note that the scission rate falls back to zero instantaneously when the flow is switched off. For this to be valid, the scission rate should be smaller than the relaxation rate of mechanical tension. If the Zimm time, 
$\tau _{Z}$, is considered as the longest relaxation time, then tension relaxation driven by segmental diffusion at short times should happen at a fraction of 
$\tau _{Z}$. More specifically, because the number of segments scales as 
$\xi ^{2}$, segmental relaxation time should scale as 
$\xi ^{-3}\tau _{Z}$ for a Zimm chain in theta solvent, or be even smaller for free-draining chains. We will see in the following experimental section that, for the range of molecular weights studied here, 
$\dot {\varepsilon }_{c}\tau _{Z}\sim 10$. Therefore, our assumption holds provided that 
$\xi ^{3}\gg 10$, which is straightforward for long flexible molecules (here, 
$\xi \geq 19$).
In addition, this model is distinct from other flow-induced scission models reported in the literature. First, it contrasts with the original simulation work of López Cascales & García de la Torre (Reference López Cascales and García de la Torre1992) on short polymers, and current models of flow-induced scission of worm-like micelles such as the Vasquez–Cook–McKinley (VCM) model (Vasquez et al. Reference Vasquez, McKinley and Cook2007; Dutta & Graham Reference Dutta and Graham2018), where the strain rate dependence is linear. We can explain this difference by noting that the strain rate plays a double role in long coiled molecule dynamics, first by setting the friction force along the body as in the case of micelles, but also as the rate of unravelling and growth of long straight segments – negligible for micelles and short polymers. The present model also contrasts with the thermally activated bond scission (TABS) theory, where the first-order scission rate can be cast in the form
\begin{equation} k_{{TABS}}=k_{0}\exp\left(\frac{\dot{\varepsilon}}{k_{1}}\right), \end{equation}
where 
$k_{0}$ is the degradation rate without flow, 
$\dot {\varepsilon }$ is the strain rate and 
$k_{1}$ is a parameter depending on polymer properties (Odell, Keller & Muller Reference Odell, Keller and Muller1992). In the TABS model the mechanical tension, which is assumed to be proportional to the strain rate, acts to reduce the activation barrier of thermally induced bond scission. Nevertheless, the averaging of the TABS kinetics, motivated by physical arguments at the bond scale, to a population of unravelling chains experiencing different tensions due to their own folding states, is not trivial until all chains are completely unravelled. Stretching dynamics can be partly introduced by letting 
$k_{1}$ be proportional to 
$\lambda ^{-1}$, as internal tension would then be assumed proportional to 
$\langle \boldsymbol {R}^{2}\rangle \dot {\varepsilon }$. Yet, for a single-pass constriction flow, 
$k_{0}$ would be very small compared with the inverse characteristic residence time. Therefore, to observe any scission, the multiplication factor 
$\exp (\dot {\varepsilon }/k_{1})$ would need to be large in some parts of the flow, which because of its exponential nature, would impart a dramatic change in scission rate as the strain rate smoothly increases. The resulting overall kinetics is that of a thresholding, where one part of the flow experiences negligible scission while all chains break up instantaneously in the remaining part of higher strain rates.
We will assess the differences between these models in the discussion section. To summarise, the flow-induced scission model proposed in this study has the following expected properties.
(i) Straining time before rupture depends on initial configuration of individual chains, where chains that are already aligned with the elongation axis break first while chains having the most complex kinks and unravelling configuration take more time. This property is given by the first-order kinetics of (2.3).
(ii) Scission occurs only above a certain strain rate threshold corresponding to a critical tension in an fully stretched chain ((2.5), condition 1).
(iii) The scission is faster when the strain rate is larger (influence of
$\boldsymbol {\nabla }\boldsymbol {u}$) or when the chains are on average unravelled (influence of $\boldsymbol{\mathsf{A}}$) as given by (2.5), condition 2.(iv) We also assume that there is no recombination of broken chains.
2.2. Viscoelastic model
Since the present experiments are carried out at dilute but finite concentrations, viscoelastic effects are expected to play a role in elongational flow patterns. The total fluid stress, 
$\boldsymbol {\sigma }$, is decomposed into
\begin{equation} \boldsymbol{\sigma}={-}p\,\boldsymbol{\mathsf{I}}+\boldsymbol{\tau}_{s}+\boldsymbol{\tau}_{p}, \end{equation}
where 
$p$ is the pressure, 
$\boldsymbol {\tau }_{s}$ is the viscous stress due to the solvent and 
$\boldsymbol {\tau }_{p}$ is the viscoelastic stress due to the polymer. The solvent is Newtonian and assumed incompressible, so that
\begin{equation} \boldsymbol{\tau}_{s}=\eta_{s}\left(\boldsymbol{\nabla}\boldsymbol{u}+\left(\boldsymbol{\nabla}\boldsymbol{u}\right)^{\textrm{T}}\right), \end{equation}
where 
$\eta _{s}$ is the solvent viscosity. As for the polymer stress, the model has to be consistent with the evolution equation (2.2) of the conformation tensor, where the finite extensibility and non-affine deformation of the chains is expressed through the friction term 
$\boldsymbol {\nabla }\boldsymbol {u}:\boldsymbol{\mathsf{A}}$. Although this approach is common for suspensions of rigid rods in strong flows, few options have been studied for long flexible chains (Larson Reference Larson1990; Hinch Reference Hinch1994; Rallison Reference Rallison1997; Verhoef, van den Brule & Hulsen Reference Verhoef, van den Brule and Hulsen1999, see also § 5.5.3 in Larson Reference Larson1988). Here we select the form
\begin{equation} \boldsymbol{\tau}_{p}=\left(\eta_{p}(\varphi)\frac{3\xi^{2}}{\tau_{Z}}+\eta_{E}(\varphi)\left\lfloor \boldsymbol{\nabla}\boldsymbol{u}:\boldsymbol{\mathsf{A}}\right\rfloor \right)\boldsymbol{\mathsf{A}},\end{equation}
where 
$\eta _{p}(\varphi )$ is the additional zero-shear viscosity due to the polymer and 
$\eta _{E}(\varphi )$ is a maximum extensional viscosity. Both 
$\eta _{p}$ and 
$\eta _{E}$ are functions of the local degradation to account for the fact that chain scission induces a decrease in polymeric viscoelasticity. The selected functional forms will be described in the next subsection. The first term of the stress accounts for the viscoelasticity at small shear rate where polymer chains are close to equilibrium, while the second term accounts for the dissipative dynamics of chains being unravelled. This second term is in fact a viscous term of fourth rank tensorial viscosity 
$2\eta _{E}{\mathsf{A}}_{ij}{\mathsf{A}}_{kl}$. With this approach, and contrary to most finite extensible nonlinear elastic (FENE) models, the geometric extensibility (
$\xi$) is decoupled from the extensional viscosity. The behaviour of this viscoelastic stress model (2.9) and the FENE-P model (FENE with Peterlin closure) commonly used for dilute polymer solutions (Larson & Desai Reference Larson and Desai2015) are compared in more detail in § 1 of the supplementary material (SM) available at https://doi.org/10.1017/jfm.2021.646. We will discuss further the influence of viscoelastic stresses on polymer degradation in § 5.6.
2.3. Mixture properties
The stress model takes into account degraded polymeric viscosities due to chain scission. Regarding the zero-shear viscosity, we assume that the Mark–Houwink–Sakurada law applies, which for the initial (i.e. non-degraded) solution gives
\begin{equation} \eta_{p0}=c_{0}[\eta]_{0}\eta_{s}=c_{0}KM_{0}^{a}\eta_{s}, \end{equation}
where 
$\eta _{p0}$ is the added polymeric viscosity of the initial solution, 
$c_{0}$ is the mass concentration of polymers, 
$[\eta ]_{0}$ is the initial intrinsic viscosity, 
$M_{0}$ is the initial polymer molecular weight, and 
$K$ and 
$a$ are parameters depending on the polymer–solvent system and temperature. Assuming perfect halving of the chains, the Mark–Houwink–Sakurada law would give an added viscosity divided by 
$2^{a}$ for a completely degraded solution. Therefore, for a mixture of both initial and degraded polymer, we have
\begin{equation} \eta_{p}(\varphi)=\varphi\frac{\eta_{p0}}{2^{a}}+\left(1-\varphi\right)\eta_{p0}. \end{equation}In addition, according to a simple free-draining bead-rod model, the maximum extensional viscosity should be divided by four if the chains are halved (molecular contribution divided by eight, molar concentration doubled); therefore,
\begin{equation} \eta_{E}(\varphi)=\varphi\frac{\eta_{E0}}{4}+\left(1-\varphi\right)\eta_{E0}, \end{equation}
where 
$\eta _{E0}$ is the maximum extensional viscosity of the initial solution.
A complete mixture model would also account for changes in relaxation time and extensibility. Indeed, shorter chains should relax faster than the original ones, and they would require higher strain rates to undergo coil-stretch transition. However, we assume that these changes would not be significant here, as we are modelling single stretching events (single-pass contraction flow).
The next section describes the experimental system used to validate this choice of modelling approach and equations.
3. Experiments
3.1. Fluidic system
In this study a single pass of polymer solutions through narrow constrictions is considered. The fluidic system is presented in figure 1. A glass syringe (SGE Gas Tight) mounted on a syringe pump (KDS) is used to push the polymer solution at a specified flow rate inside stainless steel tubings (Valco fittings), and through one of the following two kinds of constriction.
(i) The first kind of constriction is a fused glass capillary nozzle (hereinafter referred to as fuse capillary, made from a straight capillary (Microcaps 1-000-0090, 480
$\mathrm {\mu }\textrm {m}$ internal diameter) and fused in a capillary puller (Narishige PC-10) for 75 s. The capillary is mounted on a stainless steel filter (Valco fittings and 2 microns screen). The flow through this constriction can be monitored by a camera mounted on a microscope objective.(ii) The second kind of constriction is a sharp bore through a thin plate (Edmund Optics stainless steel pinholes). In this study two different nominal diameters 25
$\mathrm {\mu }\textrm {m}$ and 50 $\mathrm {\mu }\textrm {m}$ are used (referred to as pinhole 25 and pinhole 50). The pinhole is stacked together with an inlet filter (Valco 2 micron stainless steel frit), a support washer and custom made PTFE rings, inside a filter holder (Millipore 13 mm diameter). The effect of in-line filters on degradation measurements is negligible for the vast majority of flow rates, as analysed in SM.
Figure 1. Experimental set-up and nozzle assembly.
A pressure sensor (Honeywell MLH series, wetted parts: stainless steel 304L and Haynes 214 alloy) is used to monitor the pressure upstream of the constriction. For each specified flow rate, the steady-state pressure is recorded for both pure water and polymer solutions, and used to calibrate the simulated geometries and fluid properties, as we will see in § 4.
3.2. Sample preparation and production
Poly(ethylene oxide) in solid beads form, of three molecular weights (1 MDa, 600 kDa and 300 kDa nominal molecular weight, Sigma-Aldrich) are dissolved in water (analytical reagent grade, Fisher, conductivity 1.5 
$\mathrm {\mu }\textrm {S}\,\textrm {cm}^{-1}$). Water is poured on top of the beads and the polymer is left to dissolve during two to four weeks giving a manual swirl several times. Then the solutions are filtered through a 1.2 
$\mathrm {\mu }\textrm {m}$ cellulose ester membrane (Millipore) and their zero-shear viscosity measured with a cone-plate rheometer (Anton Paar MCR302, shear rates ranging from 50 to 500 s
$^{-1}$). For each molecular weight, the concentration of polymer is chosen to be approximately a third of the overlap concentration. In particular, this is achieved when there is a 33 % increase in viscosity from the original pure solvent due to the presence of polymer. The viscosity-average molecular weight, 
$M$, of each polymer is measured by first measuring the intrinsic viscosity of (unfiltered) solutions. Using the Mark–Houwink–Sakurada law (2.10), 
$M$ is then back calculated. The PEO standards of known molecular weight (1 MDa, 500 kDa and 200 kDa nominal value, polydispersity indices of 1.10, 1.05 and 1.08, respectively, Agilent) are used to measure the Mark–Houwink–Sakurada parameters at 
$20\,^{\circ }\textrm {C}$: 
$K=0.020\pm 0.007\ \textrm {ml}\,\textrm {g}^{-1}$ (with 
$M$ in 
$\textrm {g}\,\textrm {mol}^{-1}$), 
$a=0.695\pm 0.008$. The value of 
$a$ shows that water is a good solvent at 
$20\,^{\circ }\textrm {C}$. Sample properties are summarised in table 2.
Degradation experiments are carried out as follows. A nozzle is fitted on the fluidic system. For each molecular weight, the system is rinsed first manually with 5 ml of the solution (this includes a purge of the channel leading to the pressure sensor), then using the syringe pump at low flow rate (10 
$\textrm {ml}\,\textrm {h}^{-1}$) with 2 ml. Then several flow rates are set in increasing order. For each flow rate, the first 0.9 ml of the outlet solution is discarded in order to let the pressure stabilise and to purge the fluid from the previous flow rate, then 1.5 ml of solution is collected for viscosity measurement, and a new flow rate is set. Three series are carried out for each polymer and each nozzle. Experiments are carried out in an air-conditioned room at 
$20\,^{\circ }\textrm {C}$, without additional control of the fluidic system temperature.
3.3. Quantification of polymer degradation
In this study polymer degradation is assessed by measuring the decrease in zero-shear viscosity of solutions after they flow through the constrictions. Neglecting polymer adsorption in the system and assuming chain halving, the zero-shear added viscosity of the collected solution, 
$\eta _{p}$, can be linked to the initial added viscosity according to (2.11), so that
\begin{equation} \eta_{p}=\varPhi\frac{\eta_{p0}}{2^{a}}+\left(1-\varPhi\right)\eta_{p0},\end{equation}
where 
$\varPhi$ is the overall proportion of broken chains in the collected sample. Here 
$\varPhi$ is found to be proportional to the decrease in zero-shear added viscosity 
$\eta _{p0}-\eta _{p}$,
\begin{equation} \varPhi=\frac{2^{a}}{2^{a}-1}\,\frac{\eta_{p0}-\eta_{p}}{\eta_{p0}}.\end{equation}4. Simulations
4.1. Geometries characterisation and mesh generation
In a move to experimentally validate a model that can be generalized to arbitrary flows, we have characterised the nozzles described in this report for integration into the model. In each case, an axisymmetric geometry was assumed.
4.1.1. Fused glass capillary
To avoid refraction of visible light by the outer curvature of the capillary, X-ray imaging (ZEISS Xradia Versa 520) is used to measure the inner profile of the constriction (see figure 2(a i) and figure S4 in SM for a full size view). A surface line is interpolated with Bezier curves and used to build the CFD mesh (figures S6). The minimum radius of the constriction is left as an adjustable parameter of the mesh to fit experimental pressure losses.
Figure 2. Geometries and meshes convergence properties. (a) Fused capillary. (b) Pinhole 25. (c) Pinhole 50.
4.1.2. Pinholes
Scanning electron microscope (SEM) images of the pinholes inlet reveal a salient rim of a few microns thick; see figures 2(b i) and 2(c i). Interferometry (Veeco NT3300) is used to measure a three-dimensional shape of the nozzle inlet (figures S10a and S15a). An axisymmetric-averaged profile is extracted from the surface map (figures S10b and S15b) and used for the CFD mesh. The outlet is modelled as a sharp 90
$^{\circ }$ corner, as no significant polymer scission is expected at the outlet. The radius of the hole is an adjustable parameter of the mesh to fit experimental pressure losses.
4.1.3. Mesh resolution
To produce the final version of the meshes, the following method is applied for each geometry.
(a) A first mesh refinement study is carried out by simulating the flow of water at high flow rate (where boundary layers are expected to be the thinnest). The convergence of the steady-state pressure drop is analysed upon mesh refinement. Convergence is assumed if the pressure drop difference between two meshes falls under typical experimental uncertainty (0.5 bar at high flow rate); see figures 2(a ii), 2(b ii) and 2(c ii). Reference mesh sizes are reported in table 1.
(b) The pressure drop through the constriction is simulated for the flow of water at experimental flow rates, and the constriction radius is adjusted to match experimental measurements; see figures S8, S13 and S18 in SM. Results of fitted diameters are reported in table 1.
Table 1. Mesh sizes and diameters fit.
The convergence of the meshes regarding simulated degradation was not systematically assessed, because of the computational cost of running the full model on finer meshes. From our previous computational study on Newtonian flows (Rognin et al. Reference Rognin, Willis-Fox, Aljohani and Daly2018), we expect that mesh resolution only affects cases with low overall degradation, because polymer scission then occurs in a small number of cells. Nonetheless, a mesh resolution convergence test was done for the degradation of PEO 1000k in the capillary geometry at 
$80\ \textrm {ml}\,\textrm {h}^{-1}$ (43 % overall degradation). This test suggests that the span due to mesh resolution is not negligible but remains of the order of magnitude of experimental uncertainty. Results are reported in figure S19 in SM. A convergence test with respect to time step size was also done for this case using the reference mesh, showing no significant effect of the time step (figure S20 in SM).
4.2. Physical parameters
Solving evolution equations for the conformation tensor and the viscoelastic stress requires setting values for the Zimm relaxation time (
$\tau _{Z}$), the polymer extensibility (
$\xi$), the initial zero-shear polymeric added viscosity (
$\eta _{p0}$) and the initial maximum elongational viscosity (
$\eta _{E0}$). The latter is fitted using pressure losses as we will see in the results section. At dilute concentration, the Zimm time is usually defined by
\begin{equation} \tau_{Z}=\frac{\eta_{s}[\eta]M}{RT}, \end{equation}
where 
$R$ is the ideal gas constant and 
$T$ is the temperature. As already defined after (2.2), the equilibrium polymer extensibility is the ratio of the contour length (stretched polymer length), 
$L$, to the unperturbed root mean square end-to-end distance, 
$R_{EE}$. The polymer contour length is calculated from the polymer structure, in particular, backbone bond lengths and angles, i.e.
\begin{equation} L=0.84\,\ell\frac{M}{M_{1}}, \end{equation}
where 
$\ell =4.4$ Å is the cumulative length of backbone bonds in one PEO monomer, 
$M_{1}=44$ Da is the monomer molecular weight, and where the coefficient 
$0.84$ accounts for length reduction due to typical bond angles. Unperturbed coil sizes are assessed via the intrinsic viscosity with the notion of hydrodynamic volume (Teraoka Reference Teraoka2002),
\begin{equation} R_{EE}=2.1\left(\frac{[\eta]M}{N_{A}}\right)^{{1}/{3}}, \end{equation}
where 
$N_{A}$ is the Avogadro number. Combining with Mark–Houwink–Sakurada law (see (2.10)), we have
\begin{equation} R_{EE}=2.1\left(\frac{KM^{a+1}}{N_{A}}\right)^{{1}/{3}}. \end{equation}
The zero-shear polymeric added viscosity (
$\eta _{p0}$) is experimentally measured from initial solutions, as described in § 3.2. Physical parameters are summarised in table 2.
Table 2. Samples parameters (
$20\,^{\circ }\textrm {C}$).
4.3. Solver implementation
The viscoelastic model ((2.2) and (2.9)) and the mechanochemical model ((2.3) and (2.5)) were implemented using the OpenFOAM (version 6) package (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998), by a modification of the pimpleFoam algorithm. For each time step, the solver processes as follows.
(1) Loop until convergence of the outer loop (convergence criterion on pressure residual).
(i) Solve conformation tensor (2.2).
(ii) Update viscoelastic stress (2.9).
(iii) Update coefficient of the momentum equation.
(iv) Enter the velocity–pressure corrector loop (usually 2 iterations).
(a) Solve pressure equation.
(b) Update velocities.
(v) Update reaction rate (2.5).
(vi) Solve concentration (2.3).
A systematic study of the accuracy and consistency of the schemes used for the present study is outside of our scope, and we rather focus here on the overall stability. Even though Weissenberg numbers in excess of 500 are simulated here, a straightforward implementation of the conformation tensor gives good stability provided that a limited advection scheme is used (we employ the limitedLinear scheme for all variables, where the interpolation is linear with a Sweby limiter). Also, a limited least squares gradient scheme is employed for the computation of 
$\boldsymbol {\nabla }\boldsymbol {u}$, although the limiter is not required in every case. It is not necessary (again, in the scope of stability) to resort to change of variable techniques such as log-conformation formulation. The resulting conformation field is smooth and naturally bounded in regions of high extension since source terms depending on 
$\boldsymbol {\nabla }\boldsymbol {u}$ cancel each other in these regions.
With the version of OpenFOAM used in this study, it is not possible to handle the tensorial viscosity term of the polymeric stress in an implicit manner. Because of this, a chequerboard pattern appears in the pressure field where the polymeric stress dominates the momentum equation. This issue is largely documented for co-located finite volume implementations, and several remedies have been published (Oliveira, Pinho & Pinto Reference Oliveira, Pinho and Pinto1998; Oliveira Reference Oliveira2000; Favero et al. Reference Favero, Secchi, Cardozo and Jasak2010; Matos, Alves & Oliveira Reference Matos, Alves and Oliveira2010; Habla, Obermeier & Hinrichsen Reference Habla, Obermeier and Hinrichsen2013; Fernandes et al. Reference Fernandes, Araujo, Ferrás and Nóbrega2017; Pimenta & Alves Reference Pimenta and Alves2017; Niethammer et al. Reference Niethammer, Marschall, Kunkelmann and Bothe2017). In the present case, the only successful approach turned out to be the both-sides diffusion (BSD) technique (Fernandes et al. Reference Fernandes, Araujo, Ferrás and Nóbrega2017). A viscous term is added to both sides of the momentum equation as follows:
\begin{equation} \rho\left(\frac{\partial\boldsymbol{u}}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}\right)-\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\eta_{s}+\eta^{{\star}}\right)\boldsymbol{\nabla}\boldsymbol{u}={-}\boldsymbol{\nabla} p+\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\boldsymbol{\tau}_{p}-\eta^{{\star}}\boldsymbol{\nabla}\boldsymbol{u}\right). \end{equation}
Here 
$\rho$ is the fluid density and 
$\eta ^{\star }$ an artificial viscosity which will be discussed below. The left-hand side is treated implicitly by the solver and the right-hand side is an explicit source term. The implicit discretisation stencil of the Laplacian is smaller than its explicit counterpart, resulting in the addition of a fourth-order diffusion term which smoothes out high spatial frequency variations of the velocity, and whose action vanishes upon mesh refinement. The choice of the artificial viscosity depends on the constitutive model, however, for the BSD technique to be efficient, 
$\boldsymbol {\tau }_{p}$ and 
$\eta ^{\star }\boldsymbol {\nabla }\boldsymbol {u}$ should be of the same order of magnitude. Many formulations have been tested in the present study and the following gives the most stable result while minimising unwanted diffusivity:
\begin{equation} \eta^{{\star}}=\eta_{E}\lambda\,\boldsymbol{n}_{f}^{\textrm{T}}\boldsymbol{\cdot}\boldsymbol{\mathsf{A}}\boldsymbol{\cdot}\boldsymbol{n}_{f}. \end{equation}
Here 
$\eta ^{\star }$ is computed at cell faces, and 
$\boldsymbol {n}_{f}$ is the face normal. This expression includes three advantageous features: 
$\eta ^{\star }$ depends on 
$\boldsymbol{\mathsf{A}}$, which is, as mentioned above, a smooth and bounded field.
(i) For extended polymers,
$\eta ^{\star }$ allows diffusion of the momentum in the direction of the polymer strands only. In addition, when the face normal and the direction of the polymer strands are aligned, $\boldsymbol{\mathsf{A}}\sim \lambda \boldsymbol {n}_{f}\boldsymbol {n}_{f}$, and, therefore, the traction vector, $\boldsymbol {n}_{f}^{\textrm {T}}\boldsymbol {\cdot }\boldsymbol {\tau }_{p}$, is also normal. Since in that case, $\boldsymbol {n}_{f}^{\textrm {T}}\boldsymbol {\cdot }\boldsymbol {\tau }_{p}\boldsymbol {\cdot }\boldsymbol {n}_{f}\sim \eta ^{\star }\boldsymbol {n}_{f}^{\textrm {T}}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {n}_{f}$, the BSD terms are acting to make the traction vector fully implicit.(ii) In regions where the polymer conformation is near equilibrium,
$\lambda^2$ is very small, and so is $\eta^{\ast}$. Therefore, the BSD correction is switched on only in regions where the chequerboard pattern is likely.
The effect of the BSD terms on simulated degradation and pressure loss is assessed in figure S21 in SM.
The pressure at the inlet and the flux of reacted species at the outlet are monitored, and simulation is advanced until they reach a steady state (or fluctuates around a steady average).
4.4. Quantification of polymer degradation
For a steady constriction flow, we can define a steady-state global degradation, 
$\varPhi$, either by integrating the flux of degraded polymer at the outlet, or by integrating the scission rate over the whole simulation volume,
\begin{equation} \varPhi=\frac{\displaystyle\iint_{{outlet}}\varphi\boldsymbol{u}\boldsymbol{\cdot}\mathrm{d}\boldsymbol{S}}{Q}=\frac{\displaystyle\iiint k\left(1-\varphi\right)~\mathrm{d}\mathcal{V}}{Q}, \end{equation}
where 
$Q$ is the outlet flow rate. The second expression is preferred because it requires less time steps to converge to a steady-state value. At large flow rates, 
$\varPhi$ might be fluctuating because of flow instabilities, and a pseudo-steady-state value is computed by time averaging. This computed value will be compared with the experimental degradation defined in (3.2).
5. Results and discussion
5.1. Experimental results
Polymer degradation of the three molecular weights is shown in figure 3 for each constriction. As expected from literature, the common behaviour of the three geometries is an increase in degradation with increasing flow rate and increasing molecular weight. Degradation indices higher than 100 % can be seen in figure 3(b). They correspond to a situation where the measured zero-shear viscosity is lower than the one predicted by the Mark–Houwink–Sakurada law for a complete halving of polymer chains. In this scenario, it is likely that some of the chains undergo multiple scission events, leading to more than halving of the initial average molecular weight.
Figure 3. Polymer degradation as a function of flow rate. Experiments are run in triplicate: points are averaged values, filled areas show extrema. (a) Fused capillary. (b) Pinhole 25. (c) Pinhole 50.
For most of the data series, error bars are usually of 5 to 10 percentage points. This can be attributed to the consistency limit of the rheometer, assessed from pure water measurements throughout the whole study, as being approximately 0.02 mPa s. It can be considered as a maximum systematic error between series measured on different days. Going back to (3.2), an error of 0.02 mPa s in the viscosity loss would lead to a 15 percentage points error in the degradation index. For PEO 1000k through the capillary geometry, the error at low flow rates can be larger than 20 percentage points. This particularly large span was linked to the presence of initial contaminant particles seen during the first series. A large error is also noted for PEO 300k through the capillary at large flow rates, but the reason is unknown.
The direct consequence of large error bars is that it is not possible to define a critical flow rate (and, therefore, a global critical strain rate) for the onset of degradation, which is the value typically reported in the literature (Islam, Vanapalli & Solomon Reference Islam, Vanapalli and Solomon2004; Vanapalli et al. Reference Vanapalli, Ceccio and Solomon2006). This is less important in this reported work because the focus is to model the level of degradation rather than its onset point. Besides, the critical strain rate used in the model (2.5) has to be viewed as a molecular property, and there is no need a priori to deduce it from an experimental flow rate where the detailed flow gradient is heterogeneous.
As suggested by Nguyen & Kausch (Reference Nguyen and Kausch1991), in order to assess the influence of the geometry of the constriction, the degradation of PEO 1000k in each constriction is plotted against several characteristic variables of the flow (see figure 4): the flow rate, which is the parameter controlled during the experiments (figure 4a); the Reynolds number (figure 4b) defined by
\begin{equation} Re=\frac{\rho Q}{\dfrac{\rm \pi}{4}\eta_{s}D}, \end{equation}
where 
$Q$ is the flow rate, and 
$D$ is the constriction diameter; the average velocity at the constriction, 
$U$ (figure 4c), defined by
\begin{equation} U=\frac{Q}{\dfrac{\rm \pi}{4}D^{2}}; \end{equation}
the nominal strain rate, 
$\dot {\varepsilon }$ (figure 4d), defined by
\begin{equation} \dot{\varepsilon}=\frac{U}{D}=\frac{Q}{\dfrac{\rm \pi}{4}D^{3}}; \end{equation}
and the power loss through the system, 
$\dot {W}$ (figure 4e), defined by
\begin{equation} \dot{W}=Q\Delta p, \end{equation}
where 
$\Delta p$ is the measured pressure loss. Degradation curves for PEO 300k and PEO 600k are similar and reported in SM. The following few remarks can be made from reviewing these plots.
(i) None of these rescaling quantities are able to collapse the degradation curves into a single master curve.
(ii) Degradation as a function of rescaled quantities is higher overall with the pinhole geometries than with the fused capillary. We can explain this difference by noting that the pinholes produce essentially an extensional flow upstream of the constriction, while a shearing component is more present in the flow through the capillary. It is indeed expected that polymer chains are less prone to scission in a shear flow because of their tumbling motion (Odell et al. Reference Odell, Keller and Muller1992).
(iii) For a given flow rate (figure 4a), degradation is higher with the narrower pinhole, but similar with the fused capillary and the larger pinhole. The same observation applies for the Reynolds number (figure 4b) and the power loss (figure 4e). In particular, the suggestion that polymer scission would be essentially governed by the global energy input is not true in our case (Nguyen & Kausch Reference Nguyen and Kausch1991).
(iv) Rescaling by average velocity in the constriction (figure 4c) yields overlapping degradation curves for the two pinholes, but the curve is still lower for the fused capillary.
(v) The nominal strain rate does not recapitulate the degradation phenomenon (figure 4d), even when comparing the two pinholes (this is not true for PEO 300k where the results from the two pinholes overlap).
Figure 4. Degradation of PEO 1000k for each constriction geometry, as a function of various characteristic parameters of the flow. (a) Flow rate (
$\textrm {ml}\,\textrm {h}^{-1}$). (b) Reynolds number. (c) Velocity in constriction (
$\textrm {m}\,\textrm {s}^{-1}$). (d) Nominal strain rate (
$10^6\ \textrm {s}^{-1}$). (e) Power loss (mW).
To conclude, it is clear that the geometry of the constriction and equally the detailed pattern of the flow have a strong influence under the present experimental conditions. Therefore, a detailed CFD approach is necessary to model polymer scission. We now turn our attention to the results of the simulations and a comparison with the reported experiments.
5.2. Model calibration with the fused capillary geometry
Polymer scission is simulated in the fused capillary in order to fit the two parameters of the model, 
$\dot {\varepsilon }_{c}$ and 
$\eta _{E0}$, for each molecular weight. The choice of using the smooth geometry of the capillary for calibration is motivated by the relatively small effect of viscoelasticity on the flow pattern. Indeed, the pressure loss can be accurately described by the simulation when fitting 
$\eta _{E0}$. The result of the fit is given in table 3, and the resulting degradation and pressure loss curves are compared with experiments in figure 5. A 
$\pm 10\,\%$ variation on each parameter is also shown with dotted lines: an increase in 
$\dot {\varepsilon }_{c}$ leads to a lower degradation curve, and an increase in 
$\eta _{E0}$ leads to a higher pressure curve.
Table 3. Results of parameters fit with the fused capillary geometry.
Figure 5. Simulated polymer degradation (a) and pressure loss (b) in the fused capillary constriction, and comparison with experiments. Similar pressure curves are obtained for PEO 600k and PEO 300k. Points show experimental (Exp.) data, filled areas show experimental range, solid lines show simulation (Sim.) results. Dotted lines in the degradation plot show a 
$\pm 10\,\%$ variation of 
$\dot {\varepsilon }_{c}$ with constant 
$\eta _{E0}$; dotted lines in the pressure plot show a 
$\pm 10\,\%$ variation of 
$\eta _{E0}$ with constant 
$\dot {\varepsilon }_{c}$.
It can be seen from figure 5 that after calibration the model is able to describe both polymer degradation and pressure loss as a function of flow rate within experimental uncertainty.
One advantage of detailed CFD modelling is the ability to study the spatial distribution of scission events. The maps of pressure, velocity norm, strain rate 
$\mathsf {D}=\sqrt {\boldsymbol{\mathsf{D}}:\boldsymbol{\mathsf{D}}}$ with 
$\boldsymbol{\mathsf{D}}=(\boldsymbol {\nabla }\boldsymbol {u}+(\boldsymbol {\nabla }\boldsymbol {u})^{\textrm {T}})/2$, polymer square extension (
$\lambda$), scission rate (
$kc/c_{0}$) and degradation (
$\varphi$) are shown in figure 6 for PEO 1000k at 140 
$\textrm {ml}\,\textrm {h}^{-1}$ (flow from left to right). The map of the scission rate shows that scission events occur primarily just before the narrowest part of the constriction and mainly towards the symmetry axis. There is scission also computed downstream around the liquid jet, likely due to the shearing between the jet and the low-velocity surrounding fluid, although the physical interpretation of this effect is unclear. However, this contribution is small (a few %), and absent for flow rates below 80 
$\textrm {ml}\,\textrm {h}^{-1}$. The result of this spread in scission rate is a rather uniform jet where around 80 % of the chains have broken when exiting the constriction.
Figure 6. Simulation results for the flow of PEO 1000k solution at 140 
$\textrm {ml}\,\textrm {h}^{-1}$ through the fused glass capillary geometry. Magnification from full computation mesh. Flow from left to right.
5.3. Simulation results for the pinhole geometry
We now investigate simulation results for the pinhole geometry. The parameters fitted using the capillary geometry are employed for the two pinholes. Figure 7 shows maps of velocity, scission rate and degradation of PEO 1000k at 40 
$\textrm {ml}\,\textrm {h}^{-1}$ through the pinhole 25 (top images), and 150 
$\textrm {ml}\,\textrm {h}^{-1}$ through the pinhole 50 (bottom images). The reason for comparing this combination is that, as reported in figure 4(c), the fluid velocity in the constriction and total degradation should be similar. As shown in figure 7, the velocity field is characterised by an upstream lip vortex as expected of viscoelastic flows through abrupt constrictions (Boger Reference Boger1987). As in the case of the fused capillary discussed above, the maximum scission rate occurs mainly just before entering the constriction, but the effect of the vortices is to stretch the scission zone far (many constriction diameters) upstream. This contrasts with our Newtonian simulations reported previously where scission occurs mainly close to the sharp entrance edge (Rognin et al. Reference Rognin, Willis-Fox, Aljohani and Daly2018). There is also a scission zone downstream at the boundary of the jet, but it accounts likewise for less than a few % of the total degradation. Once again, the degradation is rather uniform in the jet and similar for both geometries at their respective flow rate.
Figure 7. Simulation results for the flow of PEO 1000k solution at 30 
$\textrm {ml}\,\textrm {h}^{-1}$ through the pinhole 25 constriction (a) and at 150 
$\textrm {ml}\,\textrm {h}^{-1}$ through the pinhole 50 (b,c). Magnification from full computation mesh. Flow from left to right.
Simulations at various flow rates are compared with experimental degradation curves in figures 8 and 9 (pinhole 25 and pinhole 50, respectively). Regarding pinhole 25 (figure 8) and PEO 1000k, the model accurately predicts polymer degradation at high flow rate (30 
$\textrm {ml}\,\textrm {h}^{-1}$ and above), but there is a non-negligible discrepancy at low flow rates (20 
$\textrm {ml}\,\textrm {h}^{-1}$ and below). Note that by construction the model is unable to describe experimental degradation over 100 % (multiple chain scission). Simulated values above 100 % are due to degradation downstream where the fully degraded jet mixes with fresh solution and part of this fresh solution is also degraded in the mixing velocity gradient, hence a sum larger than 100 % degradation with respect to the inlet flow. In practice, the fresh solution would be flushed by fully degraded fluid over time. Because the error is small (a few %), this behaviour is not simulated to spare computation time. Regarding PEO 600k, the discrepancy is large for the whole range of flow rates. Nevertheless, experimental and simulated curves still intersect in the vicinity of 50 % degradation, showing that the model is acceptable with respect to orders of magnitudes. Finally, regarding PEO 300k, the simulated degradation is certainly inconsistent, and this issue will be investigated below.
Figure 8. Simulated polymer degradation (a) and pressure loss (b) in the pinhole 25 constriction, and comparison with experiments.
Figure 9. Simulated polymer degradation (a) and pressure loss (b) in the pinhole 50 constriction, and comparison with experiments.
In addition, simulated pressure loss shown in figure 8 for the case of PEO 1000k exceeds experimental values, especially at highest flow rates. A similar trend is observed for the other molecular weights.
Regarding pinhole 50 (figure 9), degradation of PEO 1000k is slightly underestimated, while it is accurately predicted for PEO 600k. Once again, simulated degradation is inconsistent for PEO 300k. Simulated pressure loss shows this time an underestimation with respect to experimental measurements, contrasting with the simulated pressure in pinhole 25. This suggests that the issue is not in fitting 
$\eta _{E0}$, but rather lies in the stress model as a whole, which is not able to accurately describe the flow in this geometry. We will come back to this point below where we discuss the limits of the model.
To summarise, the process of fitting the critical strain rate using the smooth fused capillary geometry resulted in an accurate prediction of the degradation in the two pinholes in the cases of PEO 1000k and PEO 600k. This can be seen by plotting simulation results by molecular weight instead of by geometry, as shown in figures 10(a) and 10(b). This process failed for PEO 300k, and it could be explained by the fact that degradation in the capillary was too little (less than 20 % at highest flow rates) to provide a sound base for parameters fitting. If instead we use pinhole 25 to fit the critical strain rate for PEO 300k, the new value of 
$2.7\times 10^{6}$ s
$^{-1}$ is found, and degradation can be reasonably predicted in the pinhole 50 (figure 10c).
Figure 10. Simulated polymer degradation compared with experiments. A new fit of 
$\dot {\varepsilon }_{c}$ is used for PEO 300k. (a) PEO 1000k. (b) PEO 600k. (c) PEO 300k.
5.4. Analysis of fitted parameters
Critical strain rates (PEO 300k: value fitted using pinhole 25) are plotted against molecular weights in figure 11. In the present model, 
$\dot {\varepsilon }_{c}$ is a molecular parameter and should be distinguished from global nominal strain rates such as defined by (5.3) which are typically reported in the literature (see Garrepally et al. (Reference Garrepally, Jouenne, Olmsted and Lequeux2020) for a review of scaling laws with respect to the nominal strain rate). If the dilute free-draining molecular theory holds here, a 
$-2$ exponent is expected. A 
$-1.52\pm 0.52$ exponent if found in figure 11, which suggests that, although the present values could be consistent with a molecular view, the error is too large to provide a conclusive scaling law and a larger range of molecular weights would be needed.
Figure 11. Critical strain rates vs molecular weight.
Lastly, fitted values for the maximum elongational viscosity (
$\eta _{E0}$) are more difficult to interpret (see table 3). Values are only weakly dependent on molecular weight. From a dilute theory, we would expect a scaling of 
$\propto \eta _{p0}\xi ^{2}$ and because the added shear viscosity 
$\eta _{p0}$ is approximately the same for all samples, the scaling should be as 
$\propto \xi ^{2}\sim M^{1.1}$. The data are too scattered to be conclusive. Also, we note that even if present solutions are dilute with respect to zero-shear viscosity, this property can be questioned for unravelled polymers (Clasen et al. Reference Clasen, Plog, Kulicke, Owens, Macosko, Scriven, Verani and McKinley2006; Prabhakar et al. Reference Prabhakar, Sasmal, Nguyen, Sridhar and Prakash2017). Therefore, concentration effects can be expected to influence the scaling of the elongational viscosity. Finally, the stress model can account for the pressure loss accurately in the smooth capillary geometry, but only qualitatively in sharp constrictions. This suggests that the parameter 
$\eta _{E0}$ is an ad-hoc parameter that could be related only qualitatively to a molecular property.
5.5. Comparison with other scission models
We now compare the model of scission rate used in this study (2.5) to the two other forms of first-order rates mentioned in § 2.1. The first form is a linear function of the strain rate, still retaining a cut-off at a critical strain rate, 
$\dot {\varepsilon }_{c\,{lin}}$,
\begin{equation} k_{{linear}}=\begin{cases} 0 & \text{if }\left\lfloor \boldsymbol{\nabla}\boldsymbol{u}:\boldsymbol{\mathsf{A}} \right\rfloor /\lambda<\dot{\varepsilon}_{c\,{lin}},\\ \frac{1}{2}\left\lfloor \boldsymbol{\nabla}\boldsymbol{u}:\boldsymbol{\mathsf{A}}\right\rfloor & \text{otherwise}. \end{cases}\end{equation}The second form is a TABS model,
\begin{equation} k_{{TABS}}=k_{0}\exp\left(\frac{\left\lfloor \boldsymbol{\nabla}\boldsymbol{u}:\boldsymbol{\mathsf{A}}\right\rfloor }{k_{1}}\right). \end{equation}
Because we observed no significant change in viscosity over several months for samples stored at room temperature, the no-flow degradation rate can be assessed to be lower than 
$10^{-7}$ s
$^{-1}$. The thresholding effect of the TABS model described in § 2.1 implies that an exact value for 
$k_{0}$ is not required as long as it is much lower than the inverse residence time in the flow, which here is larger than the reciprocal second. Thus, setting 
$k_{0}=10^{-7}\,\textrm {s}^{-1}$, only 
$k_{1}$ needs to be fitted.
Parameters 
$\dot {\varepsilon }_{c\,{lin}}$ and 
$k_{1}$ are therefore fitted using polymer degradation of PEO 1000k flowing through the capillary geometry. Results are reported in figure 12. The linear model can account for the degradation relatively well, although not as accurately as the model used in this study. On the other hand, the TABS model produces a sharp transition to degraded polymer above 80 
$\textrm {ml}\,\textrm {h}^{-1}$, which does not reflect experimental data. This shows that the model used in this study, which is based on the unravelling dynamics of polymer chains, is more appropriate to predict scission of long and flexible molecules in a constriction flow.
Figure 12. Simulated polymer degradation of PEO 1000k in the capillary geometry, comparing different scission models. Fitted parameters: this paper (2.5): 
$\dot {\varepsilon }_{c}=4.5\times 10^{4}$ s
$^{-1}$, linear model: 
$\dot {\varepsilon }_{c\,{lin}}=9.0\times 10^{3}$ s
$^{-1}$, TABS model: 
$k_{1}=2.7\times 10^{3}$ s
$^{-1}$.
5.6. Limits of the model and possible improvements
Although the model provides an accurate prediction of polymer scission for a large range of experimental values, it gives only a reasonable order of magnitude in some cases. We now investigate the limits of the model and suggest some possible improvements.
The first limit is probably given by the very large extensibility of the molecules. The molecular weights investigated here are molecules that are several microns long when fully stretched. With respect to the several tens of microns for the characteristic dimensions of the constriction, the scission rate might cease to be a local variable, but could depend on the velocity field over an extent of several microns. One option is to allow the scission rate to depend on higher-order derivatives of the velocity field, similar to drag models of immersed objects.
Another challenge relevant to both scission and stress models is how to provide an accurate polymer physics in mixed flows, i.e. in flows that are neither purely Couette shear nor purely extensional. A broad literature exists on these two types of standard flows, but only a few studies focus on stretched polymers in mixed-type flows (Jain et al. Reference Jain, Sasmal, Hartkamp, Todd and Prakash2015; Prakash Reference Prakash2019). The models used here are relevant to purely extensional flow but could be improved to be accurate in a broader range of flow fields. This would require additional integration of molecular dynamics studies into continuum scale constitutive equations.
To further emphasise the importance of an accurate stress model, degradation of PEO 1000k was simulated using the critical strain rate of 
$4.5\times 10^{4}$ s
$^{-1}$, but assuming a Newtonian stress (limit of the infinite dilution). Results are shown in figure 13. The change is small and within experimental uncertainty in the case of the capillary (smooth constriction). By contrast, the shift is dramatic in the case of the sharp pinhole 50, and can be related to the absence of lip vortices at the entrance of the constriction in Newtonian simulations.
Figure 13. Simulated polymer degradation of PEO 1000k, comparing a viscoelastic to a Newtonian stress model.
Although concentration effects are already visible through the presence of a viscoelastic stress, intermolecular interactions are not accounted for in this model. Yet it is clear that even at dilution below the overlap concentration, polymer chains that are unravelled far from equilibrium are likely to interact. This not only affects the macroscopic stress as discussed above, but presumably influences the growth dynamics of molecular internal tension and eventually bond rupture kinetics. Experimental studies of contraction flows have shown that an increase in polymer concentration leads to a decrease in the global critical strain rate (Nghe, Tabeling & Ajdari Reference Nghe, Tabeling and Ajdari2010). Yet it is challenging to decouple the effect of intermolecular interactions from that of a changing flow pattern due to increased elastic stress. Finite concentration molecular dynamics could help us understand better the role of concentration in degradation kinetics.
Another limit of the present model is that it can account for only one scission event per chain, even though multiple scissions can occur at highest strain rates in experimental settings. This could be simulated by solving an additional concentration field (first scission products) and solving scission for this lower molecular weight. More generally, a highly polydispersed polymer population could be accounted for by using a discrete binning of the molecular weights distribution, each bin having its set of physico-chemical properties and evolution equations (Sorbie & Roberts Reference Sorbie and Roberts1984). Another approach following Pereira et al. (Reference Pereira, Mompean and Soares2018) could be to simulate the evolution of the local average molecular weight, allowing its decrease to a lower value than half of the original polymer. In both cases, a scaling law for the critical strain rate with respect to molecular weight would have to be assumed beforehand.
Finally, the model has been validated in laminar flow conditions with Reynolds numbers below 2000 (figure 4b). Nevertheless, the constitutive equations are not bound to any laminar assumption, and there is no presumption that the model would fail in turbulent conditions provided the flow is fully resolved, for example, in DNS simulations. This point applies to unsteady flows in general, and more work would be valuable to study the model in transient flows such as those encountered in inkjet, spray or sputtering.
6. Conclusion
A continuum model for the flow-induced degradation kinetics of flexible polymers was presented. The model was implemented in a finite volume CFD software and tested against degradation experiments of dilute PEO solutions flowing through narrow constrictions. Alongside typical near-equilibrium properties such as polymer extensibility, relaxation time and zero-shear polymeric viscosity, the model requires two far-from-equilibrium parameters: a critical strain rate above which chain scission can occur and a maximum extensional viscosity. The approach followed in this study was to calibrate the model by fitting those two parameters using experimental degradation data obtained from a smooth constriction flow. The model was then tested against experiments using two sharp constrictions of different diameters. The approach provided an accurate prediction of the polymer degradation, except for the lowest molecular weight for which only little degradation had been observed in the first instance in the smooth capillary. A recalibration for this low molecular weight using one sharp constriction led to a better fit for the other constriction.
The model could be used to study the influence of flow and process designs on degradation in fields where preservation of macromolecules is a concern. On the other hand, because this degradation kinetics model is based on mechanical tension in polymer chains, it is also relevant to the activation of mechanochemical compounds in fluid flows. It should provide an efficient tool to design more adequate flow systems for these novel flow-activated materials.
Supplementary material
Supplementary material and source code are available at https://doi.org/10.1017/jfm.2021.646.
Funding
This work was partially funded by the King Abdulaziz City for Science and Technology (KACST), and the Engineering and Physical Sciences Research Council (EPSRC) (grant no. EP/S009000/1).
Declaration of interests
The authors report no conflict of interest.
Author contributions
E.R. conceived the model, developed the software and carried out the simulations. E.R., N.W-F. and T.Z. carried out the experiments. E.R. wrote the original draft. T.A. and R.D. supervised the project. All authors reviewed and edited the manuscript.
Open access
