Viscous-elastic dynamics of power-law fluids within an elastic cylinder

TL;DR

This paper investigates the complex interaction between non-Newtonian power-law fluids and elastic cylinders in microfluidic channels, deriving exact solutions and revealing how viscosity influences fluid front propagation, with implications for measuring fluid properties.

Contribution

It introduces a novel analytical framework for viscous-elastic dynamics of power-law fluids in elastic cylinders, including exact solutions and insights into front propagation behavior.

Findings

Exact solutions for pressure and deformation fields for shear-thinning and shear-thickening fluids.

Inversion of viscosity role compared to Stokes' problem, with shear-thinning fluids exhibiting front propagation.

The front propagation rate scales as t^{n/(n+1)}, enabling indirect measurement of the fluid's power-law index.

Abstract

In a wide range of applications, microfluidic channels are implemented in soft substrates. In such configurations, where fluidic inertia and compressibility are negligible, the propagation of fluids in channels is governed by a balance between fluid viscosity and elasticity of the surrounding solid. The viscous-elastic interactions between elastic substrates and non-Newtonian fluids are particularly of interest due to the dependence of viscosity on the state of the system. In this work, we study the fluid-structure interaction dynamics between an incompressible non-Newtonian fluid and a slender linearly elastic cylinder under the creeping flow regime. Considering power-law fluids and applying the thin shell approximation for the elastic cylinder, we obtain a non-homogeneous p-Laplacian equation governing the viscous-elastic dynamics. We present exact solutions for the pressure and deformation fields for various initial and boundary conditions for both shear-thinning and shear-thickening fluids. We show that in contrast to Stokes' problem where a compactly supported front is obtained for shear-thickening fluids, here the role of viscosity is inversed and such fronts are obtained for shear-thinning fluids. Furthermore, we demonstrate that for the case of a step in inlet pressure, the propagation rate of the front has a $t^{\frac{n}{n+1}}$ dependence on time ($t$), suggesting the ability to indirectly measure the power-law index ($n$) of shear-thinning liquids through measurements of elastic deformation.

Figures

Tables

Table 1. Table 1: Characteristic values of physical parameters corresponding to the flow of Xanthan gum suspension in a 500 μ ​ m 𝜇 𝑚 \mu m radius PDMS channel, subjected to external pressure of 1 × 10 4 1 superscript 10 4 1\times 10^{4} Pa. Resulting values for the non-dimensional parameters (in bold) indicate that the model holds well in this regime.

Property	Symbol	Value	Units	Based on
Wall thickness	${\tilde{h}}_m$	50	$\mu m$	input
Inner radius	${\tilde{r}}_c$	500	$\mu m$	input
Young’s modulus (PDMS)	$\tilde{E}$	$1\times {10}^6$	Pa	input
Poisson’s ratio	$\nu$	0.3	–	input
Power-law index (Xanthan suspension)	$n$	$0.55$	–	input
Effective viscosity (Xanthan suspension)	${\tilde{\mu }}_{eff}$	0.57	Pa sn	input
Characteristic time scale	${\tilde{t}}^{\ast }$	0.5	s	input
Characteristic pressure	${\tilde{p}}^{\ast }$	$1\times {10}^4$	Pa	input
Axial length scale	${\tilde{l}}^{\ast }$	$14$	cm	(36a)
Typical axial velocity	${\tilde{u}}_z^{\ast }$	0.91	$m{s}^{-1}$	(36b)
Typical radial deformation	${\tilde{d}}_r^{\ast }$	50	$\mu m$	(36c)
Typical axial deformation	${\tilde{d}}_z^{\ast }$	14	mm	(36d)
Slenderness	${ϵ}_1$	$3.5\times {10}^{-3}$	–	(37a)
Wall thinness	${ϵ}_2$	0.1	–	(37b)
Smallness of deformations	${ϵ}_3$	0.1	–	(37c)
Womersley number	$Wo$	$2.5\times {10}^{-2}$	–	(37d)

Equations

$${ϵ}_1=\frac{{\tilde{r}}_c}{{\tilde{l}}^{\ast }}\ll 1,{ϵ}_2=\frac{{\tilde{h}}_m}{{\tilde{r}}_c}\ll 1,Wo\ll 1,{ϵ}_3=\frac{{\tilde{d}}_r^{\ast }}{{\tilde{r}}_c}\sim \frac{{\tilde{d}}_z^{\ast }}{{\tilde{l}}^{\ast }}\ll 1,$$

$$\tilde{\eta }(\dot{\gamma })={\tilde{\mu }}_{eff}{\dot{\gamma }}^{n-1},$$

$$\frac{1}{\tilde{r}}\frac{\partial }{\partial \tilde{r}}\left(\tilde{r}{\tilde{u}}_r\right)+\frac{\partial {\tilde{u}}_z}{\partial \tilde{z}}=0,$$

$$\tilde{\rho }\left(\frac{\partial {\tilde{u}}_z}{\partial \tilde{t}}+{\tilde{u}}_r\frac{\partial {\tilde{u}}_z}{\partial \tilde{r}}+{\tilde{u}}_z\frac{\partial {\tilde{u}}_z}{\partial \tilde{z}}\right)=-\frac{\partial \tilde{p}}{\partial \tilde{z}}+\frac{\partial }{\partial \tilde{z}}\left({\tilde{\tau }}_{zz}\right)+\frac{1}{\tilde{r}}\frac{\partial }{\partial \tilde{r}}\left(\tilde{r}{\tilde{\tau }}_{rz}\right),$$

$$\tilde{\rho }\left(\frac{\partial {\tilde{u}}_r}{\partial t}+{\tilde{u}}_r\frac{\partial {\tilde{u}}_r}{\partial \tilde{r}}+{\tilde{u}}_z\frac{\partial {\tilde{u}}_r}{\partial \tilde{z}}\right)=-\frac{\partial \tilde{p}}{\partial \tilde{r}}+\frac{1}{\tilde{r}}\frac{\partial }{\partial \tilde{r}}\left(\tilde{r}{\tilde{\tau }}_{rr}\right)+\frac{\partial }{\partial \tilde{z}}\left({\tilde{\tau }}_{zr}\right)-\frac{{\tilde{\tau }}_{\theta \theta }}{\tilde{r}},$$

$$\tilde{\mathbf{\tau }}=2{\tilde{\mu }}_{eff}{\dot{\gamma }}^{n-1}\tilde{\mathbf{D}},$$

$${\tilde{u}}_z{|}_{\tilde{r}={\tilde{r}}_c+{\tilde{d}}_r}=\frac{\partial \widetilde{d_z}}{\partial \tilde{t}},{\tilde{u}}_r{|}_{\tilde{r}={\tilde{r}}_c+{\tilde{d}}_r}=\frac{\partial \widetilde{d_r}}{\partial \tilde{t}}+{\tilde{u}}_z{|}_{\tilde{r}={\tilde{r}}_c+{\tilde{d}}_r}\frac{\partial d_r}{\partial z},$$

$$\frac{{\tilde{u}}_r^{\ast }}{{\tilde{u}}_z^{\ast }}\sim \frac{{\tilde{r}}_c}{{\tilde{l}}^{\ast }}={ϵ}_1,$$

$$\dot{\gamma }\sim \frac{{\tilde{u}}_z^{\ast }}{{\tilde{r}}_c}\left(\left|\frac{\partial u_z}{\partial r}\right|+O({ϵ}_1^2)\right).$$

$${\tilde{t}}^{\ast }\sim \frac{{\tilde{d}}_z^{\ast }}{{\tilde{u}}_z^{\ast }}\sim \frac{{\tilde{d}}_r^{\ast }}{{\tilde{u}}_r^{\ast }}.$$

$${\tilde{u}}_z^{\ast }=\frac{{ϵ}_1^{1+1/n}{\tilde{p}}^{\ast 1/n}{\tilde{l}}^{\ast }}{{\tilde{\mu }}_{eff}^{1/n}},$$

$$Wo=\frac{\tilde{\rho }{\tilde{r}}_c^2}{\tilde{\eta }(\dot{\gamma }){\tilde{t}}^{\ast }}=\frac{{ϵ}_1}{{ϵ}_3}\frac{\tilde{\rho }{\tilde{u}}_z^{\ast 2-n}{\tilde{r}}_c^n}{{\tilde{\mu }}_{eff}}=\frac{\tilde{\rho }{\tilde{u}}_z^{\ast 2-n}{\tilde{r}}_c^{n+1}}{{\tilde{\mu }}_{eff}{d}_z^{\ast }}\ll 1,$$

$$\frac{1}{r}\frac{\partial }{\partial r}\left(r{u}_r\right)+\frac{\partial u_z}{\partial z}=0,$$

$$\frac{\partial p}{\partial z}=\frac{1}{r}\frac{\partial }{\partial r}\left(r{\left|\frac{\partial u_z}{\partial r}\right|}^{n-1}\frac{\partial u_z}{\partial r}\right)+O(Wo,{ϵ}_1^2),$$

$$\frac{\partial p}{\partial r}=O({ϵ}_1^{2}Wo,{ϵ}_1^2),$$

$$u_z{|}_{r=1+{ϵ}_3{d}_r}=\frac{\partial d{}_z}{\partial t},u_r{|}_{r=1+{ϵ}_3{d}_r}=\frac{\partial d_r}{\partial t}+{ϵ}_3\frac{\partial d{}_z}{\partial t}\frac{\partial d_r}{\partial z}.$$

$$\frac{\partial u_z}{\partial r}={\left|\frac{\partial p}{\partial z}\right|}^{\frac{1}{n}-1}\frac{\partial p}{\partial z}{\left(\frac{r}{2}\right)}^{\frac{1}{n}}.$$

$$u_z=\frac{1}{2^{\frac{1}{n}}}\frac{n}{n+1}{\left|\frac{\partial p}{\partial z}\right|}^{\frac{1}{n}-1}\frac{\partial p}{\partial z}\left(r^{1+\frac{1}{n}}-(1+{ϵ}_3{d}_r(z,t){)}^{1+\frac{1}{n}}\right)+\frac{\partial d_z}{\partial t}.$$

$$\frac{q}{2\pi }={\int }_{r=0}^{r=1+{ϵ}_3{d}_r}r{u}_{z}dr=-\frac{1}{2^{\frac{1+n}{n}}}\frac{n}{3n+1}{\left|\frac{\partial p}{\partial z}\right|}^{\frac{1}{n}-1}\frac{\partial p}{\partial z}+\frac{1}{2}\frac{\partial d_z}{\partial t}+O({ϵ}_3).$$

$$\frac{\partial d_r}{\partial t}+\frac{\partial }{\partial z}\left(-\frac{1}{2^{\frac{1+n}{n}}}\frac{n}{3n+1}{\left|\frac{\partial p}{\partial z}\right|}^{\frac{1}{n}-1}\frac{\partial p}{\partial z}+\frac{1}{2}\frac{\partial d_z}{\partial t}\right)=0.$$

$$\frac{\partial }{\partial \tilde{r}}\left(\tilde{r}{\tilde{\sigma }}_{rr}\right)+\frac{\partial }{\partial \tilde{z}}\left(\tilde{r}{\tilde{\sigma }}_{rz}\right)-{\tilde{\sigma }}_{\theta \theta }=0,\frac{\partial }{\partial \tilde{r}}\left(\tilde{r}{\tilde{\sigma }}_{rz}\right)+\frac{\partial }{\partial \tilde{z}}\left(\tilde{r}{\tilde{\sigma }}_{zz}\right)=0,$$

$${\tilde{e}}_{rr}=\frac{\partial {\tilde{d}}_r}{\partial \tilde{r}},{\tilde{e}}_{\theta \theta }=\frac{{\tilde{d}}_r}{\tilde{r}},{\tilde{e}}_{zz}=\frac{\partial {\tilde{d}}_z}{\partial \tilde{z}},{\tilde{e}}_{rz}=\frac{1}{2}\left(\frac{\partial {\tilde{d}}_r}{\partial \tilde{z}}+\frac{\partial {\tilde{d}}_z}{\partial \tilde{r}}\right),$$

$${\tilde{e}}_{rr}=\frac{1}{\tilde{E}}\left({\tilde{\sigma }}_{rr}-\nu \left({\tilde{\sigma }}_{zz}+{\tilde{\sigma }}_{\theta \theta }\right)\right),{\tilde{e}}_{\theta \theta }=\frac{1}{\tilde{E}}\left({\tilde{\sigma }}_{\theta \theta }-\nu \left({\tilde{\sigma }}_{zz}+{\tilde{\sigma }}_{rr}\right)\right),$$

$${\tilde{e}}_{zz}=\frac{1}{\tilde{E}}\left({\tilde{\sigma }}_{zz}-\nu \left({\tilde{\sigma }}_{\theta \theta }+{\tilde{\sigma }}_{rr}\right)\right),{\tilde{e}}_{rz}=\frac{1+\nu }{\tilde{E}}{\tilde{\sigma }}_{rz}.$$

$${\tilde{\sigma }}_{rr}{|}_{\tilde{r}={\tilde{r}}_c+{\tilde{d}}_r}=-\tilde{p}+2{\tilde{\mu }}_{eff}{\dot{\gamma }}^{n-1}\frac{\partial {\tilde{u}}_r}{\partial \tilde{r}},{\tilde{\sigma }}_{rz}{|}_{\tilde{r}={\tilde{r}}_c+{\tilde{d}}_r}={\tilde{\mu }}_{eff}{\dot{\gamma }}^{n-1}\left(\frac{\partial {\tilde{u}}_z}{\partial \tilde{r}}+\frac{\partial {\tilde{u}}_r}{\partial \tilde{z}}\right),$$

$${\tilde{\sigma }}_{rr}{|}_{\tilde{r}={\tilde{r}}_c+{\tilde{d}}_r+{\tilde{h}}_m}=-{\tilde{p}}_e,{\tilde{\sigma }}_{rz}{|}_{\tilde{r}={\tilde{r}}_c+{\tilde{d}}_r+{\tilde{h}}_m}=0.$$

$$d_r(z,t)=\left(1-\frac{\nu }{2}\right)p(z,t)-p_e,d_z(z,t)={\int }_0^z\left[(\frac{1}{2}-\nu )p(z,t)+\nu p_e\right]dz,$$

$$\frac{{\tilde{d}}_r^{\ast }}{{\tilde{r}}_c}=\frac{{\tilde{p}}^{\ast }}{{ϵ}_2\tilde{E}}\ll 1.$$

$$\frac{\partial p}{\partial t}-\frac{\partial }{\partial z}\left({\left|\frac{\partial p}{\partial z}\right|}^{\frac{1}{n}-1}\frac{\partial p}{\partial z}\right)=\frac{4-2\nu }{5-4\nu }\frac{\partial p_e}{\partial t},$$

$$p(z,t=0)=p_{initial}(z),p(z=0,t)=p_{inlet}(t),$$

$$t\to \frac{t}{\alpha },\alpha \triangleq \frac{1}{5-4\nu }\frac{1}{2^{\frac{1-n}{n}}}\frac{n}{3n+1}.$$

Full text

I Introduction

Viscous-elastic interaction (VEI) is a sub-field of fluid-structure interaction focused on fluid motion at low-Reynolds numbers DupratStone . VEI received considerable attention over recent years due to its relevance to a wide spectrum of natural processes and applications. These include biological flows Ku ; Nichols , geophysical and geological flows Hewitt2015 , suppression of viscous fingering instabilities Pihler-Puzovic2012 ; Pihler-Puzovic2013 ; AlHousseiny ; Pihler-Puzovic2014 ; Pihler-Puzovic2015 ; Peng2015 , manufacturing of micro-electro-mechanical systems Hosoi2004 ; Lister2013 , development of flexible microfluidic valves Holmes , and soft robotics Shepherd ; Martinez ; Marchese ; Matia .

In particular, the case of Newtonian viscous flow through elastic channels has been extensively studied as a model problem for internal VEI flows such as those found in veins and arteries Canic , collapsible tubes (see Grotberg ; Heil and DupratStone , ch. 8, and references therein), and soft actuators Elbaz2014 . In the field of micro and nano-fluidics, channels are often fabricated from flexible materials (e.g. PDMS), which deform under internal or external pressures Hardy ; deRutte . While many microfluidic applications involve non-Newtonian fluids, very few studies examined the deformation field and flow of non-Newtonian fluids within elastic configurations. In addition to an academic interest in such systems, we believe they hold potential for novel applications, particularly in soft robotics and soft lab-on-a-chip systems, where the choice of fluid properties could allow new modes of actuation and flow control.

In this work, we extend the recent study of Elbaz and Gat Elbaz2014 which focused on Newtonian fluids, and analyze the viscous-elastic dynamics of non-Newtonian power-law fluids flowing through a slender linearly elastic cylinder. In Sec. II, we present the problem formulation as well as the key assumptions we use in the derivation of the model. In Sec. III we apply the lubrication approximation to the flow equations and the thin shell approximation to the elasticity equations and derive a non-linear p-Laplacian governing equation for the viscous-elastic dynamics. We provide appropriate scaling for physical variables and clarify the range of validity of the model. In Sec. IV we examine the viscous-elastic dynamics of the system, subjected to various initial and boundary conditions, and obtain both exact and numerical solutions for the pressure and deformation fields for both shear-thinning and shear-thickening fluids. In Sec. V we demonstrate and discuss the inverse role of viscosity in the p-Laplacian equation governing the viscous-elastic problem, as compared to Stokes’ problems Pascal1992 ; Ai ; Pritchard . We conclude with a discussion of the results in Sec. VI.

II Problem formulation

We study the low-Reynolds number fluid-structure interaction dynamics of non-Newtonian, incompressible, axisymmetric flow through a slender linearly elastic cylinder. Figure 1 presents a schematic illustration of the configuration and coordinate system. We hereafter denote dimensional variables by tildes, normalized variables without tildes and characteristic values by an asterisk superscript.

We employ a cylindrical coordinate system $(\tilde{r},\tilde{\theta},\tilde{z})$ whose $\tilde{z}$ axes coincides with the symmetry line of the cylinder. The elastic cylinder is clamped at $\tilde{z}=0$, and has an initial uniform radius $\tilde{r}_{c}$, constant wall thickness $\tilde{h}_{m}$, Young’s modulus $\tilde{E}$ and Poisson’s ratio $\nu$. The fluid’s density is $\tilde{\rho}$ and is assumed to be constant, its velocity is $\mathbf{\tilde{\mathrm{\mathit{\bm{u}}}}}=\left(\tilde{u}_{r},\tilde{u}_{z}\right)$ and its pressure is $\tilde{p}$. The externally applied pressures $\tilde{p}_{inlet}$ and $\tilde{p}_{e}$ at the inlet and the outer shell boundary, result in viscous-elastic interaction between the non-Newtonian fluid and the elastic cylinder, which in turn leads to both radial, $\tilde{d}_{r}$, and axial, $\tilde{d}_{z}$, deformation fields. The characteristic length scale in the $\bm{\hat{z}}$ direction is denoted by $\tilde{l}^{*}$, while $\tilde{d}_{r}^{*}$ and $\tilde{d}_{z}^{*}$ are the characteristic radial and axial deformations, respectively.

We further assume a slender geometry, thin walls, negligible fluidic inertia effects (small Womersley number), and small deformations in both the radial and axial directions. These can be summarized as

[TABLE]

allowing the use of lubrication and thin elastic shell approximations. The exact expression for the Womersley number appearing in (1) depends on the viscosity model, which we will define in (12). Furthermore, defining both ratios, $\tilde{d}_{r}^{*}/\tilde{r}_{c}$ and $\tilde{d}_{z}^{*}/\tilde{l}^{*}$, to be of scale $\epsilon_{3}$ emanates from mass conservation of the fluid and the constitutive law of the elastic structure, as shown in the following section.

III Viscous-elastic governing equations for power-law fluids

III.1 Fluidic problem

Description of the non-Newtonian fluid rheology requires the use of a specific constitutive model. Of all models, the commonly used power-law model approximately describes the behavior of shear-dependent non-Newtonian fluids over an intermediate range of shear rates, and at the same time is sufficiently simple to allow analytical treatment. Therefore, throughout this work we consider the following constitutive model for the viscosity $\tilde{\eta}$ Bird

[TABLE]

where $n$ is the dimensionless power-law index, $\tilde{\mu}_{eff}$ is the effective constant viscosity with units $\mathrm{Pa\,s^{n}}$ and $\dot{\gamma}$ is the shear rate defined as $\dot{\gamma}=\sqrt{2\bm{\tilde{D}}:\bm{\tilde{D}}},$ where $\bm{\tilde{D}}$ is the rate of deformation tensor given by $\bm{\tilde{D}}=(\bm{\tilde{\nabla}\tilde{u}}+\bm{\tilde{\nabla}\tilde{u}^{\dagger}})/2.$ The power-law index lies in the range $0<n<1$ for shear-thinning fluids, and $n>1$ for shear-thickening fluids, while $n=1$ represents the case of a Newtonian fluid.

Based on the assumptions mentioned above, the relevant governing equations and boundary conditions are the continuity equation

[TABLE]

the momentum equation in the $\bm{\hat{z}}$ direction

[TABLE]

the momentum equation in the $\bm{\hat{r}}$ direction

[TABLE]

the constitutive equation for the stress tensor $\bm{\tilde{\tau}}$ of power-law fluid

[TABLE]

as well as the no-slip and the no-penetration conditions at the solid-liquid interface

[TABLE]

where $\tilde{t}$ is the time.

Scaling by the characteristic dimensions, we define the normalized coordinates, $(r,z)=(\tilde{r}/\tilde{r}_{c},\tilde{z}/\tilde{l}^{*})$, normalized velocity, $(u_{r},u_{z})=(\tilde{u}_{r}/\tilde{u}_{r}^{*},\tilde{u}_{z}/\tilde{u}_{z}^{*})$, normalized pressure, $p=\tilde{p}/\tilde{p}^{*}$, and normalized time, $t=\tilde{t}/\tilde{t}^{*}$, where $\tilde{t}^{*}$ is a characteristic time scale, $\tilde{p}^{*}$ is a characteristic pressure, $\tilde{u}_{r}^{*}$ and $\tilde{u}_{r}^{*}$ are radial and axial characteristic velocities, respectively.

Substituting the normalized parameters into the continuity equation (3), order of magnitude analysis yields the balance between radial and axial velocities

[TABLE]

and suggests the scaling for shear rate, $\dot{\gamma}$ ,

[TABLE]

Order of magnitude analysis of the boundary conditions (7) yields the characteristic viscous-elastic times scale

[TABLE]

Substituting (8) into (4) and performing order of magnitude analysis, we obtain the axial characteristic velocity in terms of the characteristic pressure

[TABLE]

as well as the condition for a negligible Womersley number in terms of relevant physical quantities

[TABLE]

representing the ratio between the inertial-viscous time scale $\tilde{\rho}\tilde{r}_{c}^{2}/\tilde{\eta}(\dot{\gamma})=\tilde{\rho}\tilde{r}_{c}^{2}/\tilde{\mu}_{eff}(\tilde{u}_{z}^{*}/\tilde{r}_{c})^{n-1}$ and the viscous-elastic time scale $\tilde{t}^{*}=\tilde{d}_{z}^{*}/\tilde{u}_{z}^{*}.$

Substituting (8)-(12) into (3)-(7) and applying the lubrication approximation results in the following set of normalized equations and boundary conditions,

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Integrating (14) with respect to $r$ and applying the symmetry condition at $r=0$ yields

[TABLE]

Integrating (17) again with respect to $r$ and applying the no-slip boundary conditions (16), results in

[TABLE]

Defining volume flux as $q=2\pi\int_{r=0}^{r=1+\epsilon_{3}d_{r}}ru_{z}dr$, from (18) it follows

[TABLE]

Utilizing (16) and (19), and then integrating (13) with respect to $r$, we relate the pressure and deformation fields

[TABLE]

III.2 Elastic problem

A similar elastic configuration was considered by Elbaz and Gat Elbaz2014 , who provided relations between the pressure and deformation fields. For completeness, we here present the governing equations and these key relations, adapted to our notation.

Focusing on viscous-elastic time scales which are typically significantly longer than elastic inertial time scales, and further neglecting body forces, the deformation field of an axisymmetric linearly elastic material is described by the momentum equations in the $\bm{\hat{r}}$ and $\bm{\hat{z}}$ directions (see Timoshenko, ; Howell, ),

[TABLE]

the strain-displacement relations for small deformations,

[TABLE]

and Hooke’s law,

[TABLE]

The governing equations (21) are subjected to the following boundary conditions applied by the liquid

[TABLE]

and the external stress

[TABLE]

We follow the derivation of Elbaz and Gat Elbaz2014 and obtain the normalized governing equations for the radial and axial deformations in terms of the pressure

[TABLE]

as well as the characteristic radial deformation, $\tilde{d}_{r}^{*}$, in terms of characteristic pressure

[TABLE]

In the Supplemental Material SI , we show that early scaling of the radial differential, $d\tilde{r}$, in (21) by $\tilde{h}_{m}$ (instead of $\tilde{r}_{c}$ (see Elbaz2014, )) is a more natural choice which leads to significant simplification in the derivation of (26) and (27) in the current configuration.

III.3 P-Laplacian equation governing the pressure field and its asymptotic approximation

To combine the fluidic and elastic problems, we substitute elastic relation (26) into the fluidic relation (20), resulting in a single equation in terms of the pressure

[TABLE]

subjected to the following initial and boundary conditions

[TABLE]

where, for the convenience, we have rescaled the time, $t$, according to

[TABLE]

Equation (28) is a non-homogeneous one dimensional p-Laplacian equation describing the pressure in an elastic cylinder containing a power-law fluid subjected to external forcing.

Similar p-Laplacian equations are encountered in other problems involving the use of non-Newtonian power-law fluids. For example, the p-Laplacian equation represents balance between compressibility and viscosity for the flow of a weakly compressible non-Newtonian fluid in a porous medium Pascal1985 ; Pascal1991 ; DiFederico , as well as a balance between viscosity and inertia in Stokes’ first Pascal1992 ; Duffy and second problems Ai ; Pritchard , where the operator is applied to the fluid velocity. We here rely on the previous mathematical analysis of the p-Laplacian equation to provide insight into our viscous-elastic problem, such as self-similarity and compact support of the solutions for pressure and deformation fields.

For the case of a non-Newtonian power-law fluid that exhibits a weak shear-thinning or shear-thickening behavior, we follow Ross, Wilson and Duffy Ross and define the auxiliary small parameter $\varepsilon$,

[TABLE]

and the asymptotic expansion

[TABLE]

where $\left|\varepsilon\right|\ll 1$ is a small parameter which is positive for shear-thinning and negative for shear-thickening behaviors. Using (31) and (32) the expansion $1/n\sim 1+\varepsilon+O(\varepsilon^{2})$, yields

[TABLE]

We note that the expansion (33) is valid only provided the leading-order pressure gradient is non-zero and bounded throughout the domain. Utilizing (31)-(33), we obtain the leading-order and first-order corrections of the governing equation (28)

[TABLE]

and

[TABLE]

Both (34) and (35) are parabolic equations where the inhomogeneous part of (34) is related to the external forcing (see Elbaz2014, ), and the inhomogeneous part of (35) emanates from the non-Newtonian response of the fluid to the leading-order pressure gradients.

III.4 Characteristic physical parameters and model range of validity

We here summarize the characteristic values of relevant physical parameters and the range of validity of the current model. From order of magnitude analysis, the characteristic length scale $\tilde{l}^{*}$, the characteristic axial velocity $\tilde{u}_{z}^{*}$, and the characteristic deformations $\tilde{d}_{r}^{*}$ and $\tilde{d}_{z}^{*}$ are given by

[TABLE]

Since in this work we focus on semi-infinite configuration, there is no inherent axial length scale to the problem and thus scaling analysis provides only a relation between $\tilde{l}^{*}$ and $\tilde{t}^{*}$. We define the characteristic pressure $\tilde{p}^{*}$ and the characteristic time scale $\tilde{t}^{*}$ from the boundary condition $\tilde{p}(0,\tilde{t})=\tilde{p}_{inlet}(\tilde{t})$ with characteristic amplitude $\tilde{p}_{0}$ and time scale $\tilde{t}_{0}$, by setting $\tilde{p}^{*}=\tilde{p}_{0}$ and $\tilde{t}^{*}=\tilde{t}_{0}$.

Substituting (36a)-(36d) into (1) provides formulation of the assumptions made in the model in terms of known physical and geometrical parameters:

[TABLE]

To illustrate the applicability of these assumptions in soft-robotic and lab-on a-chip applications, we examine a typical micro-configuration made of PDMS (polydimethylsiloxane, $\tilde{E}=10^{6}\,\mathrm{Pa}$, Cheng, ) and containing a suspension of Xanthan gum behaving as power-law viscous fluid (n=0.55,\,\tilde{\mu}_{eff}=0.57\,$$\mathrm{Pa\,s^{n}}, Longo, ). Table 1 presents the characteristic values of physical and geometric parameters, readily satisfying the assumptions used in the current analysis.

IV Results

In this section, we examine the viscous-elastic dynamics resulting from time-varying inlet pressure. We first obtain the pressure field by solving the governing equation (28) for arbitrary values of $n$. Then, using (19) and (26), we provide the corresponding volume flux and deformation field.

IV.1 Instantaneous injection of mass at the inlet

We examine the viscous-elastic dynamics due to an initial localized pressure at $z=0$, which may be represented by

[TABLE]

where $\delta$ is Dirac’s delta function. We note that (38) physically results from instantaneous release of a fluid mass, $\tilde{m}$, into the elastic cylinder at the inlet at $t=0$ after which ($t>0$), the inlet is sealed and thus

[TABLE]

Mass is thus conserved globally according to

[TABLE]

which using (26) can be rewritten as

[TABLE]

In the following, we assume that the fluid mass, $\tilde{m}$, is such that the right-hand side of (41) is unity, $\int_{0}^{\infty}pdz=1$. Note that the conservation law, (41), can be also derived by integrating (28) with respect to $z$ from [math] to $\infty$, using (38). Following Di Federico and Ciriello DiFederico , who focused on shear-thinning fluids and considered an instantaneous mass injection of a weakly compressible power-law fluid into porous medium, we construct self-similar solutions by defining a similarity variable $\xi$

[TABLE]

and self-similar profile $f$

[TABLE]

For a Newtonian fluid ($n=1$) the solution is well-known

[TABLE]

For non-Newtonian fluids, substituting (43) into (28) and (41) yields an ordinary differential equation for the self-similar profile and the corresponding conservation law,

[TABLE]

and

[TABLE]

Integrating (45) with respect to $\xi$ and applying the far-field condition of vanishing pressure and its derivatives $f=df/d\xi\rightarrow 0$ as $\xi\rightarrow\infty$, yield

[TABLE]

Since, the pressure should be positive and the pressure gradient should be negative, i.e. $f\geq 0$ and $df/d\xi\leq 0$, we can drop the absolute value in (47) and obtain

[TABLE]

Integrating (48) with respect to $\xi$ leads to

[TABLE]

where the minus and plus signs stand for $n<1$ and $n>1$, respectively, $k$ is defined as

[TABLE]

$x_{+}=\max\left\{x,0\right\}$ and $\xi_{*}$ is a positive constant uniquely determined from (46) and given by

[TABLE]

where $\Gamma$ is the Gamma function Abramowitz .

For a shear-thinning fluid, $n<1$, the self-similar profile $f(\xi)$ becomes uniformly zero for $\xi>\xi_{front}$ thus indicating the existence of a front which moves according to $z_{front}(t)=\xi_{front}t^{\frac{n}{2}}$. The condition $f=df/d\xi\rightarrow 0$ holds now at $\xi=\xi_{front}$. In contrast to the shear-thinning fluid, the self-similar profile $f(\xi)$ of the shear-thickening fluid, $n>1$, approaches zero only as $\xi\rightarrow\infty$, similarly to the case of a Newtonian fluid.

The closed-form solutions, (49), provide the associated radial deformation

[TABLE]

which admits the same self-similar profile $f$. Substituting (44) and (49) into (26), we obtain the corresponding self-similar axial deformation in the form $d_{z}(z,t)=\left(0.5-\nu\right)\int_{0}^{z}p(z,t)dz=\left(0.5-\nu\right)D_{z}(\xi)$, where $D_{z}(\xi)$ is expressed in terms of generalized hyper-geometric functions

[TABLE]

where ${}_{2}F_{1}\left(a,b;c;z\right)$ is the Gaussian hyper-geometric function Abramowitz ; Slater . The plus and minus signs correspond to $n<1$ and $n>1$, respectively. Rescaling the time in the expression for the volume flux, (19), according to (30), we find

[TABLE]

Expressing (54) in terms of the similarity variable $\xi$ yields

[TABLE]

where $f(\xi)$ is given in (44) and in (49) for the Newtonian fluid and the non-Newtonian fluid, respectively, and $Q(\xi)=\xi f(\xi)$.

Figure (2) summarizes the results of viscous-elastic interaction due to a pressure impulse applied at the inlet for $n=0.5,1,1.5$. All solutions corresponding to shear-thinning fluids are compactly supported, in contrast to solutions of Newtonian and shear-thickening fluids. Figures 2(a) and 2(b) present the pressure distribution and the radial deformation versus the axial coordinate $z$, for different times. From (43) and (52), the decay rate in this case scales as $t^{-n/2}$, thus indicating that for long times the pressure and the radial deformation of shear-thickening fluids decay (and propagate) faster than in the case of Newtonian and shear-thinning fluids, consistent with the results of Figs. 2(a) and 2(b). Figures 2(c) and 2(d) show the axial deformation and volume flux versus axial coordinate. In this case, the decay rate is similar for Newtonian and non-Newtonian fluids (i.e. is independent of $n$).

It worth to note that at early-times, $t\ll 1$, the shear-thickening fluids have the lowest propagation rate and are thus more localized (but not compactly supported) than the Newtonian and shear-thinning fluids. This is inverted at longer times as propagation rate of the shear-thickening fluids increases, and they are least localized.

Furthermore, we note that the current analysis is limited to an intermediate range of shear rates in accordance the validity range of the power-law model. For sufficiently long times, as $\partial p/\partial z$ diminishes, the shear rate will also diminish, driving the effective viscosity outside the validity range of the power-law model. Thus, the current analysis will become invalid at long times. In practice, in such cases the effective viscosity would approach a constant value, and the dynamics could be described by a linear diffusion equation.

IV.2 Sudden change of inlet pressure

We here examine the viscous-elastic dynamics associated with a sudden change of inlet pressure

[TABLE]

from an initially uniform pressure within the cylinder, $p(z,t=0)=p_{i}$. Here $H$ stands for the Heaviside function, whereas $p_{i}$ and $p_{0}$ are the magnitudes of the initial and the applied pressure, respectively. The pressure gradient is positive for $p_{i}>p_{0}$ and negative for $p_{i}<p_{0}$, corresponding to filling or discharging of the elastic tube, respectively.

For convenience, we define $\beta=\mathrm{sgn}\left(p_{i}-p_{0}\right)$ and rewrite (28) as

[TABLE]

Following Pascal and Pascal Pascal1985 , we define a similarity variable $\xi$

[TABLE]

and require

[TABLE]

The corresponding solution for a Newtonian fluid solution is

[TABLE]

For the non-Newtonian case, we substitute (58) and (59) into (57), obtaining an ordinary differential equation for the pressure

[TABLE]

with boundary conditions for $t>0$

[TABLE]

Integrating once (61) with respect to $\xi$ yields

[TABLE]

where the upper and lower signs correspond to $n<1$ and $n>1$, respectively. The parameter $k$ is defined as

[TABLE]

and $\xi_{*}$ is a constant to be determined from boundary conditions on the pressure, in contrast to Sec. IV A where a similar constant is found from conservation of mass.

Integrating (63) with respect to $\xi$ and applying the boundary condition (62) at $\xi=0$, leads to

[TABLE]

To determine the constant $\xi_{*}$ for different values of $n$, we utilize the boundary condition $p(\xi\rightarrow\infty)=0$, yielding

[TABLE]

From (63) we obtain that for shear-thinning fluids $dp/d\xi=0$ at $\xi=\xi_{*}=\xi_{front}$, thus indicating the existence of a front moving according to $z_{front}(t)=\xi_{front}t^{n/(n+1)}$. On the other hand, for shear-thickening fluids it follows that $dp/d\xi=0$ only as $\xi\rightarrow\infty$, and therefore no front exists.

Figures 3(a) and 3(b) present the resulting pressure distribution, (65), and the radial deformation, $d_{r}(z,t)=\left(1-0.5\nu\right)p(z,t)=\left(1-0.5\nu\right)p(\xi)$, corresponding to the case of $p_{0}=1,\,p_{i}=0,\,\beta=-1$, for different times. Figure 3(c) shows the axial deformation, $d_{z}(z,t)$, given as $d_{z}(z,t)=\left(0.5-\nu\right)\int_{0}^{z}p(z,t)dz=\left(0.5-\nu\right)zD_{z}(\xi)/\xi$, where for the Newtonian fluid ($n=1$) $D_{z}(\xi)$ is the following self-similar profile,

[TABLE]

and for a non-Newtonian fluid is given by

[TABLE]

where the upper signs in (68) correspond to $n<1$ and the lower signs correspond to $n>1$, respectively. The solution for volume flux is presented in Fig. 3(d) and is obtained using (54),

[TABLE]

where the expressions of $p(\xi)$ and $D_{z}(\xi)$ are given in (60) and (67) for the Newtonian fluid, and in (65) and (68) for the non-Newtonian fluid.

While the Newtonian and non-Newtonian fluids have the same decay rate in pressure and radial deformation, the shear-thinning fluids exhibit the fastest decay in volume flux and the shear-thickening fluids exhibit the fastest growth in axial deformation.

IV.3 Steady-state dynamics due to an oscillatory inlet pressure

In this section, we examine the steady non-Newtonian viscous-elastic interaction resulting from an oscillatory pressure at the inlet

[TABLE]

where $\omega$ is the non-dimensional magnitude of the base frequency $1/\tilde{t}^{*}$, and the far-field boundary conditions is $p(z\rightarrow\infty,t)=0$.

The corresponding steady-state solution for a Newtonian fluid is

[TABLE]

For a non-Newtonian fluid, an exact solution of (28) together with the boundary condition (70) is difficult to obtain. In the Supplemental Material SI we present an asymptotic solution in terms of Green’s functions, but it is difficult to obtain direct insight from this solution as some of the resulting expressions can be evaluated only numerically. Nevertheless, some physical insight on this problem can be obtained by relying on the first-order asymptotic correction without the need for an exact solution. As opposed to Newtonian fluids where the resulting pressure oscillates only with the imposed frequency $\omega$ (see (71)), we expect non-Newtonian fluids to respond at various frequencies due to the non-linearity of (28). To explore these active frequencies $\omega_{ai}$, we examine the source term of the first-order equation (35) that reflects the non-Newtonian response of the fluid

[TABLE]

By decomposing the right-hand side of (72) into three contributions:

[TABLE]

and adopting for convenience the following notations

[TABLE]

where $m=1,2,3...$, it is evident that the first two contributions (73a) and (73b) act with the imposed frequency $\omega$.

Expanding the logarithmic term in (73c) in a Taylor series yields

[TABLE]

and thus we conclude that the third contribution (73c) adds additional forcing frequencies which are the multiplication of the imposed frequency with odd natural numbers ($3\omega,\,5\omega,\,...$).

To validate the Taylor-expansion approximation made in (75), we performed fast Fourier transform (FFT) analysis on a numerical calculation for the source term of (72). The results presented in Fig. 4(a) show full agreement with analytically predicted frequencies of the first-order source term. To further explore the viscous-elastic dynamics of the problem, we solved the governing equation (28) numerically on a fixed domain $(0,L)$ using a second-order explicit finite-difference scheme. The detailed description of the numerical scheme is presented in the Supplemental Material SI . In all simulations we used values of $L=10$ and $\omega=2\pi$.

In Figs. 4(b) and 4(c) we present the resulting frequencies of FFT analysis applied on the results of numerical simulations for shear-thinning and shear-thickening fluids, showing excellent agreement with the analytically predicted frequencies obtained from the first-order source term. From Figs. 4(b) and 4(c) we obtain that for actuation at the viscous-elastic time scale ($\omega\sim O(1)$), the dominant frequencies are

[TABLE]

and therefore the steady-state time variation of the pressure at any location $z$ can be approximately described by the superposition of only three sinusoidal waveforms

[TABLE]

where $A_{1}\left(z\right)$, $A_{2}\left(z\right)$ and $A_{3}\left(z\right)$ are complex functions solely depending on $z$ and ”c.c.” represents the complex conjugate of preceding terms.

Figure (5) summarizes the steady-state results of viscous-elastic interaction due to an oscillatory pressure (70). Figures 5(a,b,c) present the steady-state pressure distribution versus axial coordinate for one period of oscillations, obtained from numerical solution of (28) and corresponding to shear-thinning, Newtonian and shear-thickening fluids, respectively. Clearly, the shear-thinning fluid is characterized by the most rapid pressure decay rate similarly to the results in previous sections, thus suggesting the existence of a propagation front (which we however did not prove for this case).

In Fig. 5(d), we compare the numerical solutions (solid lines) and approximate solutions (dashed lines), (77), for the pressure variation as a function of time at $z=1.4$, showing good agreement. As can be inferred from Fig. 5(d), in the steady-state the shear-thinning fluid exhibits a non-sinusoidal ”saw-tooth” waveform, whereas the shear-thickening is characterized by a non-sinusoidal ”shark-tooth” waveform.

V The inverse role of viscosity between viscous-elastic and Stokes’ problems

In Secs. IV B and IV C, we studied analytically and numerically the viscous-elastic interaction of power-law fluids through an elastic cylinder resulting from a constant or an oscillating pressure at the inlet. Specifically, we solved a homogeneous p-Laplacian equation that can be written in dimensional form as

[TABLE]

We note that for $n=1$ (78) reduces to an unforced version of equation (2.53) derived by Elbaz and Gat Elbaz2014 .

This p-Laplacian equation is found in various fluid mechanics problems including Stokes’ first and second problems for power-law fluids, which involve the solution of the following p-Laplacian equation for $\tilde{u}_{z}$

[TABLE]

subjected to a constant or an oscillatory velocity at $\tilde{z}=0$. The Stokes’ first problem was considered by Pascal Pascal1992 who showed that the velocity profile is compactly supported for shear-thickening fluids. Pritchard et al. Pritchard presented semi-analytical, self-similar and periodic solutions for the flow of a power-law fluid driven by a non-sinusoidal oscillating waveform at the wall and showed that in this case these solutions predict a finite penetration length $\tilde{l}_{p}$ for shear-thickening fluids, beyond which the fluid is essentially unaffected by external forcing and is motionless. To confirm these solutions, Pritchard et al. Pritchard studied numerically the Stokes’ second problem (on a finite numerical domain) with the sinusoidal boundary condition and the numerical results corresponding to shear-thickening fluids decay rapidly with the distance from the wall thus suggesting the existence of a finite penetration length. In addition, Pritchard et al. Pritchard obtained that the velocity of shear-thickening fluids adopt a ”saw-tooth” oscillating form, while shear-thinning fluids have a ”shark-tooth” oscillating form.

Clearly, the results we obtained in Secs. IV B and IV C are opposite to the results of Pascal Pascal1992 and Pritchard et al. Pritchard . To explain this difference we should examine the effect of viscosity in the p-Laplacian equations, (78) and (79), on the fluid penetration. While in Stokes’ equation, (79), viscosity contributes to diffusivity and thus to fluid penetration, in the viscous-elastic equation, (78), viscosity reduces the diffusivity and fluid penetration. This can be observed more clearly from the expressions of the penetration length, $\tilde{l}_{p}$, which can be obtained from scaling arguments as

[TABLE]

[TABLE]

where subscript VE indicates a viscous-elastic problem parameter, S indicates a Stokes’ problem parameter. Noting that for shear-thickening fluids ($n>1$), both viscosities (S) and (VE) tend to zero as the gradient of the velocity (or pressure) goes to zero, it follows from (80) and (81) that while for the Stokes’ problem the penetration length, $\tilde{l}_{p,S}$, goes to zero, for the viscous-elastic problem the penetration length tends to infinity, $\tilde{l}_{p,VE}\rightarrow\infty$, as $\tilde{\eta}(\dot{\gamma})\rightarrow 0$.

On the other hand, for shear-thinning fluids, both viscosities (S) and (VE) diverge as the gradient of the velocity (or pressure) goes to zero, and thus Stokes’ problem yields an infinite penetration length, $\tilde{l}_{p,S}\rightarrow\infty$, whereas the viscous-elastic problem exhibits a finite penetration length indicating the existence of the front in this case.

These conclusions are fully consistent with our results from Secs. IV B and IV C and with the results of Pascal Pascal1992 and Pritchard et al. Pritchard . Shear-thinning and shear-thickening fluids play opposite roles in the Stokes’ and the viscous-elastic problems.

VI Discussion and concluding remarks

In this work, we studied the viscous-elastic interaction and the flow of power-law fluids through an elastic cylinder. Applying the lubrication and the elastic shell approximations, we derived a p-Laplacian governing equation relating the fluidic pressure and the external forces acting on the cylinder. We showed that in this equation viscosity acts as an inverse diffusion coefficient, such that the familiar roles of shear-thinning and shear-thickening are inversed relative to similar equations that appear in Stokes’ problems. For example, for instantaneous mass injection and a step pressure applied at the inlet, our analysis revealed that a compactly supported propagation front exists solely for shear-thinning fluids, in contrast to Stokes’ problem where a compactly supported front is obtained for shear-thickening fluids. Furthermore, the scaling for penetration length, (80), based on a power-law viscosity model, further suggests that any Dirichlet-type boundary condition imposed at $z=0$, would result in compactly supported propagation for shear-thinning fluids, consistent with our numerical observations in Sec. IV C.

Throughout this work we have mainly considered the viscous-elastic dynamics due to imposed pressure at the inlet. However, we note that the solutions obtained in Sec. IV can be used to construct new solutions resulting from an external pressure distribution $p_{e}(z,t)$. Applying the Laplace transform on (28) we obtain that the solution for the pressure in Sec. IV A for impulse pressure applied at the inlet of a semi-infinite cylinder also holds for the case of a sudden localized pressure, $p_{e}=p_{e}H(t)\delta(z)$, applied to an infinitely long cylinder. Similarly, the solution for the pressure in the case of a discharging elastic tube ((65) with $\beta>0$, $p_{0}=0$) also holds for discharge due to suddenly applied constant pressure $p_{e}=p_{e}H(t)$. While the solutions for the pressure in these case are identical, the solutions for the corresponding deformation fields and the volume fluxes are different, and can be obtained from (26) and (54).

For weakly non-Newtonian fluids, we obtained the leading- and first-order asymptotic approximations, and showed (see Sec. IV C) that the source term in first-order equation, (35), is important as it enables physical insight on the non-Newtonian behavior without a direct solution.

The regime of an incompressible fluid with negligible inertia is commonly encountered in lab-on-chip and soft robotics applications. For the former, controlling the propagation of mass flux within microfluidic channels could be a mechanism the timing of simple delivery on chip, or controlling the extent of mixing between reagents. For soft robotics, control with non-Newtonian liquids holds to potential for additional degrees of freedom in modulation, and particularly with the use of compactly supported deformations where part of the structure can be guaranteed to stay stationary while another is translating. Furthermore, we believe that the dependence of the front on time and on specific value of $n$, can be used to indirectly measure the rheological properties of power-law shear-thinning liquids through measurements of deformation and front location.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme, grant agreement no. 678734 (MetamorphChip). We gratefully acknowledge support by the Israel Science Foundation (grant no. 818/13). E.B. is supported by the Adams Fellowship Program of the Israel Academy of Sciences and Humanities.