Analysis and modeling of localized invariant solutions in pipe flow
Paul Ritter, Stefan Zammert, Bruno Eckhardt, Marc Avila

TL;DR
This paper investigates the spatial structure of localized turbulence in pipe flow, modeling decay rates and velocities of turbulent fronts using linearized Navier-Stokes equations, and validates the model against turbulent puffs.
Contribution
It introduces a linearized modeling approach for the decay and velocity profiles of localized turbulent structures in pipe flow, linking kinematics with spatial structure.
Findings
Exponential decay of turbulent fronts towards laminar flow
Model accurately predicts decay rates and front velocities
Validates model with turbulent puff data
Abstract
Turbulent spots surrounded by laminar flow are a landmark of transitional shear flows, but the dependence of their kinematic properties on spatial structure is poorly understood. We here investigate this dependence in pipe flow for Reynolds numbers between 1500 and 5000. We compute spatially localized relative periodic orbits in long pipes and show that their upstream and downstream fronts decay exponentially towards the laminar profile. This allows to model the fronts by employing the linearized Navier-Stokes equations, and the resulting model yields the spatial decay rate and the front velocity profiles of the periodic orbits as a function of Reynolds number, azimuthal wave number and propagation speed. In addition, when applied to a localized turbulent puff, the model is shown to accurately approximate the spatial decay rate of its upstream and downstream tails. Our study provides…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Analysis and modeling of localized invariant solutions in
pipe flow
Paul Ritter
Center of Applied Space Technology and Microgravity (ZARM), University of Bremen, 28359 Bremen, Germany
Institute of Fluid Mechanics, Friedrich-Alexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany
Stefan Zammert
Laboratory for Aero- and Hydrodynamics, TU Delft, 2682 Delft, The Netherlands
Bruno Eckhardt
Fachbereich Physik, Philipps-Universität Marburg, 35032 Marburg, Germany
JM Burgerscentrum, TU Delft, 2682 Delft, The Netherlands
Marc Avila
Center of Applied Space Technology and Microgravity (ZARM), University of Bremen, 28359 Bremen, Germany
Institute of Fluid Mechanics, Friedrich-Alexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany
Abstract
Turbulent spots surrounded by laminar flow are a landmark of transitional shear flows, but the dependence of their kinematic properties on spatial structure is poorly understood. We here investigate this dependence in pipe flow for Reynolds numbers between 1500 and 5000. We compute spatially localized relative periodic orbits in long pipes and show that their upstream and downstream fronts decay exponentially towards the laminar profile. This allows to model the fronts by employing the linearized Navier–Stokes equations, and the resulting model yields the spatial decay rate and the front velocity profiles of the periodic orbits as a function of Reynolds number, azimuthal wave number and propagation speed. In addition, when applied to a localized turbulent puff, the model is shown to accurately approximate the spatial decay rate of its upstream and downstream tails. Our study provides insight into the relationship between the kinematics and spatial structure of localized turbulence and more generally into the physics of localization.
I Introduction
Due to its stochastic and fluctuating nature, the classical approach towards understanding turbulent fluids has been a statistical one, which dates back to Osbourne Reynolds Reynolds (1895). In recent years an alternative approach has emerged, in which the (discretised) Navier-Stokes equations are viewed as a high-dimensional dynamical system and the tools of bifurcation and chaos theory are applied to describe turbulent motions Kerswell (2005); Eckhardt et al. (2007); Eckhardt (2008). The key idea of this approach is that the turbulent dynamics is shaped by simple exact invariant solutions to the governing equations such as traveling waves Faisst and Eckhardt (2003); Wedin and Kerswell (2004) and relative periodic orbits Duguet et al. (2008); Willis et al. (2013) in pipe flow. The dynamically most relevant solutions have a relatively small number of unstable directions so that a generic turbulent trajectory spends a significant amount of time in their vicinity Gibson et al. (2008); Kawahara et al. (2012); Schneider et al. (2007). Turbulent trajectories depart from the vicinity of the solutions along their unstable manifolds, which can subsequently govern the flow evolution for a considerable period of time Suri et al. (2017). In principle, all properties of the turbulent flow can be derived by an appropriate weighted average over the fundamental solutions Cvitanović (1988); Cvitanović and Eckhardt (1989); Cvitanovic and Eckhardt (1991); Cvitanović (2013); Kreilos and Eckhardt (2012), but deploying this approach is extremely challenging even at low (transitional) Reynolds numbers Schneider et al. (2007); Willis et al. (2016); Kreilos and Eckhardt (2012).
Transitional shear flows are characterized by localized chaotic spots surrounded by laminar flow Reynolds (1883); Emmons (1951); Tillmark and Alfredsson (1992); Lemoult et al. (2013). These spots already contain all the salient features of fully turbulent flow Wygnanski and Champagne (1973); Barkley et al. (2015); Cerbus et al. (2017) and hence pose an ideal prototype for a bottom-up study of turbulent dynamics. Because of prevalence of intermittency at the onset of turbulent shear flow Rotta (1956); Moxey and Barkley (2010); Avila et al. (2011), spatially localized invariant solutions are indispensable for its successful description as a dynamical system. The first such solutions were discovered in plane Couette flow by Schneider et al. Schneider et al. (2010), who computed spanwise-localized equilibria and traveling waves in wide but streamwise short domains.
The first streamwise-localized simple invariant solution was found by Avila et al. in pipe flow Avila et al. (2013). It is a relative periodic orbit with reflectional and two-fold rotational symmetry appearing at a saddle-node bifurcation. In this symmetry subspace the lower branch solution (shown in fig. 1a and referred to as LB2 in the following) has a single unstable direction, whereas the upper branch solution is stable close to the saddle-node bifurcation. As the Reynolds number increases, a bifurcation cascade culminating at a boundary crisis gives rise to transient chaotic dynamics Avila et al. (2013); Ritter et al. (2016), and subsequent changes in the phase-space progressively enhance the spatio-temporal complexity of the flow Mellibovsky et al. (2009); Ritter et al. (2016). The bifurcations of the coherent structures in pipe flow follow the same pattern as observed in small computational cells in plane Couette Kreilos and Eckhardt (2012); Kreilos et al. (2014) and plane Poiseuille flows Zammert and Eckhardt (2015). The finding of spanwise- and streamwise-localized solutions in both flows Zammert and Eckhardt (2014); Gibson and Brand (2014) suggests that similar scenarios may occur in spatially extended domains.
Gibson and Brand Gibson and Brand (2014) observed that the amplitude of spanwise localized equilibria in Couette flow decays exponentially far enough from their energetic core. Hence they proposed to model their spatial decay by using the linearized Navier–Stokes equations and solving the arising eigenvalue problem. Interestingly, their model gave with high accuracy the observed spatial decay rates and their dependence on streamwise wave number and the Reynolds number. A simplified version of their model was shown by the same authors to accurately reproduce the spatial decay of their doubly-localized solutions in the streamwise direction Brand and Gibson (2014).
Recently, Zammert and Eckhardt Zammert and Eckhardt (2016) and Barnett et al. Barnett et al. (2016) extended this approach to streamwise relative periodic orbits in channel flow, where the decay rates also depend on their group velocity, in addition to the spanwise wave number and the Reynolds number. In this paper, we use a similar approach to pipe flow. For this purpose, we compute LB2 and its three-fold cousin Chantry et al. (2014) (LB3, see fig. 1b) for a wide range of Reynolds numbers. As in channel flows, we find that the tails of the states decay exponentially in the streamwise direction and that the decay rates can be deduced from the linearized Navier–Stokes equations. However, pipe flow does not permit the simplified modeling approaches used by Brand and Gibson Brand and Gibson (2014), and Zammert and Eckhardt Zammert and Eckhardt (2016). Furthermore, we extend the analysis to the upstream and downstream tails of localized turbulent puffs and verify that the correct decay rates are obtained for such chaotically evolving states as well.
II Numerical method
We consider the incompressible, isothermal flow of a fluid with constant density and kinematic viscosity in a cylindrical pipe of radius driven at a constant average speed . This flow is governed by the Navier-Stokes equations (NSE):
[TABLE]
where is the pressure and is the fluid velocity field in cylindrical coordinates. It satisfies the no-slip boundary condition at the pipe wall and periodic boundary conditions in the azimuthal and axial directions. The length of the computational domain was chosen sufficiently large in order to avoid interaction of the two fronts via the axial periodicity. All results shown in this paper were obtained in pipes of (LB2) and (LB3) in length.
The Hagen-Poiseuille profile is the steady, parabolic laminar solution and reads (subscript “b” is for base flow):
[TABLE]
where is the maximum velocity at the centerline, and denotes the axial unit vector. To facilitate both numerical and theoretical treatment, the NSE are rendered dimensionless by using , and as reference scales for the velocity, pressure and length, respectively. As a consequence, the dimensionless NSE are identical to eq. (1) but setting and replacing the viscosity with the inverse of the Reynolds number , which is the sole control parameter of the problem. The velocity and pressure gradient of the dimensionless laminar flow then take the form and , respectively. Throughout the paper the velocity disturbance is used to visualize the structures.
The direct numerical simulations of the Navier–Stokes equations (1) have been carried out using openpipeflow.org Willis (2016), a hybrid spectral finite-difference Navier-Stokes solver, which uses primitive variables and a PPE-formulation with correct pressure boundary conditions via the influence-matrix method Rempfer (2006); Guseva et al. (2016). In order to compute the localized structures, a two-step approach was employed. First, the edge-tracking technique Itano and Toh (2001); Skufca et al. (2006) was used to bracket the relative periodic orbits to a sufficient degree so that it could be converged in a second step with a Newton–Krylov-hookstep algorithm Viswanath (2007) to relative error . The necessary spatial resolution of the periodic directions depends on the enforced rotational symmetry. We used an axial resolution of Fourier modes for a pipe of length in case of two-fold symmetry and the same amount of modes for a pipe in the three-fold case. The spanwise resolution was () Fourier modes for LB3 (LB2) capturing up to 36th (32nd) wave number (of which a third/half has the same amplitude due to symmetry). The code uses the -rule for dealiasing (i.e. padding) resulting in a physical grid which has three times as many points as there are wave numbers (). The radial direction has been discretized with a minimum of and maximum of finite-difference nodes depending on Re and pipe length. The time step was fixed at a value of .
III Exponential localisation of solutions
The approach of the fields to the asymptotic parabolic flow is best visualized by the deviation from the Hagen-Poiseuille profile, since they have to decay to zero. The isosurfaces of streamwise velocity deviation from laminar flow shown in fig. 1 illustrate the spatial arrangement of streaks of the spatially localized relative periodic orbits LB2 and LB3 at . Far from the active core, all three velocity components decay quickly with respect to the streamwise direction . The semilogarithmic representation in fig. 2 shows that the decay is predominantly exponential, with approximately two orders of magnitude larger than the cross-stream velocities and , thus dominating the decay toward laminar flow. Interestingly, the decay rate of the azimuthal velocity at the downstream tail differs from that of the other two components. To shed light on the origin of this difference, isosurfaces of all three velocity components are shown in fig. 3 for LB3. In the upstream tail all three velocity components feature a predominant sixfold rotational symmetry, whereas in the downstream tail and are predominantly axisymmetric and features a threefold symmetric structure. LB2 exhibits the same features but with fourfold and twofold symmetry, instead of sixfold and threefold, respectively, and hence it is not shown here.
The length of the core of LB2 and LB3 remains nearly constant, whereas its amplitude decreases as Re increases (see fig. 4). This is not surprising because LB2 and LB3 are edge states and can thus be seen as minimal seeds to trigger turbulence Pringle et al. (2012). The decay rate of both their upstream and downstream tails decreases with Re, i.e. the axial velocity profile gradually “opens up”. There is a marked asymmetry, however. While the decay rates of the downstream tails change little with Re, the decay rates of the upstream tails decrease rapidly with increasing Re. This behaviour is similar to that for the relative periodic orbits in plane Poiseuille flow Zammert and Eckhardt (2016); Barnett et al. (2016). Overall, the localization becomes weaker as Re increases.
IV Linear model of spatial decay
IV.1 Mathematical formulation of the model
The exponential decay observed at the tails suggests that these can be modelled with the linearised Navier-Stokes equations (LNSE). Following Gibson & Brand Brand and Gibson (2014), we look for normal mode solutions of the form
[TABLE]
where is the azimuthal wave number dominating at the tail, the spatial decay rate at the tail and the group velocity at which the localized solution (wave packet) travels in the axial direction. The latter is not to be confused with the phase velocity of individual waves within the solution (envelope). Note that is generally complex in a spatial setting. Its real part describes the spatial attenuation (decay rate) and its imaginary part the spatial modulation of the localized solution fronts. Moreover, note that equation (4) describes the tails in a reference frame moving with the group velocity . Although strictly speaking this equation is only valid for relative equilibria (see Barnett et al. (2016)), the temporal variation is negligible at the tails of our relative periodic orbits. Hence equation (4) is used here to model the spatial decay rates far away from their core.
Inserting ansatz (4) into the dimensionless LNSE gives:
[TABLE]
where
[TABLE]
with
[TABLE]
Rearranging (5) with respect to one obtains the following quadratic eigenvalue problem (EVP):
[TABLE]
where
[TABLE]
[TABLE]
and
[TABLE]
with
[TABLE]
The radial derivatives were discretized with a spectral method at Chebyshev collocation points. In order to reduce clustering of grid points near the origin (where the solution is smoother), the differentiation matrices were computed over the interval [-1,1] using points. The derivatives on (0,1] were obtained by “quotienting” out the symmetry of the 2-to-1 map from () to in this representation using the appropriate parities, respectively Trefethen (2000).
The quadratic EVP (6) can be linearized analogous to the reduction of a second-order ODE to first-order, namely by replacing it with a linear system with twice as many unknowns and equations Tisseur and Meerbergen (2001). Here, we choose the so-called “first companion form” (see Tisseur and Meerbergen (2001)) by making the substitution . This yields the generalized eigenvalue problem
[TABLE]
which is subsequently solved with QZ-factorization (generalized Schur decomposition, see Moler and Stewart (1973)).
IV.2 Model results
For given Reynolds number Re, azimuthal wave number and group velocity , positive (negative) eigenvalues of EVP (6) give an approximation for the decay rate at the upstream (downstream) tail of a localized solution. The associated eigenvectors approximate the velocity profiles at the tails.
The model predictions for the tails are computed as follows. First, the Reynolds number is fixed and the group velocity is determined from the DNS. Fig. 5 shows the evolution of as a function of Reynolds number. Close to their saddle-node bifurcation points, LB2 and LB3 travel slightly faster than the mean flow speed 0.5\text{,}\mathrm{\mathit{U}_{\mathrm{cl}}}$$ and the differences grow slowly as Re increases.
Second, the azimuthal wave number is determined from the velocity profiles. As shown in fig. 3, the upstream tail of LB3 (LB2) is dominated by a () rotational symmetry for all components, whereas their downstream tails are predominantly axisymmetric for and , but feature () in . Thus, one needs to consider both and () separately in the model for the downstream tail. The case actually decouples the azimuthal velocity from the other equations in EVP (6) leaving only the diagonal block of the matrices . Since their sum is not singular, is obtained, which is not observed (fig. 3c) but consistent with the fact that decays with a different rate than predicted by the axisymmetric mode for and . Moreover, the mean azimuthal velocity (which is the mode) has to vanish because our localized solutions are reflection symmetric, which precludes a mean rotation.
Figure 6 compares the decay rate of the streamwise velocity disturbance obtained by exponential fits to at the tails of the solutions (square/circular markers), to the prediction based on the LNSE (dashed lines). The agreement is excellent, which confirms the validity of the model.
The radial velocity profiles obtained from the model are compared to the DNS data in fig. 7 for LB3. The agreement of the upstream eigenvectors (model ) with DNS is very good. At the downstream tail, the axisymmetric () model result for and match the DNS result very well, too. To obtain a model prediction for the azimuthal velocity, we solve the LNSE with as suggested by fig. 3c. This yields -0.195\text{,}\mathrm{\mathit{R}}^{-1}$$, whereas the decay rate obtained from DNS is in both cases. However, the magnitude of is very small. Note also that its decay is modulated in space (see figs. 2b, 3c) and this is correctly predicted by the model with . In all other cases, the imaginary part of is zero, consistent with the absence of modulations in the spatial decay. The results for LB2 are very similar to those of LB3 and hence not shown here.
IV.3 Contribution of terms to the spatial decay at tails
The decay rate of the streamwise tails of certain localized solutions in Couette Brand and Gibson (2014) and channel Zammert and Eckhardt (2016) flows can be accurately modeled with a single equation for the streamwise velocity component of the disturbance. These authors compare the contributions of all terms of the streamwise momentum conservation equation at the solution tails, and find that three terms dominate: linear advection of the disturbance by the basic laminar flow, and diffusion of momentum in the spanwise and wall-normal directions. The model resulting from consideration of only these three terms is an advection-diffusion equation, which in pipe flow takes the following form
[TABLE]
We computed the decay rates from (8) and found that they disagree with those of the DNS and full LNSE, as also observed for spatially localized modulated Tollmien-Schlichting waves in channel flow Barnett et al. (2016). However, the upstream rates and eigenvectors computed from the advection-diffusion model are at least comparable to the full simulation (fig. 6), whereas the downstream rates are utterly false (hence not shown in fig. 6). We assessed the reasons underlying the failure of (8), by analyzing individual contributions of all terms in the streamwise momentum conservation equation for LB3 (see Fig. 8, the relative contributions are similar for LB2 and hence not shown here). The diffusive term (yellow) is much smaller than the in-plane contributions to and the other terms. In fact, setting and solving the resulting linear eigenvalue problem does not have any effect on the decay rates, confirming that axial diffusion can be neglected. The pressure gradient (violet) in the upstream and downstream tails approaches its small and constant value in the surrounding base flow. This value is a consequence of the periodic boundary conditions and further decreases in longer pipes and with higher Re.
The relative contribution of the lift-up term (magenta) was found to be larger here than for the solutions of Zammert and Eckhardt Zammert and Eckhardt (2016) in channel and Brand and Gibson Brand and Gibson (2014) in Couette flows. In our solutions, the lift-up term is of similar magnitude as the in-plane diffusion. This suggests that the absence of the lift-up term in the simple advection-diffusion model (8) is responsible for its failure. We gauged the role of the lift-up term by solving EVP (6) without it (i.e. by omitting in ). The corresponding upstream decay rates are nearly identical to those from the advection-diffusion equation, whereas the downstream rates deviate strongly from the DNS (to a similar degree as the advection-diffusion equation).
These findings indicate that the lift-up term plays a key role in the decay of the tails in the pipe. For certain solutions of Couette and channel flows, the main coupling of the streamwise velocity with the other components is via the mass-conservation equation resulting in a very small wall-normal velocity, which does not influence the axial decay rate (see the vanishing value of in the tails in fig. a of Zammert and Eckhardt (2016)). In pipe flow, however, one cannot neglect in the tails and the radial momentum equation is strongly coupled to the axial one via . Hence it is not possible to formulate an accurate single-equation model for the decay of the streamwise velocity at the tails of LB2 and LB3, exactly as for spatially localized modulated Tollmien-Schlichting waves in channel flow Barnett et al. (2016).
IV.4 Application to a turbulent puff
Models based on LNSE have been so far applied to describe the tails of exact coherent solutions. Mellibovsky et al. Mellibovsky et al. (2009) showed that the tails of an edge state and a turbulent puff at decay exponentially, suggesting that LNSE may correctly describe their spatial decay. We here applied the LNSE to a turbulent puff at (see fig. 9), which propagates at exactly the mean speed Avila et al. (2011). Using this value of yields 0.13\text{,}\mathrm{\mathit{R}}^{-1} and $\mu=$-0.0395\text{\,}\mathrm{\mathit{R}}^{-1}, for the upstream () and downstream () front, respectively, which is in excellent agreement with the results from exponential fits to at the tails of the puff. Similarly the velocity profiles from the model match those of the DNS at the tails as shown in fig. 10.
V Discussion
Localized exact coherent structures in pipe flow exhibit exponential localization far away from their active core. This allows for accurate models of the decay based on the LNSE, as in the cases of Couette Gibson and Brand (2014); Brand and Gibson (2014) and channel flow Zammert and Eckhardt (2016); Barnett et al. (2016). The solution of the resulting spatial eigenvalue problem yields two decay rates of different sign for the velocity disturbance at the upstream () and downstream tails (), as a function of Reynolds number Re, azimuthal wave number and group velocity .
The localized solutions investigated here are relative periodic orbits and have either two- or three-fold rotational symmetry and are reflectional symmetric. Their upstream tails feature four/six streaks and all three components of the velocity disturbance decay at the rates as in the model for . Their downstream tails are predominantly axisymmetric and consist of large scale meridional circulation , whereas presents a two-/three-fold symmetry, is much smaller and decays at a much faster rate. The model accurately predicts the decay rates for using at the downstream tails, but the decay of is only qualitatively recovered even when are used in the model.
The decay rate of some solutions in Couette Brand and Gibson (2014) and channel Zammert and Eckhardt (2016) flows can be accurately modeled with a single advection- diffusion equation for the streamwise velocity disturbance. Interestingly, in this equation the Reynolds number and decay rate appear only through the combination , so that one may expect to decrease as . However, the propagation speed of the structures was found to increase with Re in channel flow leading to a faster decrease of with Re upstream and a slower, nearly constant increase downstream. For spatially localized modulated Tollmien-Schlichting waves in channel flow Barnett et al. (2016), and for our solutions, a simple advection-diffusion equation cannot reproduce the results from the full LNSE. We here showed that the lift-up term is significant enough so that it cannot be neglected. This term couples the radial and axial momentum equations and so it is no longer possible to retain a single equation model for the decay rate of the streamwise velocity disturbance. Despite the key role of the lift-up term, we still found that the scaling of the decay rates with Re corresponds to what one would expect from the advection-diffusion equation (8) once the dependence on is taken into account, similar to channel flow Zammert and Eckhardt (2016); Barnett et al. (2016).
The LNSE were also applied to data from DNS of localized turbulence. Here a turbulent puff at was analyzed. The spatial decay rates obtained from the model with , and for the downstream and upstream tails, were found to be in excellent agreement with the DNS data. This confirms the validity of the model for real turbulent patches and emphasizes the crucial interdependence between propagation speed of a localized structure and the spatial decay rate at its tails. Note that the dependence of on the group velocity suggests that in the transition from localized puffs to expanding slugs the localization rate must change accordingly.
Acknowledgements.
Support from the Deutsche Forschungsgemeinschaft (DFG) through grant FOR 1182 and computing time from the “Regionales Rechenzentrum Erlangen (RRZE)” are acknowledged. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915. S. Z. acknowledges financial support by Stichting FOM/NWO-I.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Reynolds (1895) O. Reynolds, “On the dynamical theory of incompressible viscous fluids and the determination of the criterion,” Phil. Trans. Roy. Soc. Lond. A 186 , 123–164 (1895) .
- 2Kerswell (2005) R. R. Kerswell, “Recent progress in understanding the transition to turbulence in a pipe,” Nonlinearity 18 , R 17 (2005) . · doi ↗
- 3Eckhardt et al. (2007) B. Eckhardt, T. M. Schneider, B. Hof, and J. Westerweel, “Turbulence transition in pipe flow,” Ann. Rev. Fluid Mech. 39 , 447–468 (2007) . · doi ↗
- 4Eckhardt (2008) B. Eckhardt, “Turbulence transition in pipe flow: some open questions,” Nonlinearity 21 , T 1–T 11 (2008) . · doi ↗
- 5Faisst and Eckhardt (2003) H. Faisst and B. Eckhardt, “Traveling waves in pipe flow,” Phys. Rev. Lett. 91 , 224502 (2003) . · doi ↗
- 6Wedin and Kerswell (2004) H. Wedin and R. R. Kerswell, “Exact coherent structures in pipe flow: travelling wave solutions,” J. Fluid Mech. 508 , 333–371 (2004) . · doi ↗
- 7Duguet et al. (2008) Y. Duguet, C. C. T. Pringle, and R. R. Kerswell, “Relative periodic orbits in transitional pipe flow,” Phys. Fluids 20 , 114102 (2008) . · doi ↗
- 8Willis et al. (2013) A. P. Willis, P. Cvitanović, and M. Avila, “Revealing the state space of turbulent pipe flow by symmetry reduction,” J. Fluid Mech. 721 , 514–540 (2013) . · doi ↗
