2. Experiments

Two types of commercial filter paper (Whatman grades 5 and 6) were used in the experiments. To quantify the volume of water absorbed by the papers, paper strips measuring 5 mm in width and 5 mm in length were prepared. Each strip was vertically mounted on a support and then brought into contact with a water reservoir, allowing the strip to imbibe water via capillary action. The volume of water absorbed by a paper strip was determined by measuring the change in mass of the strip before and after saturation with water. To eliminate the influence of paper expansion on the measured volume of absorbed liquid, a similar experiment was conducted using silicone oil, which can be absorbed by paper without inducing volume expansion (Chang et al. Reference Chang, Seo, Hong, Lee and Kim2018). The masses and volumes of water and silicone oil absorbed by the paper strips, denoted as $m_w$ and $m_o$ , respectively, are listed in table 1. For all paper strips tested, the volume of absorbed water exceeded that of silicone oil. The densities of silicone oil and water were 0.91 and 1.0 g $\textrm {cm}^{-3}$ , respectively.

Table 1. Measured parameters for multilayered channels composed of grade 5 and grade 6 paper.

Figure 1. Schematic illustrations of (a) the fabrication process of a multilayered paper channel, and (b) the experimental set-up for measuring imbibition length through multilayered paper channels. (c) Sequential images showing water imbibition through a multilayered paper channel composed of grade 6 paper. The wet region appears darker between the black boundaries, indicating the advancing wetting front.

Figure 2. (a) Longitudinal section images of multilayered paper channels before and after water infiltration. Temporal variations in the volume of a single paper layer composed of (b) grade 5 paper and (c) grade 6 paper. In (b) and (c), $V_p$ and $V_p^\infty$ denote the increased volume of a single paper layer and its final increased volume after complete expansion, respectively. The circle and star symbols indicate the times when the paper layer reaches 10 % and 99 % of its final increased volume, respectively, and $t_c$ and $t_s$ are the times required to reach these values, respectively. The solid lines represent (3.1).

We fabricated multilayered paper channels using two types of commercial filter paper (Whatman grades 5 and 6), as illustrated in figure 1(a). To pattern the paper, we applied wax using a commercial wax printer (ColorQube 8870, Xerox), leaving a 3 mm wide line uncoated. The wax-patterned paper was then baked on a hot plate at $150\,^{\circ }\rm {C}$ for 2 minutes to allow the wax to penetrate the paper structure. Two layers of paper were laminated using double-sided adhesive tape (ST-930H, Daehyun ST) with height $H \approx 230$ ${\unicode{x03BC}}$ m, which served as a spacer (Channon et al. Reference Channon, Nguyen, Henry and Dandy2019). The thickness of the paper was 200 ${\unicode{x03BC}}$ m for grade 5, and 180 ${\unicode{x03BC}}$ m for grade 6.

We constructed an experimental set-up to visualise water flow through the multilayered paper channels, as depicted in figure 1(b). The multilayered paper channel was placed on a horizontal support to eliminate gravitational effects on water imbibition. Water imbibition was initiated by bringing the water reservoir into contact with the channel. The progression of the wetting front was recorded using a high-speed camera (Mini AX200, Photron). Figure 1(c) presents sequential images of the advancing wetting front during water infiltration through a grade 6 paper channel.

3. Experimental observation

We investigated water infiltration through multilayered paper channels by cutting the channels along their centreline in the longitudinal direction, and observing the flow using a high-speed camera (Mini AX200, Photron), as shown in figure 2(a). As water progressed through the channel, it rapidly infiltrated the pores of the paper layers above and below the channel. For both types of paper used, local infiltration occurred within 0.01 s. The absorption of water into the compressed cellulose fibre network induced expansion of the paper layers. The expansion followed shortly after the initial water absorption, proceeding rapidly thereafter. The swelling of the paper layers was predominantly inwards, with the final channel height reduced to approximately 10 % of its initial height for both types of paper.

Figure 2(b,c) show the temporal evolution of the volume of the paper channels. We define $V_p$ as the increased volume of a single paper layer, with $V_p^\infty$ representing the final increased volume after complete expansion. The volume ratio $V_p/V_p^\infty$ thus ranges from 0 to 1 over time. Initially, the volume ratio increased gradually, taking approximately 0.2 s to reach 0.1 for both cases. Following this initial phase, the rate of volume increase accelerated rapidly, approaching 1 within a short period.

Figure 3. Dependence of $l^2$ on $t$ for multilayered paper channels composed of (a) grade 5 and (b) grade 6 paper. The solid straight lines are fits to the data in the early stage ( $t \lt 0.4$ s). Insets show the experimental data on a logarithmic scale. At early times, the data follow the lower asymptotic line (dotted, slope $1/2$ ), while at later times, the imbibition length approaches the upper asymptotic line with the same slope.

Based on phenomenological observations of water absorption and paper swelling, we propose simple empirical equations to describe the time-dependent evolution of channel volume changes resulting from deformation of the cellulose fibre network. We define $t_c$ as the time required for the paper layer to reach 10 % of its final increased volume, indicated by circles in figure 2(b,c). For simplicity, we assume that water uptake and swelling proceed concurrently after the initial local infiltration, which occurs within approximately 0.01 s. Neglecting the initial 10 % swelling during the period $0 \lt t \lt t_c$ , the rapid volume increase for $t \gt t_c$ can be described by the exponential function (Lee et al. Reference Lee, Kim, Kim and Mahadevan2016)

(3.1) \begin{equation} \frac {V_p}{V_p^\infty }=1-\textrm{e}^{-(t-t_c)/\lambda } \quad (t\gt t_c), \end{equation}

where $\lambda$ is the time constant associated with paper swelling. Similarly, the volume of water absorbed by the expanded paper layer, $V_w$ , can be expressed as

(3.2) \begin{equation} \frac {V_w}{V_w^\infty }=1-\textrm{e}^{-(t-t_c)/\lambda } \quad (t\gt t_c), \end{equation}

where $V_w^\infty$ represents the total volume of water absorbed when the paper reaches full expansion. We further define $t_s$ as the time required for the volume of the paper layer to reach 99 % of its final value, marked by stars in figure 2(b,c), measured from the moment $t = t_c$ . For all paper grades, $t_s$ was approximately 0.1 s. The values of $t_c$ and $t_s$ for each paper type are listed in table 1. In figure 2(b,c), (3.1) is represented by solid lines, which reasonably capture the characteristic rapid increase in the volume of the paper layers (see § 6 for the determination of $\lambda$ ).

Figure 3 presents our experimental measurements of water imbibition length through the multilayered paper channels over a period of 3 s. During the early stage ( $t \lt 0.4$ s), the experimental data exhibit an approximately linear increase with respect to time, aligning well with the straight data-fitting lines. However, for $t \gt 0.4$ s, the experimental data deviate from the straight lines and curve upwards for all paper channels.

4. Theoretical analysis

We begin the theoretical analysis by constructing a geometric model that incorporates both absorption and swelling of the paper layers, as illustrated in figure 4. The side surfaces of the multilayered paper channel are neglected, as the channel width ( $\approx3$ mm) is much greater than its height ( $\approx230$ ${\unicode{x03BC}}$ m). As shown in figure 4(a), capillary-driven flow between two horizontal, impermeable plates can be analysed analogously to capillary flow in a cylindrical tube, the basis of the classical Washburn equation, since both are governed by the balance between capillary forces and viscous resistance. Balancing the capillary force per unit depth, $2\sigma \cos \theta$ , with the viscous resistance per unit depth, $12\,{\mu} l l' / H$ , yields the momentum equation $l' = \sigma\!H \cos \theta / (6\,{\mu} l)$ , where $H$ is the initial channel height. Accordingly, the imbibition length can be expressed as $l = k\sqrt {t}$ , where $k = (\sigma\!H \cos \theta / (3\,{\mu} ))^{1/2}$ .

Figure 4. Model geometries for a channel with (a) an impermeable boundary, (b) a boundary that absorbs water without swelling, and (c) the present model geometry accounting for both absorption and swelling of the paper layers.

Compared to the flow in figure 4(a), the channel boundaries in figure 4(c) can absorb water and swell, resulting in a channel height that varies along the longitudinal direction. We consider an infinitesimal segment of the channel of length $\textrm {d}x$ located at a distance $\xi$ from the channel inlet, as indicated by the dashed box in figure 4(c). Within this segment, the continuity equation for water is given by $h(\textrm {d}u/\textrm {d}x) = -2(v_{\textit{ab}} - v_{w\textit{all}})$ , where $u$ is the local average velocity of water in the longitudinal direction, $v_{\textit{ab}}$ is the velocity of water moving into the paper boundary due to absorption, and $v_{w\textit{all}}$ is the inward velocity of the channel boundary resulting from paper swelling. Integrating the continuity equation from $\xi$ to $l$ yields

(4.1) \begin{equation} u(\xi )= l'(t)+\int ^{l}_{\xi } \frac {2(v_{\textit{ab}}-v_{w\textit{all}})}{h}\, \textrm {d}x, \end{equation}

with $l'$ being $\textrm {d}l/\textrm {d}t$ , and $u(l)=l'(t)$ .

The height of the channel, $h$ , depends on the local volume expansion of the paper and is bounded below by the minimum height $H - \Delta H$ . The expanded volume of the paper layer, given by $w\,\textrm {d}x\,(H - h)$ , is equal to $2V_p^\infty (1 - \textrm{e}^{-(t - \tau - t_c)/\lambda })$ from (3.1), where $w$ is the channel width. Substituting $V_p^\infty = w\textrm {d}x\,(\Delta H/2)$ , the channel height $h$ is expressed as

(4.2) \begin{equation} h(\xi )=H-\Delta H(1-\textrm{e}^{-(t-\tau -t_c)/\lambda }), \end{equation}

where $\tau$ is the time satisfying $l(\tau )=\xi$ , thus implying that $t-\tau$ is the local swelling time. Since $\textrm {d}h/\textrm {d}t=-2v_{w\textit{all}}$ , $v_{w\textit{all}}$ is expressed as

(4.3) \begin{align} v_{w\textit{all}}= \begin{cases} 0 \quad (t-\tau \lt t_c), \\[3pt] \displaystyle \frac {\Delta H}{2\lambda }\, \textrm{e}^{-(t-\tau -t_c)/\lambda } \quad (t-\tau \gt t_c). \end{cases} \end{align}

Figure 5. Schematics of the water imbibition phases for (a) $t \lt t_a$ , (b) $t_a \lt t \lt t_c$ and (c) $t \gt t_c$ .

We now derive an expression for $v_{\textit{ab}}$ . To do so, we introduce the time scale $t_a$ , which characterises the initial local infiltration of water into the pores of the paper layer. Upon contact with water, the adjacent paper layers are first saturated as water fills their pore spaces. During this initial phase ( $t \lt t_a$ ), the velocity of water entering the paper is given by $a\phi / t_a$ , where $a$ and $\phi$ denote the thickness and porosity of the paper, respectively. After $t_a$ , additional water absorption occurs, accompanied by swelling of the paper layer. The rate of water volume absorbed within the segment is given by $(w\,\textrm {d}x) v_{\textit{ab}} = \textrm {d}V_w/\textrm {d}t$ . Combining this relation with (3.2) and $V_w^\infty = (w\,\textrm {d}x)(\Delta H/2)$ , we obtain the following expression for $v_{\textit{ab}}$ :

(4.4) \begin{align} v_{\textit{ab}}= \begin{cases} a\phi /t_a \quad (t-\tau \lt t_a), \\[3pt] \displaystyle 0 \quad (t_a\lt t-\tau \lt t_c), \\[3pt] \displaystyle \frac {\Delta H}{2\lambda } \frac {V_w^\infty }{V_p^\infty }\, \textrm{e}^{-(t-\tau -t_c)/\lambda } \quad (t-\tau \gt t_c). \end{cases} \end{align}

The momentum equation for water flow in the longitudinal direction is derived from the balance between capillary force and viscous resistance. Due to water absorption and the associated volume expansion of the paper layers, the flow speed within the channel varies along the longitudinal direction. The momentum equation is given as $2 \sigma \cos \theta = \int ^l_0 2 \varGamma\, \textrm {d} x$ , where $\varGamma =6\,{\mu} u/h$ is the local shear stress in the paper channel. One can simplify the momentum equation to

(4.5) \begin{equation} \int ^l_0 \frac {u}{h}\, \textrm {d}x = \frac {\sigma \cos \theta }{6\,{\mu} }. \end{equation}

A single governing equation for $l(t)$ is obtained by substituting (4.1) into (4.5):

(4.6) \begin{equation} l'(t)\int ^{l}_{0}\frac {H}{h}\, \textrm {d}x+\int ^{l}_{0}\frac {H}{h}\left (\int ^{l}_{\xi }\frac {2(v_{\textit{ab}}-v_{w\textit{all}})}{h}\,\textrm {d}x\right )\textrm {d}x=\frac {1}{2}k^2, \end{equation}

where $k=(\sigma\!H \cos \theta /(3\,{\mu} ))^{1/2}$ . Since the water flow involves two characteristic time scales, $t_a$ and $t_c$ , distinct forms of the coupled equations arise in each corresponding time regime. We first consider the case $t \lt t_a$ (see figure 5 a). Substituting $v_{\textit{ab}} = a\phi / t_a$ , $v_{w\textit{all}} = 0$ and $h = H$ into (4.6) yields

(4.7) \begin{equation} l(t)\,l'(t)+\frac {a\phi }{H t_a}\{l(t)\}^2=\frac {1}{2}k^2, \end{equation}

where $\textrm {d}x=l'(\tau )\,\textrm {d}\tau$ is used.

Figure 5(b) illustrates the time regime $t_a \lt t \lt t_c$ . In this period, the flow speed in the channel, $u(\xi )$ , must be described by two separate expressions for the regions $0 \lt \xi \lt \alpha$ and $\xi \gt \alpha$ , where $\alpha = l(t - t_a)$ represents the length of the paper that is fully saturated with water but has not yet undergone expansion. In the region $0 \lt \xi \lt \alpha$ , we have $v_{\textit{ab}} = 0$ , $v_{w\textit{all}} = 0$ and $h = H$ . For $\xi \gt \alpha$ , the conditions are $v_{\textit{ab}} = a\phi / t_a$ , $v_{w\textit{all}} = 0$ and $h = H$ . Substituting these expressions into (4.6) yields

(4.8) \begin{equation} l(t)\,l'(t)+\frac {a\phi }{H t_a}\left[\{l(t)\}^2-\{l(t-t_a)\}^2\right]=\frac {1}{2}k^2. \end{equation}

Finally, we consider the time regime $t \gt t_c$ , as shown in figure 5(c). In this domain, the flow speed in the channel, $u(\xi )$ , must be described by three different expressions over the intervals $0 \lt \xi \lt \beta$ , $\beta \lt \xi \lt \alpha$ and $\xi \gt \alpha$ , where $\beta = l(t - t_c)$ denotes the length of the region where the paper has begun to swell. In the region $0 \lt \xi \lt \beta$ , the governing quantities are given by $v_{\textit{ab}} = (\Delta H / (2\lambda )) (V_w^\infty / V_p^\infty )\, \textrm{e}^{-(t - \tau - t_c)/\lambda }$ , $v_{w\textit{all}} = (\Delta H / (2\lambda ))\, \textrm{e}^{-(t - \tau - t_c)/\lambda }$ and $h = H - \Delta H\,(1 - \textrm{e}^{-(t - \tau - t_c)/\lambda })$ . For $\beta \lt \xi \lt \alpha$ , we use $v_{\textit{ab}} = 0$ , $v_{w\textit{all}} = 0$ and $h = H$ . For $\xi \gt \alpha$ , the conditions are $v_{\textit{ab}} = a\phi / t_a$ , $v_{w\textit{all}} = 0$ and $h = H$ . Substituting these expressions into (4.6) yields

(4.9) \begin{align} &l'(t) \left \{ l(t)-l(t-t_c) + \int ^{t-t_c}_{0}\frac {H\, l'(\tau )}{H-\Delta H\,(1-\textrm{e}^{-(t-\tau -t_c)/\lambda })}\, \textrm {d}\tau \right \} \\ \nonumber &\quad{}+\frac {a\phi }{H t_a} \left \{l(t)-l(t-t_a)\right \} \\ \nonumber & \quad{}\times \left [ l(t)+l(t-t_a) - 2 l(t-t_c)+ \int ^{t-t_c}_{0} \frac {H\, l'(\tau )}{H-\Delta H\,(1-\textrm{e}^{-(t-\tau -t_c)/\lambda })}\, \textrm {d}\tau \right ] \\ \nonumber &\quad{}+\frac {\Delta H}{\lambda } \left (\frac {V_w^\infty }{V_p^\infty }-1\right ) \\ \nonumber & \quad{}\times \int ^{t-t_c}_{0} \frac {H\, l'(\tau )}{H-\Delta H\,(1-\textrm{e}^{-(t-\tau -t_c)/\lambda })} \left ( \int ^{t-t_c}_{\tau } \frac {l'(\tau )\, \textrm{e}^{-(t-\tau -t_c)/\lambda }}{H-\Delta H\,(1-\textrm{e}^{-(t-\tau -t_c)/\lambda })}\, \textrm {d}\tau \right ) \textrm {d}\tau \\ \nonumber &{}=\frac {1}{2}k^2. \end{align}

In our experiments, the local infiltration time $t_a$ , defined as the time required for the paper to become locally saturated with water prior to swelling, was measured to be less than 0.01 s. This value is approximately an order of magnitude smaller than $t_c \approx 0.2$ s for both types of paper used. Given this clear separation of time scales, we simplify the forms of (4.7), (4.8) and (4.9). In the limit $t_a \rightarrow 0$ , (4.7) vanishes, while (4.8) and (4.9) reduce to

(4.10) \begin{equation} l(t)\,l'(t) \left(1 + \frac {2a\phi }{H}\right) = \frac {1}{2}k^2 \quad (t\lt t_c) \end{equation}

and

(4.11) \begin{align} &l'(t) \left (1+ \frac {2 a \phi }{H} \right ) \left \{ l(t)-l(t-t_c) + \int ^{t-t_c}_{0}\frac {H\, l'(\tau )}{H-\Delta H\,(1-\textrm{e}^{-(t-\tau -t_c)/\lambda })}\, \textrm {d}\tau \right \} \nonumber\\ \nonumber &\quad{}+\frac {\Delta H}{\lambda } \left (\frac {V_w^\infty }{V_p^\infty }-1\right ) \nonumber\\ \nonumber & \quad{}\times \int ^{t-t_c}_{0} \frac {H\, l'(\tau )}{H-\Delta H\,(1-\textrm{e}^{-(t-\tau -t_c)/\lambda })} \left ( \int ^{t-t_c}_{\tau } \frac {l'(\tau )\, \textrm{e}^{-(t-\tau -t_c)/\lambda }}{H-\Delta H\,(1-\textrm{e}^{-(t-\tau -t_c)/\lambda })}\, \textrm {d}\tau \right ) \textrm {d}\tau \nonumber\\ &{}=\frac {1}{2}k^2 \quad (t\gt t_c), \end{align}

respectively.

Table 2. Dimensionless parameters used in the model.

We proceed by non-dimensionalising (4.10) and (4.11) using the characteristic time scale $t_c$ and length scale $l_c = k \sqrt {t_c}$ . Expressed in terms of the dimensionless variables $t^\ast = t / t_c$ and $l^\ast = l / l_c$ , the governing equation becomes

(4.12) \begin{equation} l^\ast (t^\ast )\,{l^\ast }'(t^\ast )\, (1 + 2\psi ) = \frac {1}{2} \quad (t^\ast \lt 1) \end{equation}

and

(4.13) \begin{equation} {l^\ast }'(t^\ast )\,(1+2\psi )\,f(t^\ast ) + \varepsilon (\omega -1)\, g(t^\ast )= \frac {1}{2} \quad (t^\ast \gt 1), \end{equation}

where

(4.14) \begin{equation} f(t^\ast ) = l^\ast (t^\ast )-l^\ast (t^\ast -1) + \int ^{t^\ast -1}_{0} \frac {{l^\ast }'(\tau ^\ast )}{1-\varepsilon (1-\textrm{e}^{-(t^\ast -\tau ^\ast -1)/\eta })}\textrm {d}\tau ^\ast \end{equation}

and

(4.15) \begin{equation} g(t^\ast ) = \int ^{t^\ast -1}_{0} \frac {{l^\ast }'(\tau ^\ast )}{1-\varepsilon (1-\textrm{e}^{-(t^\ast -\tau ^\ast -1)/\eta })} \left ( \int ^{t^\ast -1}_{\tau ^\ast } \frac {{l^\ast }'(\tau ^\ast ) \,\textrm{e}^{-(t^\ast -\tau ^\ast -1)/\eta } / \eta }{1-\varepsilon (1-\textrm{e}^{-(t^\ast -\tau ^\ast -1)/\eta })}\textrm {d}\tau ^\ast \right ) \textrm {d}\tau ^\ast. \end{equation}

Here, $\psi = a\phi / H$ represents the ratio of the void volume in the paper boundary to the volume of the main channel, $\varepsilon = \Delta H / H$ is the ratio of the reduced channel height to the initial height, $\omega = V_w^\infty / V_p^\infty$ denotes the ratio of the absorbed water volume in the expanded region of the paper layer to its total expansion volume, and $\eta = \lambda / t_c$ is the ratio of the swelling time constant to the characteristic delay time. The dimensionless parameters are summarised in table 2.

In the time regime $t^\ast \lt 1$ , (4.12) indicates that imbibition through the channel is governed solely by water loss through the channel boundaries. When $\psi = 0$ , corresponding to the case where no water is absorbed into the boundaries, (4.12) reduces to the classical Washburn equation, $l^\ast (t^\ast )\,{l^\ast }'(t^\ast ) = 1/2$ .

When wall swelling occurs for $t^\ast \gt 1$ , the paper boundaries begin to expand locally. Because water fills the pore space within the paper boundaries as it flows, the absorbed water volume cannot exceed the expanded volume of the boundaries. This is consistent with our observation that the paper boundaries expand more than the volume of water that they retain, resulting in $0 \lt \omega \lt 1$ . Comparing (4.13) with (4.12), we observe that $l^\ast (t^\ast )$ in the first term is replaced by $f(t^\ast )$ , and an additional term, $\varepsilon (\omega - 1)\, g(t^\ast )$ , is introduced. As the channel walls expand inwards, two competing effects influence the imbibition speed. First, the reduction in channel height due to swelling increases the viscous shear stress at the boundaries, thereby slowing the flow. Second, because the boundaries expand more than the volume of water they absorb, they compress the fluid within the channel and push it forwards, enhancing the imbibition speed. Therefore, the two effects represented by $f(t^\ast )$ and $\varepsilon (\omega - 1)\, g(t^\ast )$ respectively act in opposition to one another in determining the flow speed within the main channel.

It is worth noting several limiting cases of (4.13). When $\omega = 1$ , the boundary expands by exactly the same volume as the absorbed water, causing the term $\varepsilon (\omega - 1)\, g(t^\ast )$ to vanish. In this case, only the first effect remains, represented by $f(t^\ast )$ . The condition $\varepsilon = 0$ corresponds to a situation where the channel walls do not swell at all, as shown in figure 4(b), reducing (4.13) to the simpler form of (4.12). Likewise, as $\eta$ increases, the swelling of the walls slows down. In the limit $\eta \rightarrow \infty$ , (4.13) also converges to (4.12), representing a no-swelling condition (see Appendix A).

5. Numerical solution

Since the analytical solution to (4.12) is readily obtained as $l^\ast (t^\ast ) = ({t^\ast / (1 + 2\psi )})^{1/2}$ , we numerically solve the dimensionless governing equation (4.13) for given values of $\psi$ , $\varepsilon$ , $\omega$ and $\eta$ , as detailed in the following numerical procedure.

The dimensionless time domain $(0, t_L^\ast )$ is discretised into $n + 1$ nodes: $t_0^\ast, t_1^\ast, \ldots,t_n^\ast{} (= t_L^\ast )$ . We numerically solve the equation ${l^\ast }'(t^\ast )\,(1 + 2\psi )\,f(t^\ast ) + \varepsilon (\omega - 1)\,g(t^\ast ) = 1/2$ at $t_m^\ast$ , where $r + 1 \leqslant m \leqslant n$ , and $r$ is defined such that $t_r^\ast \lt 1$ and $t_{r+1}^\ast \gt 1$ . Given that ${l^\ast }'(t^\ast ) = \{1/2-\varepsilon (\omega -1)\,g(t^\ast )\} / \{(1+2\psi )\, f(t^\ast )\}$ , integration of the equation from $t_{m-1}^\ast$ to $t_m^\ast$ yields.

(5.1) \begin{equation} \int ^{l^\ast (t_m^\ast )}_{l^\ast (t_{m-1}^\ast )} \textrm {d}l^\ast = \int ^{t_m^\ast }_{t_{m-1}^\ast } \left[ \frac{\{1/2-\varepsilon (\omega -1)\,g(t^\ast )\} } {\{(1+2\psi )\, f(t^\ast )\}} \right]\, \textrm {d}t^\ast \end{equation}

The numerical integration results in a simple recursive relation $l^\ast (t_m^\ast ) \approx l^\ast (t_{m-1}^\ast ) + [ \{1/2-\varepsilon (\omega -1)\,g(t_m^\ast )\} / \{(1+2\psi )\, f(t_m^\ast )\} ]\, \delta t^\ast$ , where $\delta t^\ast = t_L^\ast / n$ . The values of $f(t_m^\ast )$ and $g(t_m^\ast )$ are computed via numerical integration of the discretised governing equations, given by

(5.2) \begin{equation} f(t_m^\ast ) = l^\ast (t_{m-1}^\ast ) + l^\ast (t_{m-r}^\ast ) + \left \{\sum ^{m-r}_{j=0} \frac {{l^\ast }'(t_j^\ast )}{1-\varepsilon (1-\textrm{e}^{-(t_m^\ast -t_j^\ast -t_r^\ast )/\eta })}\right \} \times (\delta t^\ast ) \end{equation}

and

(5.3) \begin{equation} g(t_m^\ast ) = \left [\sum ^{m-r}_{j=0} \left \{ \frac {{l^\ast }'(t_j^\ast )}{1-\varepsilon (1-\textrm{e}^{-(t_m^\ast -t_j^\ast -t_r^\ast )/\eta })} \left ( \sum ^{m-r}_{i=j} \frac {{l^\ast }'(t_i^\ast )\, \textrm{e}^{-(t_m^\ast -t_i^\ast -t_r^\ast )/\eta } / \eta }{1-\varepsilon (1-\textrm{e}^{-(t_m^\ast -t_i^\ast -t_r^\ast )/\eta })} \!\right ) \!\right \}\! \right ]\! \times (\delta t^\ast )^2. \end{equation}

Sequential calculations can be initiated from $t_{r+1}^\ast$ , yielding results for $l^\ast (t_{r+1}^\ast )$ to $l^\ast (t_L^\ast )$ .