2.1 The flow system

In our analysis, we consider a horizontal capillary channel of cross-sectional area $A$ and perimeter $s$ . It is initially prefilled with a liquid column of length $h_{0}$ and in contact with a drop of the same substance. The empty part of the conduit is totally unwetted as long as it is not being occupied by the intruding fluid. The cross-sectional dimension of the tube is assumed to be much smaller than the size of the drop so that it can be approximated as an infinite reservoir with a free surface. The contact point between the drop and the channel is very near to the free surface so that the flow inside the vessel is solely driven by capillary action without any gravitational influence. The system is schematically described in figure 1. It also shows side by side the arrangement considered in earlier works where pressure at the entry point has to be evaluated by modelling the outside flow. Such a calculation is not needed in our case as the atmospheric condition prevails at the entrance of our considered system.

Figure 1. (a) The system considered in our analysis where a capillary channel transports fluid from a drop with a free surface due to the surface tension effect. (b) Earlier works have considered a slightly different problem where the supplying fluid body is confined in an infinite chamber requiring a mathematical model for entry loss before the entrance.

The conduit is assumed to be prefilled by a liquid column of initial length $h_{0}$ before the start of the dynamics. We anticipate the need for such partial filling for accurate rheological measurements. This is because the discrepancy between the theoretical prediction and the experimental data can be unnecessarily high due to entry effects if the penetrated distance is of the order of the cross-sectional dimensions. In the rheometric set-up, $h_{0}$ can be introduced by either using a retractable barrier or dissolvable obstacle in the path of the fluid. As a result, the motion in the system vanishes for some time after reaching a distance, and resumes from a stationary condition after removal of the blockage. Accordingly, the additional parameter $h_{0}$ in the analysis will provide us with the mathematical flexibility to accommodate future experimental necessity. If this feature is not needed, $h_{0}$ can simply be taken as 0.

The encroaching medium is a complex fluid whose dynamics is governed by a linear viscoelastic equation. As a result, the instantaneous stress $\unicode[STIX]{x1D748}$ in the liquid domain is dependent on a weighted history of the velocity field $\boldsymbol{v}$ so that

Here, superscript T represents transpose, $\unicode[STIX]{x1D707}_{0}$ is the steady-state viscosity and $\hat{L}_{t}$ is an integral operator in time $t$ and $\unicode[STIX]{x1D70F}$

The dimensionless weight function $M$ manifests the viscoelastic effect by including the cumulative influence of past strain rates on the instantaneous stress at a point in the medium (Rousse Jr Reference Rousse1953). Different mildly viscoelastic liquids can be characterized by the memory function $M$ . For example, a Maxwell fluid is defined by an exponential $M$ . For purely viscous case, $M$ becomes a Dirac delta function, so that stress depends only on the strain rate at that instant and $\hat{L}_{t}$ transforms to an identity operator. Fourier transform of $M$ yields the complex viscosity normalized by $\unicode[STIX]{x1D707}_{0}$ . Consequently, it would be also possible to evaluate the frequency-dependent viscoelastic coefficients $G^{\prime }$ and $G^{\prime \prime }$ of the medium by relating them to the real and imaginary parts of the Fourier transform of $M$ .

The definition of steady-state viscosity $\unicode[STIX]{x1D707}_{0}$ implies the zeroth moment of the memory function $M$ is unity. The convolution in (2.2) explicitly involves the effective time scale $\unicode[STIX]{x1D70F}_{v}$ for viscoelastic relaxation. It can be uniquely obtained by ensuring the first moment of $M$ to be unity. In that case, even if the medium involves multiple scales for relaxation, their weighted average would be expressed by $\unicode[STIX]{x1D70F}_{v}$ . In other words, the definitions of $\unicode[STIX]{x1D707}_{0}$ and $\unicode[STIX]{x1D70F}_{v}$ impose two constraints on the functional distribution of the non-dimensional weight function $M$ by specifying two of its moments

Properties of any mildly viscoelastic liquid can be represented by $\unicode[STIX]{x1D707}_{0}$ , $\unicode[STIX]{x1D70F}_{v}$ and $M(\unicode[STIX]{x1D701})$ with $\unicode[STIX]{x1D701}$ being the normalized time interval. Even if multiple time scales of relaxation are present, these can be described by combining $\unicode[STIX]{x1D70F}_{v}$ and $M(\unicode[STIX]{x1D701})$ , where the latter would reveal non-trivial non-monotonic variations.

The fluid domain mainly contains unidirectional but transient flow as long as $h_{0}\gg s$ . So apart from the two ends of the liquid column, the dynamics is represented by the velocity component $v_{z}$ along the direction $z$ of the channel length. Accordingly, the system is governed by the following momentum equation involving $v_{z}$ and pressure $p$

where $\unicode[STIX]{x1D70C}$ is the density of the encroaching medium and $\unicode[STIX]{x1D735}_{\Vert }$ is the gradient in cross-sectional plane. When (2.4) is combined with the integral mass balance equation involving average over cross-sectional area $A$

the penetrated length $h$ can be evaluated as function of $t$ .

2.2 Scaling and non-dimensionalization

We solve (2.4) and (2.5) in dimensionless form to provide a concise and general description of the transport phenomenon. This requires identification of proper scales for the variables relevant in the dynamics of the system.

The scales for the cross-sectional coordinates and pressure field are readily available from the geometric and parametric properties. The former is defined as $b_{c}$ , and identified as the area-to-perimeter ratio. The latter, denoted $P_{s}$ , is determined by equating the average forces due to pressure and capillary action. This means that $P_{s}$ is the cross-sectional mean of the Laplace pressure near the propagating interface. Thus, we conclude:

Here, the constant $\unicode[STIX]{x1D6FE}$ corresponds to the net capillary force along $z$ if multiplied by the cross-sectional projection of the perimeter. It is given by the difference between surface-tension coefficients $\unicode[STIX]{x1D6FE}_{sa}$ and $\unicode[STIX]{x1D6FE}_{sl}$ for solid–air and solid–liquid interfaces, respectively. Equivalently, due to equilibrium of the contact line, $\unicode[STIX]{x1D6FE}$ is the same as the product of the liquid–air interfacial tension $\unicode[STIX]{x1D6FE}_{la}$ and the cosine of the contact angle $\unicode[STIX]{x1D703}$ . This means that

For low capillary number, $\unicode[STIX]{x1D6FE}$ does not vary, as the propagating front does not change its shape considerably during the transport process keeping the contact angle $\unicode[STIX]{x1D703}$ fixed.

The remaining scales correspond to the penetration length, the phenomenal duration and the velocity field. These are not only dependent on geometric and material parameters but also dictated by the period of the relevant transport process. Kinematic consideration concludes the velocity scale $v_{s}$ to be the ratio of characteristic penetration $h_{s}$ and time $t_{s}$ . However, we need two more constraints to ultimately determine $v_{s}$ , $h_{s}$ , $t_{s}$ . At this point, we recall that for rheometry, one has to resolve the timeframe in which transient effects are comparable to the dissipative influence of viscous friction. As a result, the characteristic values for both should be equated to the same for the pressure term induced by the capillary action. So one can derive the following conditions

From these it is inferred that

This implies that $t_{s}$ is defined as the relevant period for viscous dissipation in (2.9). The infiltration and its rate are characterized by $h_{s}$ and $v_{s}$ accordingly. For rheometric purposes, $t_{s}$ and $\unicode[STIX]{x1D70F}_{v}$ should be comparable. Under such conditions, however, the integral approximation calculating resistance from the steady-state profile would not be accurate enough.

The scales defined in (2.6) and (2.9) are used to normalize all quantities in (2.4) and (2.5). As a result, the non-dimensional governing relations become

where $\bar{v}_{z}=v_{z}/v_{s}$ , $\bar{t}=t/t_{s}$ , $\bar{p}=p/P_{s}$ , $\bar{h}=h/h_{s}$ are all dimensionless and $\bar{\unicode[STIX]{x1D735}}_{\Vert }$ is the cross-sectional gradient normalized by $b_{c}$ . The corresponding dimensionless interpretation of the operator $\hat{L}_{t}$ in (2.2) redefines it as

Here, $Wi$ is the Weissenberg number

which is a system-defining non-dimensional parameter denoting the ratio of the time scales for viscoelastic relaxation and the transient transport process.

2.3 Estimation of the pressure gradient

The capillary effect on the flow is manifested by suction at the propagating front. This creates a driving pressure gradient across the region of unidirectional streamlines

where the dimensionless suction pressure between two ends of the domain is $\unicode[STIX]{x0394}\bar{p}$ . We provide a proper model for $\unicode[STIX]{x0394}\bar{p}$ , and solve (2.10) to describe the transport process.

In past studies (Stange, Dreyer & Rath Reference Stange, Dreyer and Rath2003; Waghmare & Mitra Reference Waghmare and Mitra2010a ; Das et al. Reference Das, Waghmare and Mitra2012) on this topic, $\unicode[STIX]{x0394}\bar{p}$ is related to two separate pressure drops. The first of these appears outside the inlet of the channel where the pressure has to be slightly less than the ambient far away from the entry point to sustain a converging flow towards the conduit. Our problem, however, is fundamentally different from the systems considered in the earlier works where the capillary vessel sucks from an infinite body of fluid without exposure to the atmosphere in close proximity. In contrast, we have a free-surface reservoir in the form of an infinite drop with its atmospheric interface situated very near to the entrance. As mentioned before, this distinction is well represented in the two diagrams in figure 1. Thus, in the present problem, there is no difference in elevation between the inlet and the atmosphere. Moreover, the large size of the drop ensures negligible Laplace pressure in its interior. So one can safely assume an ambient condition at the entrance of the channel. As a result, we can find $\unicode[STIX]{x0394}\bar{p}$ by only taking into account the second contribution in $\unicode[STIX]{x0394}\bar{p}$ coming from the region attached to the propagating front.

The region near the front contains axially varying flow with three-dimensional structure transforming the unidirectional profile into a flow near a propagating front. Thus, the three-dimensional unsteady Navier–Stokes equation should govern the hydrodynamic fields here. Still, the scales for time, cross-sectional dimension, velocity and pressure must be identical between the two adjacent regions with one-dimensional and three-dimensional streamlines. The only scale which changes between these two subdomains is the axial length scale, which is $b_{c}$ for the latter instead of $h_{s}$ defined as the characteristic penetrated distance in § 2.2. When this change is included in the dimensional analysis, it makes the pressure gradient and the convective acceleration predominant over temporal and viscous terms. Consequently, one finds a new dimensionless momentum equation for velocity $\bar{\boldsymbol{v}}$

The constant non-dimensional parameter $\unicode[STIX]{x1D6FF}$ represents the effective capillary number which is the ratio of the characteristic viscous and surface-tension forces relevant to the system

If $b_{c}$ is much greater than 100 nanometres, $\unicode[STIX]{x1D6FF}$ is a very small quantity for a wide range of fluids. We are typically interested in cross-sectional dimensions of a few millimetres. Then, the value of $\unicode[STIX]{x1D6FF}$ becomes less than $10^{-3}$ . Hence, we can safely identify $\unicode[STIX]{x1D6FF}$ to be a very small quantity. Such a fact does not affect the transient nature of the unidirectional flow, as the pressure gradient there becomes comparable to the viscous and temporal terms due to the rescaling of the axial scale from $b_{c}$ to $h_{s}$ . This, however, simplifies the analysis near the front in two ways. Firstly, it makes the viscous resistance at the boundary layer attached to any interface negligible, especially in the calculation of the leading-order pressure field. Secondly, we can assume the shape of the propagating front to be fixed, as a small value of $\unicode[STIX]{x1D6FF}$ implies a small capillary number and a much stronger capillary force compared to the viscous stress. Consequently, for leading-order analysis in $\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x0394}\bar{p}$ near the front should be obtained from a steady inviscid flow equation near an undeformed free surface. This is why all the past studies (Stange et al. Reference Stange, Dreyer and Rath2003; Waghmare & Mitra Reference Waghmare and Mitra2010a ; Das et al. Reference Das, Waghmare and Mitra2012) have correctly obtained the pressure drop across the front region by using conservation of linear momentum without considering any boundary layer dissipation or interfacial deformation. We import their conclusion into our present analysis.

Accordingly, we note that the cross-sectional average $\unicode[STIX]{x0394}\bar{p}_{f}$ of the non-dimensional Laplace pressure at the front must be $-1$ due to the normalizing scale. This nominal value should be modified according to linear momentum conservation at the subdomain attached to the front to find the pressure $\unicode[STIX]{x0394}\bar{p}$ at the transition point between the unidirectional and three-dimensional flow. The difference between $\unicode[STIX]{x0394}\bar{p}$ and $\unicode[STIX]{x0394}\bar{p}_{f}$ should be equated to a change in linear momentum flux across the region. When we substitute $\unicode[STIX]{x0394}\bar{p}_{f}$ as $-1$ in this integral momentum relation, we find the well-known expression (Waghmare & Mitra Reference Waghmare and Mitra2010a )

where $\text{d}\bar{A}=\text{d}A/A$ . When the expression for $\unicode[STIX]{x0394}\bar{p}$ is combined with (2.10) and (2.13), the final form of the governing relation can be derived

which has to be solved to compute the flow field, the transport rate and the time-dependent penetration length.

2.4 Eigenfunction expansion

To obtain the flow solution, we consider a transient eigenfunction expansion of the velocity $\bar{v}_{z}$ , so that the unsteady flow is expressed in a separable form:

Here, time-dependent amplitudes are $\unicode[STIX]{x1D6FC}_{i}$ , whereas time-invariant eigenfunctions $u_{i}$ only vary with $\bar{\boldsymbol{r}}_{\Vert }$ which is the cross-sectional coordinates normalized by $b_{c}$ . A Helmholtz equation corresponding to the following eigen value equation defines $u_{i}$

where $-\unicode[STIX]{x1D6FD}_{i}^{2}$ are the allowable eigen values for which $u_{i}$ can satisfy no-slip boundary conditions at the periphery of the vessel. Sturm–Liouville theorem confirms that the eigenfunctions $u_{i}$ for any arbitrary geometry should be orthogonal to each other. We normalize $u_{i}$ in such a way that they become orthonormal as well:

where $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta. These functions can be constructed for any cross-section of interest using standard techniques such as the finite volume technique (Patankar Reference Patankar1980).

Ideally, the summation in (2.18) is an infinite one. However, we truncate it at a finite but reasonably high limit so that the solution of $\bar{v}_{z}$ is properly represented in the spectral space. Then, a complete set of integro-differential equations is derived involving all relevant $\unicode[STIX]{x1D6FC}_{i}(\bar{t})$ as well as $\bar{h}(\bar{t})$ by considering two facts. Firstly, the linear momentum conservation relation in (2.17) is exploited to obtain the number, $N$ , of independent equations for the same number of eigen modes included in the analysis. Secondly, the integral mass conservation also provides a single relation involving $\bar{h}(\bar{t})$ and all $\unicode[STIX]{x1D6FC}_{i}(\bar{t})$ . This means that there are $(N+1)$ unknowns including $\unicode[STIX]{x1D6FC}_{1},\unicode[STIX]{x1D6FC}_{2},\ldots ,\unicode[STIX]{x1D6FC}_{N}$ as well as $\bar{h}(\bar{t})$ . These have to be evaluated by simultaneously solving the $(N+1)$ coupled equations.

In our derivation, we first derive the expression of $\unicode[STIX]{x0394}\bar{p}$ in terms of eigen amplitudes $\unicode[STIX]{x1D6FC}_{i}(\bar{t})$ by applying Parseval theorem in (2.16)

where $\unicode[STIX]{x1D702}_{i}$ is the cross-sectional integral of $u_{i}$

Then, $\bar{v}_{z}$ is replaced in the linear momentum conservation relation of (2.17) with the expansion (2.18). The resulting equality is multiplied by the $i$ th eigenfunction $u_{i}$ , and integrated over the cross-section. Then, equations (2.19)–(2.22) are used to obtain the following:

for $i=1,2,\ldots ,N$ . The remaining equation is obtained from mass conservation so that the rate of penetration can be expressed in terms of $\unicode[STIX]{x1D6FC}_{i}$ and $\unicode[STIX]{x1D702}_{i}$

We solve (2.23) and (2.24) simultaneously to find $\bar{h}(\bar{t})$ and relevant $\unicode[STIX]{x1D6FC}_{i}(\bar{t})$ .

3.1 Details of the numerical solution scheme

The computation requires determination of the geometric parameters representing the shape of the conduit. These are the eigenfunctions $u_{i}$ , eigen values $\unicode[STIX]{x1D6FD}_{i}^{2}$ and cross-sectional integrals $\unicode[STIX]{x1D702}_{i}$ . However, we recognize that we need two eigen parameters associated with the two-dimensional feature of a rectangular cross-section. Accordingly, we replace the subscript $i$ by two indices $n$ and $m$ . Thus, from now on, $u_{i}$ , $\unicode[STIX]{x1D6FD}_{i}$ and $\unicode[STIX]{x1D702}_{i}$ will be referred to as $u_{nm}$ , $\unicode[STIX]{x1D6FD}_{nm}$ and $\unicode[STIX]{x1D702}_{nm}$ , respectively. For a rectangular tube, these quantities are

and

Here, $\bar{x}$ and $\bar{y}$ are the $x$ and $y$ coordinates normalized by $b_{c}=A/s$ , whereas $\bar{L}_{x}$ and $\bar{L}_{y}$ are the non-dimensional lengths of the sides along the $x$ and $y$ axes. The aspect ratio of the cross-section is defined by $\unicode[STIX]{x1D706}$ which is the ratio of $\bar{L}_{x}$ and $\bar{L}_{y}$ . Although ideally an infinite number of eigenfunctions should be considered, we take the maximum value for both $m$ and $n$ to be 40. This ensures less than 0.01 % relative error due to the spectral convergence.

Subsequent computation considers the flow to be driven by the capillary action from an initially static condition with an initial prefilled length $h_{0}$ . So the simulation essentially solves a start-up problem where the initial conditions are taken as $\bar{h}(0)=\bar{h}_{0}=h_{0}/h_{s}$ and $\unicode[STIX]{x1D6FC}_{nm}(0)=0$ . Then, a fourth-order Runge–Kutta scheme is used to simultaneously compute $\bar{h}(\bar{t})$ and $\unicode[STIX]{x1D6FC}_{nm}(\bar{t})$ from (2.23) and (2.24).

In each time step, we apply a Simpson technique to calculate the integral in the operator $\hat{L}_{t}$ defined by (2.2) representing the effect of viscoelasticity. In this integral, the lower limit is replaced by 0 in place of $-\infty$ , as there has been no motion before the initial time. Then stored past values of $\unicode[STIX]{x1D6FC}_{nm}$ are multiplied with the memory function $M$ , and the product is numerically integrated to determine the viscoelastic stress in (2.23). The fluid is considered to be a Maxwell fluid so that $M$ is exponential in nature

It is to be noted that (3.4) satisfies both constraints for $M$ as shown in (2.3).

3.2 Temporal variation in penetration dynamics

The outlined solution scheme finds the rate of imbibition of the viscoelastic liquid by using the computed values of $\unicode[STIX]{x1D6FC}_{nm}$ at each time step in (2.24). The intrusion rate is presented as a temporal function in figure 2. The result is obtained for both square channels and slit pores. Moreover, the simulation considers two different prefilled lengths for each geometry, as well as two values of the Weissenberg number to illustrate the impact of viscoelasticity. For mildly viscoelastic fluids like dilute polymeric solutions, a typical Weissenberg number should be less than 1. If the Weissenberg number becomes more than 1, the medium might not be governed by a linear viscoelastic law. This is the rationale behind choosing the specific $Wi$ values for the plots shown in figure 2. For proper comparison, the quantities are also calculated for a purely viscous fluid as a base case.

Figure 2. Dimensionless encroachment rates versus normalized time for a square channel (a,c) and slit pore (b,d) with initial prefilled length 0.5 (a,b) and 1.0 (c,d). The Weissenberg number $Wi$ is considered to be 0.5 (solid lines) and 0.2 (dash-dot line), whereas the dotted lines represent the curves for a purely viscous fluid.

In figure 2, higher values of the Weissenberg number $Wi$ imply a stronger viscoelastic effect causing significant departure in the short-time dynamics from the behaviour of a purely viscous liquid. Initially, the resistive stress in a viscoelastic fluid is less than the same for a viscous medium. So, the transport rate for the former is higher than the latter in the initial phase. At a later time, however, the difference is reversed, as the flow is pushed back in overcompensation after the relaxation of the memory effect. Such a reversal is a common feature in systems with viscoelasticity (Felderhof Reference Felderhof2009). Further explanation for the temporal oscillation is provided in § 4 where the Weissenberg number perturbation is described.

We expect the intrusion rate to decrease with the prefilled length $h_{0}$ . This is because a higher value of $h_{0}$ means a higher resistance to the flow. Such behaviour is evident in figure 2. It is to be noted, however, that $h_{0}$ does not affect the temporal period over which the relative departure between the viscoelastic and the purely viscous systems occurs.

Curves for different $Wi$ eventually approach the one for a purely viscous fluid under the long-time limit. The leading-order decay characteristic is actually independent of  $Wi$ . However, the next-order asymptotic decrease does depend on $Wi$ . Details of this feature are described in our subsequent perturbation analysis in § 5.

We also note that the flow rate is higher in a square channel compared with a slit pore. This is an expected behaviour, as the former is more conducive than the latter as long as the hydraulic radius given by $2b_{c}$ is the same for both.

The computation also provides the penetrated length at each time step. We plot this as a temporal function in figure 3 for all the parametric values considered in figure 2.

Figure 3. Dimensionless encroachment length versus normalized time for a square channel (a,c) and slit pore (b,d) with initial prefilled length 0.5 (a,b) and 1.0 (c,d). Weissenberg number $Wi$ is considered to be 0.5 (solid lines) and 0.2 (dash-dot line), the dotted lines represent curves for a purely viscous fluid.

The plots reveal two regimes of temporal variation. In the initial time, the penetration increases linearly. However, under the long-time limit, it shows Washburn’s prediction $\bar{h}\sim \sqrt{\bar{t}}$ . These results are corroborated by recent experiments (Quere Reference Quere1997; Das et al. Reference Das, Waghmare and Mitra2012).

From figure 3, one can conclude that the penetrated length $h$ is always larger in systems with higher $Wi$ at any temporal point. Although the intrusion rate is overcompensated after the initial larger value in a viscoelastic fluid, the compensation is not strong enough to reverse the difference in $h$ . Also, $h(t)$ for fixed $Wi$ is higher for a larger prefilled length as well as for a square channel compared to a slit pore.

3.3 Postprocessing algorithm for rheometry

In the actual rheometric experiments, the rheological properties defined by the memory function $M$ have to be detected from the observed transiency in the penetration rate. This is the reverse calculation of what is presented in § 3.2, where time-dependent behaviour is computed for given $M$ . Thus, an effective measurement technique requires an inversion algorithm.

The first step for this reverse analysis is to consider the temporal Fourier transform of (2.4). The resulting relation is, then, convoluted with the spatial eigenfunctions $u_{nm}$ . As a result, after rearrangement, one finds

Here, $\unicode[STIX]{x1D714}$ is the Fourier frequency, whereas the transformations of the amplitudes, memory function and force density are given as

At this point, we consider (2.24) in dimensional form and conclude the Fourier transform ${\hat{j}}$ for penetration rate to be

When (3.5) and (3.7) are combined, the inversion relation is derived

If ${\hat{j}}$ and $\hat{f}$ are known, the only unknowns remaining in (3.8) are the nominal viscosity for steady state $\unicode[STIX]{x1D707}_{0}$ and the normalized frequency-dependent complex coefficient $\hat{\unicode[STIX]{x1D707}}$ . The postprocessing algorithm exploits (3.8) to determine these rheological properties.

From the recorded experimental data, ${\hat{j}}$ can be obtained directly from the Fourier transform of the dimensional version of the plots presented in figure 2. In contrast, $\hat{f}$ has to be determined iteratively due to the presence of the subdominant but non-zero convective terms in (2.16). Still, in the first approximation, one can consider $\hat{f}$ to be a simple Fourier transform of inverse of the penetration length $h$ representing the predominant capillary effect. Then one can solve (3.8) for different $\unicode[STIX]{x1D714}$ to find $\unicode[STIX]{x1D707}_{0}$ and $\hat{\unicode[STIX]{x1D707}}(\unicode[STIX]{x1D714})$ . For $\unicode[STIX]{x1D714}=0$ , normalized $\hat{\unicode[STIX]{x1D707}}$ is $1$ by definition. Hence, for such a steady-state $\unicode[STIX]{x1D714}$ value, equation (3.8) yields the actual value for $\unicode[STIX]{x1D707}_{0}$ . Once this is evaluated, $\hat{\unicode[STIX]{x1D707}}$ can be calculated for non-zero $\unicode[STIX]{x1D714}$ . When the complete description of the rheology is available, one can repeat the flow analysis presented earlier to compute the hydrodynamic fields. As a result, the small convective contribution can be included in a more accurate $\hat{f}$ . A subsequent solution of (3.8) would, next, render new values for the rheological properties. In this way, an iterative scheme can be devised to improve $\hat{f}$ and the frequency-dependent viscosity in successive iterations until the results converge. The smallness of the convection term ensures the fast convergence of the proposed algorithm.

If $\unicode[STIX]{x1D707}_{0}$ and $\hat{\unicode[STIX]{x1D707}}(\unicode[STIX]{x1D714})$ are determined, the frequency-dependent complex viscosity would be known for a mildly viscoelastic medium. Then, the memory function $M$ can be constructed from the Fourier transform of $\hat{\unicode[STIX]{x1D707}}$ in the spectral domain and the viscoelastic time scale can be detected using the first moment normalization in (2.3). Also, the viscoelastic coefficients $G^{\prime }$ and $G^{\prime \prime }$ are related to $\unicode[STIX]{x1D707}_{0}\hat{\unicode[STIX]{x1D707}}(\unicode[STIX]{x1D714})$ by known simple expressions. Thus, the outlined procedure can determine the desired viscoelastic coefficients from easily acquired experimental data. Such technology has the potential to provide an accurate and cost effective means for rheometric measurement.

4.1 Taylor series expansion of the viscoelastic term

The first step in the perturbation analysis is to represent the following viscoelasticity term in a power series of the Weissenberg number:

Here $M$ is taken as a considerably localized distribution for small Weissenberg number, which implies $t_{s}\!\gg \!\unicode[STIX]{x1D70F}_{v}$ . As a result, we can disregard any difference in the viscoelastic dynamics caused by the initial condition of the system.

For the Weissenberg number expansion, we consider the following change in variable

The new temporal variable $\unicode[STIX]{x1D701}$ represents the time relevant to the viscoelastic scale instead of the transport scale. This transforms the viscoelastic integral

which is an equivalent reinterpretation of the memory effect.

We note that a Taylor series expansion in $Wi$ yields

When (4.4) is substituted into (4.3), one finds

where,

is the $n$ th moment of the distribution $M(\unicode[STIX]{x1D701})$ with $M_{0}$ and $M_{1}$ both being 1 as per (2.3).

4.2 Perturbation expansion in the Weissenberg number

The leading-order effect of the viscoelasticity on the dynamics is captured by expanding the dependent variables of (2.23) and (2.24) in power series of $Wi$

Here, the superscript 0 represents the quantities corresponding to the purely viscous fluid. In contrast, the first-order corrections in these variables due to viscoelasticity are denoted by the superscript 1.

We replace the dependent variables in (2.24) by their respective expansions, and form the mass conservation relation for each order in $Wi$

where $k$ is either 0 or 1. When (2.23) is similarly modified using (4.5)–(4.7) and arranged according to the exponent of $Wi$ , the zeroth-order recovers the momentum equation for the purely viscous flow (Bhattacharya et al. Reference Bhattacharya, Azese and Singha2016)

Also the leading-order correction in momentum conservation due to viscoelastic stress is represented by

In our earlier paper (Bhattacharya et al. Reference Bhattacharya, Azese and Singha2016), we solve (4.8) and (4.9) to obtain the zeroth-order temporal variations given by $\bar{h}^{(0)}$ and $\unicode[STIX]{x1D6FC}_{i}^{(0)}$ . These known functions are used in (4.10) to determine the leading-order corrections $\bar{h}^{(1)}$ and $\unicode[STIX]{x1D6FC}_{i}^{(1)}$ .

4.3 Simplified mathematical model

We attempt to describe the corrections $\bar{h}^{(1)}$ and $\unicode[STIX]{x1D6FC}_{i}^{(1)}$ in two ways. Firstly, equations (4.8) and (4.10) are solved exactly by numerical means using the Runge–Kutta method assuming a start-up problem from rest. These precise solutions, however, are not adequate for a complete physical understanding so we also pursue a second approach, where necessary simplifications are assumed so that the results can be justified more intuitively. Also such alternative analysis can verify the simulation independently.

Accordingly, we first simplify (4.10) for $i>0$ by neglecting the quadratic terms involving the zeroth-order amplitudes $\unicode[STIX]{x1D6FC}_{j}^{(0)}$

The rationale behind this simplification is derived from our prior study on purely viscous flow where $\unicode[STIX]{x1D6FC}_{j}^{(0)}$ is seen to be much less than unity.

Next, we also observe $\unicode[STIX]{x1D6FD}_{i}$ is much larger than unity for any channel when $i>0$ . This means that $\unicode[STIX]{x1D6FC}_{i}^{(1)}$ for $i>0$ saturates quickly to its long-time solution due to the fast exponential decay of the contribution from the initial condition. Hence, for $i>0$ , we approximate the long-time solution of (4.11) to be the actual $\unicode[STIX]{x1D6FC}_{i}^{(1)}$

Such an approximation immediately reduces the ( $N+1$ )-variable problem into a system with two variables only: $h^{(1)}$ and $\unicode[STIX]{x1D6FC}_{0}^{(1)}$ .

The approximate solution in (4.12) is then included into the integral mass conservation relation (4.8) which after rearrangements yields

Here, $A(\bar{t})$ , $B(\bar{t})$ , $C(\bar{t})$ are temporal functions dependent on the zeroth-order quantities:

One can explicitly compute $A(\bar{t})$ , $B(\bar{t})$ , $C(\bar{t})$ from the known transient variation of $h^{(0)}$ and $\unicode[STIX]{x1D6FC}_{i}^{(0)}$ .

Finally, the simplified momentum relation (4.11) is considered for $i=0$ and $\unicode[STIX]{x1D6FC}_{0}^{(1)}$ is replaced by using (4.13). As a result, we derive the following linear second-order equation involving only one dependent variable $h^{(1)}$

Here, $G(\bar{t})$ , $H(\bar{t})$ , $F(\bar{t})$ are three time-dependent functions representing transient damping, restoring and forcing terms, respectively. These are numerically computed from their respective expressions

We present $G(\bar{t})$ , $H(\bar{t})$ , $F(\bar{t})$ , for different geometries and initial conditions in figure 4.

Figure 4. The quantities defined in (4.16) are plotted as functions of time $\bar{t}$ for square channels (solid lines) and slit pores (dashed lines) with $\bar{h}_{0}=1.0$ (thick lines) and $\bar{h}_{0}=0.5$ (thin lines).

In figure 4, the effective dissipation in the system is represented by the function $G$ . This is essentially the ratio between the damping factor and the mass, both of which increase proportionally with the intruded length. So $G$ is expected to be approximately a constant in time. This is evident in figure 4.

In contrast, the function $H$ manifests a restoring spring term due to the viscoelasticity of the liquid. For a gradual deceleration over a long time, the elastic effect in the medium becomes irrelevant. This happens because for time scales much larger than the duration of viscoelastic relaxation, the medium acts like a purely viscous fluid. So the restoration function $H$ should monotonically decay to 0, as revealed in figure 4.

The crucial impact of viscoelasticity is captured by the forcing term $F$ which underlines the discrepancy in stress due to the memory effect. The main contribution in $F$ comes from the leading-order flow acceleration represented by the temporal derivative of $\unicode[STIX]{x1D6FC}_{i}^{(0)}$ . For a purely viscous system, the fluid first accelerates from rest under the action of the surface tension. Then, after some time, the flow starts to decelerate due to enhanced resistance coming from the increased length of the penetrated column. Accordingly, we expect $F$ to change sign after an initial period of positive value. This is precisely what can be seen in figure 4. Such a reversal is the core reason explaining the oscillatory behaviour in figure 2.

4.4 Comparative studies between reduced mathematical models and exact simulations

We compute the leading-order deviation in the transport rate for a purely viscous case due to the viscoelastic effect in three independent ways. Firstly, we consider a very small Weissenberg number ( $Wi=0.01$ ) and solve (2.23) and (2.24) directly, as described in § 3. As a result, the departure between the results for mildly viscoelastic and purely viscous systems is calculated exactly so that first-order correction ${\dot{h}}^{(1)}$ in the transport rate can be obtained by dividing the difference by $Wi$ . Then, ${\dot{h}}^{(1)}$ is also evaluated by solving the perturbation relation given by (4.8) and (4.10) using a Runge–Kutta method assuming a start-up problem from rest. Thirdly and finally, the approximate formulation given by (4.15) is exploited in another Runge–Kutta scheme to determine  ${\dot{h}}^{(1)}$ .

Figure 5. First-order deviation in penetration rate due to the viscoelastic effect is obtained as a function of time in three different ways for a square channel (a,c) and slit pore (b,d) where $\bar{h}_{0}$ is 0.5 (a,b) and $\bar{1.0}$ (c,d). We first use an approximate reduced description in (4.15) (dashed lines), and then the perturbative relations (4.8) and (4.10) (dotted lines). Also, the same results are directly computed from the departure between mildly viscoelastic ( $Wi=0.01$ ) and purely viscous ( $Wi=0$ ) systems using (2.23) and (2.24) (solid dots).

The three sets of simulations describing ${\dot{h}}^{(1)}$ as a temporal function are presented in figure 5 for a comparative study. We consider two different geometries and two separate initial conditions in figure 5. All cases show that the direct simulation and the exact perturbation corroborate each other precisely. It also reveals that the approximate model given by (4.15) predicts values which are very close to the other two results. Thus, we can assert the correctness of our numerical simulation and the validity of both the perturbative theory as well as the approximate reduced model.

The validity of the reduced model allows us to use its simplicity for an intuitive explanation of the oscillating temporal behaviour seen in both figures 2 and 5. The cause of the oscillation can be attributed to the fluctuating forcing function $F(\bar{t})$ given in figure 4 and defined in (4.15)–(4.16). A positive $F$ makes ${\dot{h}}^{(1)}>0$ at the initial time. However, the subsequent reversal in $F$ induces an eventual reversal in ${\dot{h}}^{(1)}$ too. However, the latter happens after a certain delay which can be justified by considering the phase difference due to the inertial second-order term and the restorative zeroth-order term in (4.15). The lag is also understandable from the causality perspective.

The dynamics of the system seems to be overdamped, as ${\dot{h}}^{(1)}$ approaches 0 gradually in the long-time limit. The reason for this is obvious from the functions in figure 4. The restorative term $H(\bar{t})$ , representing the elastic effect in (4.15), decays with time whereas the dissipative term $G$ remains more or less constant. This means that the system has to correspond to an overdamped situation for a large temporal scale.

5.1 Rescaling of variables and perturbation expansion

The perturbation analysis considers a transport time scale $t_{o}$ which is much larger than $t_{s}$ defined in (2.9). If time is normalized by $t_{o}$ , the corresponding non-dimensional variable is defined as $\tilde{t}$ . The rescaled time is related to $\bar{t}$ by the following renormalization:

where the ratio between $t_{s}$ and $t_{o}$ is $\unicode[STIX]{x1D716}$ which is identified as a small quantity based on which the perturbation theory is formulated.

We also take into account the corresponding rescaling of the intruded length and the amplitudes of the eigenfunctions. Accordingly, it is assumed that the leading orders of the viscous dissipation and the capillary force balance each other. Consequently, the rescaled penetration length and the $n$ spectral amplitudes are respectively defined as

and

As pressure gradient due to surface tension is inversely proportional to penetration, the factor $\sqrt{\unicode[STIX]{x1D716}}$ has to be involved in the renormalized quantities in (5.2) and (5.3).

Then we use the following expansions in $\unicode[STIX]{x1D716}$ to identify the decay features in the dynamics

and

where the $i$ th-order quantities in the small parameter $\unicode[STIX]{x1D716}$ are $\tilde{h}^{(i)}$ and $\tilde{\unicode[STIX]{x1D6FC}}_{nm}^{(i)}$ , respectively. The relation between $\tilde{h}^{(i)}$ and the corresponding $i$ th-order decaying feature in $\bar{h}$ is represented by the following transformation

If the $i$ th-order solutions are calculated properly, the transformation in (5.6) should ensure that the $i$ th-order decay in $\bar{h}$ is independent of $\unicode[STIX]{x1D716}$ . This invariance is a feature which can be used as a verification for the mathematical treatment.

Finally, the integral operator in (4.5) is recast in an expansion of $\unicode[STIX]{x1D716}$ after modifying the expression in terms of the transformed variables

Expansions of $\hat{L}_{t}$ as well as the transformed quantities are substituted into (2.23) and (2.24) to obtain a series of relations corresponding to specific orders of $\unicode[STIX]{x1D716}$ . The leading-order and the next-order equations are the focus of the subsequent analysis.

5.2 The leading-order behaviour

The leading-order momentum equation in $\unicode[STIX]{x1D716}$ provides a simple algebraic relation

where we recall $M_{0}=1$ . When this relation is used in the leading-order mass conservation equation, a simple differential equation for $\tilde{h}^{(0)}$ can be derived

Here, $\bar{k}$ is a constant dependent on the cross-sectional geometry

which represents the cumulative viscous dissipation in the system.

The solution for (5.9) is

Consequently, the transformation in (5.6) gives the leading-order trend in $\bar{h}$

as well as the long-time penetration rate

Here, $\tilde{c}_{0}$ and $\bar{c}_{0}$ are two integration constants which are related to the short-time dynamics and are irrelevant in the leading-order asymptotic feature

and

It is to be noted that $\bar{h}^{(0)}$ and $\dot{\bar{h}}^{(0)}$ are independent of $\unicode[STIX]{x1D716}$ , confirming validity of the analysis.

It is evident that the leading-order decay does not involve the Weissenberg number $Wi$ , as it is the same result also given by Washburn for a purely viscous fluid (Washburn Reference Washburn1921). This means that this asymptotic behaviour is not affected by viscoelasticity. Such an observation can be intuitively justified by the fact that a viscoelastic liquid acts as a purely viscous one in gradual transport under a long-time limit much greater than the relaxation time $\unicode[STIX]{x1D70F}_{v}$ .

5.3 First-order correction in the decay characteristic

The effect of viscoelasticity on the long-time dynamics is manifested in the next-order correction in $\unicode[STIX]{x1D716}$ , as the corresponding term in the expansion of the memory integral involves $Wi$ in (5.7). To highlight this fact, we obtain the solutions for $\tilde{\unicode[STIX]{x1D6FC}}_{nm}^{(1)}$ and $\tilde{h}^{(1)}$ .

Accordingly, the next-order momentum equation in $\unicode[STIX]{x1D716}$ is isolated

where $k^{\prime }$ is a constant dependent on the cross-sectional geometry

We substitute (5.16) in the mass conservation of similar order to derive the governing equation for $\tilde{h}^{(1)}$

where $M_{0}$ and $M_{1}$ are replaced by unity, as per (2.3).

We transform (5.18) into a convenient form by using (5.9)

Finally, the solution for $\tilde{h}^{(1)}$ is obtained as

where $\tilde{c}_{1}^{e}$ is a integration constant dependent on the short-time dynamics and $k_{p}$ is a factor defined below

Thus, the long-time decaying feature of $\bar{h}^{(1)}$ is given by the following $\unicode[STIX]{x1D716}$ -invariant expression

where $k_{1}^{e}$ is a constant related to $\tilde{c}_{1}^{e}$ .

Figure 6. Difference between the actual intruded length and the leading-order asymptotic behaviour in (5.14) normalized by the corresponding first-order decay function in (5.22). The ratio is plotted as a function of the non-dimensional time for a square channel (a,c) and slit pore (b,d), where the prefilled length is 0.5 (a,b) and 1.0 (c,d). For each plot, the Weissenberg number is 0.5 (solid line) or 0.2 (dashed line).

For mutual corroboration of the simulation and the asymptotic theory, we subtract the leading-order behaviour in (5.14) from the computed values of $\bar{h}^{(0)}$ obtained from the direct simulation. The difference is then divided by the first-order solution given in (5.22). The resulting ratio is plotted as a function of time $\bar{t}$ in figure 6 for slit pores and square channels with different values of prefilled length $h_{0}$ and Weissenberg number $Wi$ . All curves saturate to 1 in the long-time limit. This feature provides mutual corroboration between the computational data and the perturbation theory.

The higher-order perturbation analysis also shows how viscoelasticity affects the subdominant components of the asymptotic decay. This effect is represented by the factor $k_{p}$ which is dependent on both the cross-sectional geometry and the viscoelastic condition. This is why Weissenberg number appears in the expression of $k_{p}$ in addition to its variation with geometric parameters.

In figure 7, we present the constant $k_{p}$ as a function of the aspect ratio $\unicode[STIX]{x1D706}$ . We consider three values of $Wi$ to show how $k_{p}$ shifts with viscoelastic condition. For a slit pore, under the limit $Wi\rightarrow 0$ , $k_{p}$ should be $7/45$ which is matched by the computed result.

Figure 7. The geometry-dependent constant $k_{p}$ defined in (5.21) is plotted as a function of the aspect ratio $\unicode[STIX]{x1D706}$ for $Wi=0.4$ (solid line), $Wi=0.2$ (dotted line) and $Wi=0.1$ (dash-dot line).