Parameters and Branching Auto-Pulses in a Fluid Channel with Active Walls

Abstract

We present numerical solutions of the semi-phenomenological model of self-propagating fluid pulses (auto-pulses) in the channel branching into two thinner channels, which simulates branching of a hypothetical artificial artery. The model is based on the lubrication theory coupled with elasticity and has the form of a single nonlinear partial differential equation with respect to the displacement of the elastic wall as a function of the distance along the channel and time. The equation is solved numerically using the 1D integrated radial basis function network method. Using homogeneous boundary conditions on the edges of space domain and continuity condition at the branching point, we obtained and analyzed solutions in the form of auto-pulses penetrating through the branching point from the thick channel into the thin channels. We evaluated magnitudes of the phenomenological coefficients responsible for the active motion of the walls in the model.

1. Introduction

The arterial systems are characterised by branching with a network of larger arteries splitting into smaller arteries which continue to bifurcate into arterioles and then into the capillaries. Therefore, in this paper we focus on simulation of pulses in branching channels. In our recent work [1] only non-branching (single) channels were considered; this is the major difference between the two studies. In addition, in the present paper we evaluate the empirical coefficients playing the key role in our model.

We emphasize that the model simulates an artificial, not a real, artery. Note that the recent years saw a remarkable progress in design and fabrication of artificial muscles. For example, Reference [2] describes the fluid-driven and origami-inspired artificial muscles.

Arterial blood flows are usually modelled by the classical Navier-Stokes equations. However, under certain conditions, non-Newtonian models are also used [3,4,5]. While some models consider arteries with rigid walls [6], most models assume the walls to be elastic [7]. Roberts [8] presented a simple argument why the elasticity of an artery is important: if arteries were not elastic then each pump of the heart would cause an immediate rise in blood pressure throughout the body.

Our analysis is based on the model [9], which further developed the ideas of Roberts of channels with active walls [8]. While Roberts considered the case of negligible viscous forces compared to inertia, which is relevant to wider channels, we will consider the case when inertia is negligible compared to viscous forces, which is relevant to narrow channels. In the previous paper of the current authors [1], a variety of pulse solutions is obtained for a single channel. In the current paper we consider branching channels and analyse the conditions under which a pulse can successfully propagate from a thick channel into a thin one through the branching point. Ultimately we aim to make our model self-consistent that is to have characteristics that gurantee the propagation of the pulses over the entire system of branching channels.

Strunin [9] considered the flow between infinite elastic walls, assuming symmetry with respect to the plane, $z=0$; hence it will suffice to analyse only half of the flow, $0<z<H\left(t\right)$, (see Figure 1). The mean flow is assumed to be in the x-direction, driven by the pressure gradient in that direction. Adopting the lubrication theory for the flow [10] we equate the pressure gradient to the viscous friction.

$$\frac{{\partial }^{2}v}{\partial z^2}=\frac{1}{\eta }\frac{\partial p}{\partial x},$$

(1)

where $v(x,z,t)$ is the flow velocity in the x direction, $p(x,t)$ is the pressure and $\eta$ the viscosity.

Assuming p to be independent of z. we integrate Equation (1) twice with respect to z,

$$v=\frac{1}{2\eta }\frac{\partial p}{\partial x}\left(z^2-H^2\right)+v(x,H,t).$$

(2)

The mass flux Q is

$$Q={\int }_0^{H}vdz=-\frac{H^3}{3\eta }\frac{\partial p}{\partial x}+v(x,H,t)H.$$

(3)

Define the displacement, $w(x,t)$, of the wall in the z direction from the neutral position, $w=H_0$, by

$$H=H_0+w.$$

(4)

Then the continuity equation becomes

$$\frac{\partial w}{\partial t}+\frac{\partial Q}{\partial x}=0.$$

(5)

Substituting (3) into (5) gives

$$\frac{\partial w}{\partial t}=\frac{\partial }{\partial x}\left[\frac{H^3}{3\eta }\frac{\partial p}{\partial x}-v(x,H,t)H\right].$$

(6)

As we can see, the Equation (6) links the displacement of the flow boundary, coinciding with the wall position, to the flow pressure. The elasticity theory [11,12] provides the reverse link from the pressure to the displacement

$$p=D\frac{{\partial }^{4}w}{\partial x^4}-\frac{\partial }{\partial x}\left(N\frac{\partial w}{\partial x}\right),$$

(7)

where

$$N=\frac{Eh}{1-{\nu }^2}\left[\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\right].$$

(8)

In (7) and (8) $u(x,t)$ is the wall’s displacement along the flow, D is the flexural rigidity of the wall, E is Young’s modulus, h is the thickness of the wall, $\nu$ is Poisson’s ratio and N the force caused by the displacements. Strunin [9] supposed that the walls exert extra pressure relative to (7),

$$p=D\frac{{\partial }^{4}w}{\partial x^4}-\frac{\partial }{\partial x}\left(N\frac{\partial w}{\partial x}\right)+p_1,$$

(9)

where $p_1$ depends on w. He postulated that $p_1$ is proportional to the 4th power of the vertical displacement,

$$p_1=-\alpha w^4,\alpha >0.$$

(10)

After further simplifying assumptions (for more details see Reference [9]), the model has the form

$$\begin{array}{c}\hfill \frac{\partial w}{\partial t}=\frac{D}{3\eta }\frac{\partial }{\partial x}\left[H^3\frac{{\partial }^{5}w}{\partial x^5}\right]-\frac{\alpha }{3\eta }\frac{\partial }{\partial x}\left[H^3\frac{\partial }{\partial x}\left(w^4\right)\right]+\beta \frac{\partial }{\partial x}\left(H^3{w}^5\right).\end{array}$$

(11)

Dynamically the term ${\partial }_x\left(H^3{\partial }_x^{5}w\right)$ has dissipative effect, the nonlinear term, $-{\partial }_x\left[H^3{\partial }_x\left(w^4\right)\right]$, represents excitation and the nonlinear term ${\partial }_x\left(H^3{w}^5\right)$ transfers the energy from the excitation to dissipation.

Assuming $w<<H_0$ and replacing H by $H=H_0+w\approx H_0$ in Equation (11) we get

$$\begin{array}{c}\hfill \frac{\partial w}{\partial t}=\frac{D{H}_0^3}{3\eta }\frac{{\partial }^{6}w}{\partial x^6}-\frac{\alpha H_0^3}{3\eta }\frac{{\partial }^2}{\partial x^2}\left(w^4\right)+H_0^3\beta \frac{\partial }{\partial x}\left(w^5\right).\end{array}$$

(12)

Now we non-dimensionalize Equation (12) using some relevant scales to the form

$$\frac{\partial w}{\partial t}=A\frac{{\partial }^{6}w}{\partial x^6}-B\frac{{\partial }^2}{\partial x^2}\left(w^4\right)+C\frac{\partial }{\partial x}\left(w^5\right),$$

(13)

where the coefficients A, B and C are non-dimensional.

2. Numerical Experiments

To solve Equation (13) numerically we used the one-dimensional integration radial basis function network (1D-IRBFN) method in conjunction with one-step Picard iteration (PI1) scheme [13]. This method is more efficient than the original IRBFN method presented by Mai-Duy and Tran-Cong [14]. The purpose of using integration instead of conventional differeniation to construct the Radial Basis Function (RBF) approximations is to improve the stability and accuracy of the numerical solution. In the present paper, we use the following notations: $\widehat{\left[\right]}$ for a vector/matrix $\left[\right]$ that is associated with a grid line and ${\left[\right]}_{\left(n\right)}$ to denote selected components of the vector $\left[\right]$. In this method, the highest-order derivative, 6th order, is approximated by radial basis functions. The lower-order derivatives and function itself are then obtained by integration.

$$\begin{array}{cc}\hfill \frac{{\partial }^{6}w(x,t)}{\partial x^6}& =\sum _{i=1}^N{u}_i\left(t\right)G_i\left(x\right)=\sum _{i=1}^N{u}_i\left(t\right)H_6^{\left(i\right)}\left(x\right),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \frac{{\partial }^{5}w(x,t)}{\partial x^5}& =\sum _{i=1}^N{u}_i\left(t\right)H_5^{\left(i\right)}\left(x\right)+c_1,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \frac{{\partial }^{4}w(x,t)}{\partial x^4}& =\sum _{i=1}^N{u}_i\left(t\right)H_4^{\left(i\right)}\left(x\right)+c_{1}x+c_2,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \frac{{\partial }^{3}w(x,t)}{\partial x^3}& =\sum _{i=1}^N{u}_i\left(t\right)H_3^{\left(i\right)}\left(x\right)+\frac{c_1}{2}{x}^2+c_{2}x+c_3,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \frac{{\partial }^{2}w(x,t)}{\partial x^2}& =\sum _{i=1}^N{u}_i\left(t\right)H_2^{\left(i\right)}\left(x\right)+\frac{c_1}{6}{x}^3+\frac{c_2}{2}{x}^2+c_{3}x+c_4,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \frac{\partial w(x,t)}{\partial x}& =\sum _{i=1}^N{u}_i\left(t\right)H_1^{\left(i\right)}\left(x\right)+\frac{c_1}{24}{x}^4+\frac{c_2}{6}{x}^3+\frac{c_3}{2}{x}^2+c_{4}x+c_5,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill w(x,t)& =\sum _{i=1}^N{u}_i\left(t\right)H_0^{\left(i\right)}\left(x\right)+\frac{c_1}{120}{x}^5+\frac{c_2}{24}{x}^4+\frac{c_3}{6}{x}^3+\frac{c_4}{2}{x}^2+c_{5}x+c_6,\hfill \end{array}$$

(20)

where x is the input vector, N is the number of nodes on the x-axis, ${\left\{u_i\left(t\right)\right\}}_{i=1}^N$ are RBF weights to be determined; ${\left\{G_i\left(x\right)\right\}}_{i=1}^N={\left\{H_6^{\left(i\right)}\left(x\right)\right\}}_{i=1}^N$ are known RBFs, for example, for the case of multiquadrics (MQ)

$$G_i\left(x\right)=\sqrt{{(x-c_i)}^2+{a_i}^2},$$

where $c_i$ and $a_i$ are the centre and width of the $i^{th}$ MQ-RBF, respectively. The set of centres is chosen to be the same as the set of collocation points and the RBF width is determined as $a_i=b{d}_i$, $b>0$ is a factor (presently $b=1$) and $d_i$ is the distance from the i-th centre to the nearest; $H_5^{\left(i\right)}\left(x\right)=\int H_6^{\left(i\right)}\left(x\right)dx$; $H_4^{\left(i\right)}\left(x\right)=\int H_5^{\left(i\right)}\left(x\right)dx$; $H_3^{\left(i\right)}\left(x\right)=\int H_4^{\left(i\right)}\left(x\right)dx$; $H_2^{\left(i\right)}\left(x\right)=\int H_3^{\left(i\right)}\left(x\right)dx$; $H_1^{\left(i\right)}\left(x\right)=\int H_2^{\left(i\right)}\left(x\right)dx$; $H_0^{\left(i\right)}\left(x\right)=\int H_1^{\left(i\right)}\left(x\right)dx$; ${\left\{c_i\right\}}_{i=1}^6$ the set of constants arising from the integration process. The new basis functions $H_5^{\left(i\right)}\left(x\right)$, $H_4^{\left(i\right)}\left(x\right)$, $H_3^{\left(i\right)}\left(x\right)$, $H_2^{\left(i\right)}\left(x\right)$ and $H_1^{\left(i\right)}\left(x\right)$ are obtained from integrating the multiquadrics as follows,

$$\begin{array}{cc}\hfill H_5^{\left(i\right)}\left(x\right)& =\frac{r}{2}\sqrt{r^2+{a_i}^2}+\frac{{a_i}^2}{2}ln\left|r+\sqrt{r^2+{a_i}^2}\right|,\hfill \\ \hfill H_4^{\left(i\right)}\left(x\right)& =\left(\frac{r^2}{6}-\frac{{a_i}^2}{3}\right)\sqrt{r^2+{a_i}^2}+\frac{{a_i}^{2}r}{2}ln\left|r+\sqrt{r^2+{a_i}^2}\right|,\hfill \\ \hfill H_3^{\left(i\right)}\left(x\right)& =\left(\frac{r^3}{24}-\frac{13{a_i}^{2}r}{48}\right)\sqrt{r^2+{a_i}^2}+\left(\frac{{a_i}^2{r}^2}{4}-\frac{{a_i}^4}{16}\right)ln\left|r+\sqrt{r^2+{a_i}^2}\right|,\hfill \\ \hfill H_2^{\left(i\right)}\left(x\right)& =\left(\frac{r^4}{120}-\frac{83{a_i}^2{r}^2}{720}+\frac{{a_i}^4}{45}\right)\sqrt{r^2+{a_i}^2}+\left(\frac{{a_i}^2{r}^3}{12}-\frac{{a_i}^{4}r}{16})ln\left|r+\sqrt{r^2+{a_i}^2}\right|\right),\hfill \\ \hfill H_1^{\left(i\right)}\left(x\right)& =\left(\frac{r^5}{720}-\frac{97{a_i}^2{r}^3}{2880}+\frac{113{a_i}^{4}r}{5760}\right)\sqrt{r^2+{a_i}^2}\hfill \\ & +\left(\frac{{a_i}^2{r}^4}{48}-\frac{{a_i}^4{r}^2}{32}+\frac{{a_i}^6}{384}\right)ln\left|r+\sqrt{r^2+{a_i}^2}\right|,\hfill \\ \hfill H_0^{\left(i\right)}\left(x\right)& =\left(\frac{r^6}{5040}-\frac{253{a_i}^2{r}^4}{33600}+\frac{593{a_i}^4{r}^2}{67200}-\frac{{a_i}^6}{1575}\right)\sqrt{r^2+{a_i}^2}\hfill \\ & +\left(\frac{{a_i}^2{r}^5}{240}+\frac{{a_i}^{6}r}{384}-\frac{{a_i}^4{r}^3}{96}\right)ln\left|r+\sqrt{r^2+{a_i}^2}\right|,\hfill \end{array}$$

where $r=x-c_i$. After discretization, Equations (14)–(20) can be written in a compact form as

$$\widehat{\frac{{\partial }^{6}w}{\partial x^6}}=\widehat{H_6^{\left(i\right)}}\widehat{\alpha },\widehat{\frac{{\partial }^{5}w}{\partial x^5}}=\widehat{H_5^{\left(i\right)}}\widehat{\alpha },\widehat{\frac{{\partial }^{4}w}{\partial x^4}}=\widehat{H_4^{\left(i\right)}}\widehat{\alpha },$$

$$\widehat{\frac{{\partial }^{3}w}{\partial x^3}}=\widehat{H_3^{\left(i\right)}}\widehat{\alpha },\widehat{\frac{{\partial }^{2}w}{\partial x^2}}=\widehat{H_2^{\left(i\right)}}\widehat{\alpha },\widehat{\frac{\partial w}{\partial x}}=\widehat{H_1^{\left(i\right)}}\widehat{\alpha },\widehat{w}=\widehat{H_0^{\left(i\right)}}\widehat{\alpha },$$

(21)

where

$$\widehat{H_6^{\left(i\right)}}=\left[\begin{array}{cccccccccc}{H}_6^{\left(1\right)}\left(x_1\right)& H_6^{\left(2\right)}\left(x_1\right)& \cdots & H_6^{\left(N\right)}\left(x_1\right)& 0& 0& 0& 0& 0& 0\\ H_6^{\left(1\right)}\left(x_2\right)& H_6^{\left(2\right)}\left(x_2\right)& \cdots & H_6^{\left(N\right)}\left(x_2\right)& 0& 0& 0& 0& 0& 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ H_6^{\left(1\right)}\left(x_N\right)& H_6^{\left(2\right)}\left(x_N\right)& \cdots & H_6^{\left(N\right)}\left(x_N\right)& 0& 0& 0& 0& 0& 0\end{array}\right],$$

$$\widehat{H_5^{\left(i\right)}}=\left[\begin{array}{cccccccccc}{H}_5^{\left(1\right)}\left(x_1\right)& H_5^{\left(2\right)}\left(x_1\right)& \cdots & H_5^{\left(N\right)}\left(x_1\right)& 1& 0& 0& 0& 0& 0\\ H_5^{\left(1\right)}\left(x_2\right)& H_5^{\left(2\right)}\left(x_2\right)& \cdots & H_5^{\left(N\right)}\left(x_2\right)& 1& 0& 0& 0& 0& 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ H_5^{\left(1\right)}\left(x_N\right)& H_5^{\left(2\right)}\left(x_N\right)& \cdots & H_5^{\left(N\right)}\left(x_N\right)& 1& 0& 0& 0& 0& 0\end{array}\right],$$

…,

$$\widehat{H_0^{\left(i\right)}}=\left[\begin{array}{cccccccccc}{H}_0^{\left(1\right)}\left(x_1\right)& H_0^{\left(2\right)}\left(x_1\right)& \cdots & H_0^{\left(N\right)}\left(x_1\right)& x_1^5/5!& x_1^4/4!& x_1^3/3!& x_1^2/2!& x_1& 1\\ H_0^{\left(1\right)}\left(x_2\right)& H_0^{\left(2\right)}\left(x_2\right)& \cdots & H_0^{\left(N\right)}\left(x_2\right)& x_2^5/5!& x_2^4/4!& x_2^3/3!& x_2^2/2!& x_2& 1\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ H_0^{\left(1\right)}\left(x_N\right)& H_0^{\left(2\right)}\left(x_N\right)& \cdots & H_0^{\left(N\right)}\left(x_N\right)& x_N^5/5!& x_N^4/4!& x_N^3/3!& x_N^2/2!& x_N& 1\end{array}\right],$$

where ${\left\{x_i\right\}}_{i=1}^N$ is the set of nodal points, $\widehat{w}={\left(w_1,w_2,w_3,\dots ,w_N\right)}^T$, $w_i=w(x_i,t)$, $\widehat{u}={\left(u_1,u_2,u_3,\dots ,u_N\right)}^T$ and $\widehat{c}={(c_1,c_2,c_3,c_4,c_5,c_6)}^T$. We denote $\widehat{\alpha }={(u_1,u_2,u_3,\dots ,u_N,c_1,c_2,\dots ,c_6)}^T$. Equation (13) is discretized with respect to both time and space variables. Firstly, the time interval $[0,T]$ is partitioned into N subintervals $[t^{\left(n\right)},t^{(n+1)}]$ of length $\mathsf{\Delta }t=T/N$ with $t^{\left(0\right)}=1$ and $t^{(N+1)}=T$. The temporal discretization is then accomplished by a time-stepping scheme, followed by the spatial discretization based on the IRBFN method. Among the many possible time-stepping schemes, the standard $\theta$-scheme [15], $0\le \theta \le 1$ is used in this work. The extreme cases $\theta =0$ and $\theta =1$ correspond to the well-known forward (fully explicit) and backward (fully implicit) Euler schemes, respectively. The scheme with $\theta =1/2$ is known as the (semi-implicit) Crank-Nicolson method which is second-order accurate. We re-write Equation (13) as

$$\frac{\partial w}{\partial t}=f\left(w\right),$$

(22)

where

$$\begin{array}{c}\hfill f\left(w\right)=A\frac{{\partial }^{6}w}{\partial x^6}-B\left[12w^2{\left(\frac{\partial w}{\partial x}\right)}^2+4w^3\frac{{\partial }^{2}w}{\partial x^2}\right]+5C{w}^4\frac{\partial w}{\partial x}.\end{array}$$

(23)

Applying the standard $\theta$-scheme with $\theta =1/2$ Equation (22) is discretized as

$$\frac{w^{\left(n+1\right)}-w^{\left(n\right)}}{\mathsf{\Delta }t}=\frac{1}{2}{\left[f\left(w\right)\right]}^{\left(n+1\right)}+\frac{1}{2}{\left[f\left(w\right)\right]}^{\left(n\right)}$$

(24)

where $▵t=t^{(n+1)}-t^{\left(n\right)}$ is the physical time step; and the superscripts $\left(n\right)$ and $(n+1)$ denote the previous and current physical time levels, respectively. Note that ${\left[f\left(w\right)\right]}^{\left(n+1\right)}$ in Equation (24) consists of nonlinear terms. The one-step Picard iteration uses the solution at the previous time level to linearize the nonlinear terms. Equation (24) becomes

$$\begin{array}{cc}\hfill & \frac{w^{\left(n+1\right)}}{\mathsf{\Delta }t}-\frac{1}{2}\left\{A{\left(\frac{{\partial }^{6}w}{\partial x^6}\right)}^{\left(n+1\right)}+5C{\left(w^4\right)}^{\left(n\right)}{\left(\frac{\partial w}{\partial x}\right)}^{\left(n+1\right)}\right.\hfill \\ \left.-B\left[12w^{\left(n\right)}{\left({\left(\frac{\partial w}{\partial x}\right)}^2\right)}^{\left(n\right)}{w}^{\left(n+1\right)}+4{\left(w^3\right)}^{\left(n\right)}{\left(\frac{{\partial }^{2}w}{\partial x^2}\right)}^{\left(n+1\right)}\right]\right\}=\frac{w^{\left(n\right)}}{\mathsf{\Delta }t}+\frac{1}{2}{\left[f\left(w\right)\right]}^{\left(n\right)}\hfill \end{array}$$

(25)

After inserting appropriate values of w and its derivatives using (21), Equation (25) is written as

$$\begin{array}{c}\left[\begin{array}{c}\frac{{\widehat{H}}_0}{\mathsf{\Delta }t}-\frac{1}{2}\left\{\begin{array}{c}A{\widehat{H}}_6+5C{\left({\widehat{H}}_0{\widehat{\alpha }}^{\left(n\right)}\right)}^4{\widehat{H}}_1\hfill \\ -B\left[12{\widehat{H}}_0{\widehat{\alpha }}^{\left(n\right)}{\left({\widehat{H}}_1{\widehat{\alpha }}^{\left(n\right)}\right)}^2{\widehat{H}}_0+4{\left({\widehat{H}}_0{\widehat{\alpha }}^{\left(n\right)}\right)}^3{\widehat{H}}_2\right]\hfill \end{array}\right\}\hfill \end{array}\right]{\widehat{\alpha }}^{(n+1)}=\hfill \\ \left[\begin{array}{c}\frac{{\widehat{H}}_0}{\mathsf{\Delta }t}+\frac{1}{2}\left\{\begin{array}{c}A{\widehat{H}}_6+5C{\left({\widehat{H}}_0{\widehat{\alpha }}^{\left(n\right)}\right)}^4{\widehat{H}}_1\hfill \\ -B\left[12{\widehat{H}}_0{\widehat{\alpha }}^{\left(n\right)}{\left({\widehat{H}}_1{\widehat{\alpha }}^{\left(n\right)}\right)}^2{\widehat{H}}_0+4{\left({\widehat{H}}_0{\widehat{\alpha }}^{\left(n\right)}\right)}^3{\widehat{H}}_2\right]\hfill \end{array}\right\}\hfill \end{array}\right]{\widehat{\alpha }}^{\left(n\right)}\hfill \end{array}$$

(26)

For simplicity, the above equation can be written as,

$$E_1{\widehat{\alpha }}^{(n+1)}={\mathrm{RHS}}_1$$

(27)

and the boundary conditions as

$$E_2{\widehat{\alpha }}^{(n+1)}={\mathrm{RHS}}_2.$$

(28)

The system of Equations (27) and (28) is solved simultaneously at each time step for ${\widehat{\alpha }}^{(n+1)}$ until the prescribed time T is reached. Numerical experiments are conducted on a computer with an Intel i7, 3.60 GHz processor and 16 GB RAM using Matlab-R2017b software. The following subsections present the numerical results.

2.1. Two-Channel Experiment

In this numerical experiment we consider a thick channel branching into two thin channels as shown in Figure 2. In the state of rest each of the thin channels has half the width of the thick channel.

The displacement of the channel walls satisfy the equations

$$\begin{array}{cc}\hfill \frac{\partial w_1}{\partial t}& =\frac{D}{3\eta }{H}_0^3\frac{{\partial }^6{w}_1}{\partial x^6}-\frac{\alpha }{3\eta }{H}_0^3\frac{{\partial }^2}{\partial x^2}\left(w_1^4\right)+\beta H_0^3\frac{\partial }{\partial x}\left(w_1^5\right),\hfill \end{array}$$

(29)

$$\begin{array}{c}\hfill \frac{\partial w_2}{\partial t}=\frac{D}{3\eta }{\left(\frac{H_0}{2}\right)}^3\frac{{\partial }^6{w}_2}{\partial x^6}-\frac{\alpha }{3\eta }{\left(\frac{H_0}{2}\right)}^3\frac{{\partial }^2}{\partial x^2}\left(w_2^4\right)+\beta {\left(\frac{H_0}{2}\right)}^3\frac{\partial }{\partial x}\left(w_2^5\right).\end{array}$$

(30)

where we assume that the parameters ${\alpha }_i$ (and ${\beta }_i$) may be different between the channels. The boundary conditions are chosen homogeneous, namely zero values of the function and its first two derivatives at the edges, $x=x_l$ (l stands for “left”) and $x=x_r$ (r stands for “right” ),

$$w_1\left(x_r\right)=0,\frac{\partial w_1}{\partial x}\left(x_r\right)=0,\frac{{\partial }^2{w}_1}{\partial x^2}\left(x_r\right)=0,$$

$$w_2\left(x_l\right)=0,\frac{\partial w_2}{\partial x}\left(x_l\right)=0,\frac{{\partial }^2{w}_2}{\partial x^2}\left(x_l\right)=0.$$

In the experiments we used $x_r=0$ and $x_l=-45$ (see Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10). For the branching (contact) point the following boundary conditions are set. The kinematic condition expresses continuity of the wall displacement,

$$w_1\left(x_c\right)=w_2\left(x_c\right),$$

(31)

where c stands for “contact”, see Figure 2. The next conditions ensure continuity of the pressure (see (7)),

$$\begin{array}{cc}\hfill & \frac{\partial w_1}{\partial x}\left(x_c\right)=\frac{\partial w_2}{\partial x}\left(x_c\right),\frac{{\partial }^2{w}_1}{\partial x^2}\left(x_c\right)=\frac{{\partial }^2{w}_2}{\partial x^2}\left(x_c\right),\hfill \\ \frac{{\partial }^3{w}_1}{\partial x^3}\left(x_c\right)=\frac{{\partial }^3{w}_2}{\partial x^3}\left(x_c\right),\frac{{\partial }^4{w}_1}{\partial x^4}\left(x_c\right)=\frac{{\partial }^4{w}_2}{\partial x^4}\left(x_c\right).\hfill \end{array}$$

(32)

We used $x_c=-25$ (see Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10). The last condition will ensure continuity of the mass flux. According to (3), the flux for the thick channel is

$$Q_1=-\frac{H_0^3}{3\eta }\frac{\partial p_1}{\partial x}+v_1(x,H_0,t)H_0.$$

(33)

As the width of the thin channel is $H_0/2$, from Equation (33) the mass flux for the thin channel is

$$Q_2=-\frac{{\left(H_0/2\right)}^3}{3\eta }\frac{\partial p_2}{\partial x}+v_2\left(x,\frac{H_0}{2},t\right)\frac{H_0}{2}.$$

(34)

The continuity of the flux requires

$$Q_1=2Q_2,$$

(35)

Therefore, using (33), (34) we get

$$\begin{array}{c}\hfill -\frac{{\left(H_0/2\right)}^3}{3\eta }\frac{\partial p_2}{\partial x}+v_2\left(x,\frac{H_0}{2},t\right)\frac{H_0}{2}=-\frac{H_0^3}{6\eta }\frac{\partial p_1}{\partial x}+v_1\left(x,H_0,t\right)\frac{H_0}{2}.\end{array}$$

(36)

Assuming continuity of the velocity in (36) we have

$$4\frac{\partial p_1}{\partial x}=\frac{\partial p_2}{\partial x}.$$

(37)

Substituting (9) into (37), we obtain

$$\begin{array}{cc}\hfill & 4\left[D\frac{{\partial }^5{w}_1\left(x_c\right)}{\partial x^5}-\frac{{\partial }^2}{\partial x^2}\left(N\frac{\partial w_1\left(x_c\right)}{\partial x}\right)-\alpha \frac{\partial }{\partial x}{w}_1^4\left(x_c\right)\right]\hfill \\ =\left[D\frac{{\partial }^5{w}_2\left(x_c\right)}{\partial x^5}-\frac{{\partial }^2}{\partial x^2}\left(N\frac{\partial w_2\left(x_c\right)}{\partial x}\right)-\alpha \frac{\partial }{\partial x}{w}_2^4\left(x_c\right)\right],\hfill \end{array}$$

(38)

and, inserting (8) into (38),

$$\begin{array}{cc}\hfill & 4\left[D\frac{{\partial }^5{w}_1\left(x_c\right)}{\partial x^5}-\alpha \frac{\partial }{\partial x}{w}_1^4\left(x_c\right)\right.\left.-\frac{Eh}{1-{\nu }^2}\frac{{\partial }^2}{\partial x^2}\left(\frac{\partial u}{\partial x}\frac{\partial w_1\left(x_c\right)}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_1\left(x_c\right)}{\partial x}\right)}^3\right)\right]\hfill \\ =\left[D\frac{{\partial }^5{w}_2\left(x_c\right)}{\partial x^5}-\alpha \frac{\partial }{\partial x}{w}_2^4\left(x_c\right)\right.\left.-\frac{Eh}{1-{\nu }^2}\frac{{\partial }^2}{\partial x^2}\left(\frac{\partial u}{\partial x}\frac{\partial w_2\left(x_c\right)}{\partial x}+\frac{1}{2}{\left(\frac{\partial w_2\left(x_c\right)}{\partial x}\right)}^3\right)\right].\hfill \end{array}$$

(39)

Assuming u small and negligible the Equation (39) becomes

$$\begin{array}{cc}\hfill & 4D\frac{{\partial }^5{w}_1\left(x_c\right)}{\partial x^5}-4\alpha \frac{\partial }{\partial x}{w}_1^4\left(x_c\right)-\frac{4Eh}{2(1-{\nu }^2)}\frac{{\partial }^2}{\partial x^2}{\left(\frac{\partial w_1\left(x_c\right)}{\partial x}\right)}^3\hfill \\ =D\frac{{\partial }^5{w}_2\left(x_c\right)}{\partial x^5}-\frac{Eh}{2(1-{\nu }^2)}\frac{{\partial }^2}{\partial x^2}{\left(\frac{\partial w_2\left(x_c\right)}{\partial x}\right)}^3-\alpha \frac{\partial }{\partial x}{w}_2^4\left(x_c\right).\hfill \end{array}$$

(40)

Using the conditions (31) and (32) in Equation (40), we get

$$\begin{array}{c}\hfill 4D\frac{{\partial }^5{w}_1\left(x_c\right)}{\partial x^5}-3\frac{Eh}{2(1-{\nu }^2)}\frac{{\partial }^2}{\partial x^2}{\left(\frac{\partial w_1\left(x_c\right)}{\partial x}\right)}^3-3\alpha \frac{\partial }{\partial x}{w}_1^4\left(x_c\right)=D\frac{{\partial }^5{w}_2\left(x_c\right)}{\partial x^5}.\end{array}$$

(41)

Prior to the numerical experiment we non-dimensionalized all the equations (see Section 3 for details) and the values of A, B and C were chosen such that $A=1$, $B=1$ and $C=1$. The size of spatial domain in each numerical experiment is chosen such that: (1) the experiment does not last too long before a settled stage is reached and (2) the pulses are satisfactoryly resolved in space. The number of nodes is 100, the time step $0.001$ and the size of the domain $x_r-x_l=45$. In the first experiment, we assumed that ${\alpha }_2={\alpha }_1$ and ${\beta }_2={\beta }_1$. After non-dimensionalizing the thick channel equation to the form

$$\begin{array}{c}\hfill \frac{\partial w_1}{\partial t}=A\frac{{\partial }^6{w}_1}{\partial x^6}-B\frac{{\partial }^2}{\partial x^2}\left(w_1^4\right)+C\frac{\partial }{\partial x}\left(w_1^5\right)\end{array}$$

(42)

the non-dimensional thin channel equation becomes

$$\begin{array}{c}\hfill \frac{\partial w_2}{\partial t}=a\frac{{\partial }^6{w}_2}{\partial x^6}-b\frac{{\partial }^2}{\partial x^2}\left(w_2^4\right)+c\frac{\partial }{\partial x}\left(w_2^5\right),\end{array}$$

(43)

where $a=b=c=1/8$. We set the initial condition in the thick channel as $w_1(x,0)=1.2.exp\left[-0.25{(x+2.5)}^2\right]$ and in the thin channel $w_2(x,0)=0$ as shown in Figure 3.

From Figure 4, we see that in this particular experiment the pulse in the thick channel is not propagating to the thin channel. In the second experiment, in an attempt to make the pulse propagate into the thin channel we increased the value of ${\beta }_2$ (the thin channel parameter) by the factor of eight $(a=1/8,b=1/8,c=1).$

As a result, we observed the dynamics shown in Figure 5. The pulse propagates from the right end of the thick channel towards its left end as time goes. The evolution of the pulse starts from a single hump in the thick channel. After a while, the solution takes a familiar pulse-like form moving towards the contact point. The homogeneous boundary conditions were used at the left end of the thin channel and right end of the thick channel. At $t=12$, the pulse reaches the contact point $(x=-25)$.

Despite the pulse does make it into the thin channel and spends some time propagating through it, its height rapidly decreases and the pulse eventually decays without reaching the end of the channel. Looking at the mathematical structure of the governing equations, the only difference is the factor 1/8 in the right-hand side of the thin-channel equation compared to the thick-channel equation. Hence, the pulse solution definitely exists for the thin channel. But, as opposed to the way the pulse is started in the thick channel, within the thin channel we do not have freedom to choose the way the pulse is started, because it is already shaped by previous dynamics. In order to support the propagation in the thin channel we decided to increase the excitation coefficient ${\alpha }_2$ in the thin channel, so that the thin channel non-dimensional parameters $a=1/8,b=1/4,c=1$ (Figure 6). An interesting question is what is the minimum (critical) value of b in the thin channel to guarantee the pulse survival. We increased the excitation coefficient b in the thin channel from experiment to experiment by increments (Figure 6, Figure 7, Figure 8 and Figure 9) until the pulse survived. Observe from Figure 6 that at a relatively low excitation level in the thin channel the pulse only travels to the mark of just over 30. At a higher excitation rate, as illustrated in Figure 7, the pulse travels further, namely to the mark of over 35. For an even higher excitation, see Figure 9, the pulse manages to propagate through the entire length of the thin channel.

The pulse finally disappears after hitting the left boundary of the thin channel as shown in Figure 10.

We determined that the critical value of the excitation coefficient in the thin channel is approximately $b=0.6$.

2.2. Three-Channel Experiment

In this numerical experiment, we consider a three-channel configuration where a thick channel branches out into two thin channels and each of them in turn branches out into two even thinner channels. In the state of rest each of the thinner channels has half the width of the thicker channel with which it is in contact. The boundary conditions and the initial condition are the same as in the previous experiment. The displacement of the thinnest channel wall, $w_3$, satisfies the equation

$$\begin{array}{c}\hfill \frac{\partial w_3}{\partial t}=\frac{D}{3\eta }{\left(\frac{H_0}{4}\right)}^3\frac{{\partial }^6{w}_3}{\partial x^6}-\frac{{\alpha }_3}{3\eta }{\left(\frac{H_0}{4}\right)}^3\frac{{\partial }^2}{\partial x^2}\left(w_3^4\right)+{\beta }_3{\left(\frac{H_0}{4}\right)}^3\frac{\partial }{\partial x}\left(w_3^5\right).\end{array}$$

(44)

Now we non-dimensionalize Equation (44) to the form

$$\begin{array}{c}\hfill \frac{\partial w_3}{\partial t}=a_1\frac{{\partial }^6{w}_3}{\partial x^6}-b_1\frac{{\partial }^2}{\partial x^2}\left(w_3^4\right)+c_1\frac{\partial }{\partial x}\left(w_3^5\right).\end{array}$$

(45)

where the non-dimensional coefficient $a_1=1/46,b_1=4.9/46,c_1=1/8$. We run the numerical experiments in the similar fashion to the two-channel experiments, until the pulse survived in the thinnest channel as shown in Figure 11.

Figure 12 shows that the pulse decays after hitting the left boundary of the thinnest channel.

We then determined that the critical value of the non-dimensional excitation coefficient in the thinnest channel is approximately $0.38$.

3. Evaluation of the Model Coefficients

In this section, we evaluate the coefficients $\alpha$ and $\beta$, which were introduced in the original model [9] empirically. The idea is to select $\alpha$ and $\beta$ in such a way that the dimensional amplitude and dimensional velocity of the pulse have realistic values, that are values typical for a human artery (by the order of magnitude). Let us re-scale the main dimensional Equation (12) in two steps. In step one, we non-dimensionalize the equation using $H_0$ as the spatial scale and $\eta /E$ as the time scale. Then, in step two, we re-scale w, t and x using the yet-to-be-determined non-dimensional scaling factors W, T and X, respectively. Thus,

$$w=H_0{w}_3,w_3=W{w}_4,$$

$$x=H_0{x}_1,x_1=X{x}_2,$$

$$t=\eta t_1/E,t_1=T{t}_2.$$

The resulting non-dimensional equation is

$$\frac{\partial w_4}{\partial t_2}=A\frac{{\partial }^6{w}_4}{\partial x_2^6}-B\frac{{\partial }^2}{\partial x_2^2}\left(w_4^4\right)+C\frac{\partial }{\partial x_2}\left(w_4^5\right),$$

(46)

where

$$A=\frac{DT}{3E{H}_0^3{X}^6},B=\frac{\alpha H_0^{4}T{W}^3}{3E{X}^2},C=\frac{\beta \eta H_0^{6}T{W}^4}{EX}.$$

(47)

In fact, (46) is the non-dimensional Equation (13) which we solved numerically.

In the numerical experiments we used

$$A=1,B=1,C=1.$$

(48)

Using the settled single-pulse solution in say the thinnest channel (before the pulse hit the boundary), we can measure the non-dimensional amplitude of the pulse $\mathsf{\Delta }{w}_4$ and the non-dimensional pulse speed $\mathsf{\Delta }{x}_2/\mathsf{\Delta }{t}_2$. Returning to the dimensional quantities, we have, for the dimensional pulse amplitude,

$$\mathsf{\Delta }w=H_{0}W\mathsf{\Delta }{w}_4$$

(49)

and, for the dimensional pulse speed,

$$v=\frac{\mathsf{\Delta }x}{\mathsf{\Delta }t}=\frac{H_{0}XE}{\eta T}\frac{\mathsf{\Delta }{x}_2}{\mathsf{\Delta }{t}_2}.$$

(50)

We require that the coefficients $\alpha$ and $\beta$ lead to realistic orders of magnitude for $\mathsf{\Delta }w$ and v in (49) and (50). Of course we remember that our hypothetical channels do not have realistic cylindrical shape but the methodology of evaluation of $\alpha$ and $\beta$ is applicable to any version of the model including possible future versions with cylindrical configuration. From the literature, by the order of magnitude, $v=4\mathrm{m}/\mathrm{s}$ [16,17] and $\mathsf{\Delta }w=1\mathrm{mm}$. The three Equations (48) (where A, B and C are given by (47)) and Equations (49) and (50) form a system of five equations with respect to the five parameters W, X, T, $\alpha$ and $\beta$ that are to be determined. We find

$$W=\left(\frac{\mathsf{\Delta }w}{\mathsf{\Delta }{w}_4}\right)H_0^{-1},$$

$$X=D^{1/5}{\left(3\eta \right)}^{-1/5}{v}^{-1/5}{H}_0^{-2/5}{\left(\frac{\mathsf{\Delta }{x}_2}{\mathsf{\Delta }{t}_2}\right)}^{1/5},$$

$$T=E{D}^{1/5}{H}_0^{3/5}{3}^{-1/5}{\left(v\eta \right)}^{-6/5}{\left(\frac{\mathsf{\Delta }{x}_2}{\mathsf{\Delta }{t}_2}\right)}^{6/5},$$

$$\alpha =H_0^{-12/5}{\left(\frac{\mathsf{\Delta }{x}_2}{\mathsf{\Delta }{t}_2}\right)}^{-4/5}{D}^{1/5}{\left(3\eta \right)}^{4/5}{\left(\frac{\mathsf{\Delta }{w}_4}{\mathsf{\Delta }w}\right)}^3{v}^{4/5},$$

$$\beta ={\left(\frac{\mathsf{\Delta }{x}_2}{\mathsf{\Delta }{t}_2}\right)}^{-1}{\left(\frac{\mathsf{\Delta }{w}_4}{\mathsf{\Delta }w}\right)}^{4}v{H}_0^{-3}.$$

As we noted, $\mathsf{\Delta }{x}_2/\mathsf{\Delta }{t}_2$ and $\mathsf{\Delta }{w}_4$ are measured from the numerical experiments (Figure 11): $\mathsf{\Delta }{x}_2/\mathsf{\Delta }{t}_2=0.4$ and $\mathsf{\Delta }{w}_4=1$ approximately. The other parameters have the following approximate values: $H_0=0.5\mathrm{cm}$, the Young’s modulus $E=3\times {10}^5\mathrm{Pa}$ [18], the Poisson ratio $\nu =0.5$ [19,20,21,22,23,24], the blood viscosity $\eta =0.004\mathrm{Pa}·\mathrm{s}$ [22] and the wall thickness $h=0.1\mathrm{cm}$ [24]. Therefore, the flexural rigidity $D=E{h}^3/\left[12(1-{\nu }^2)\right]=0.002\mathrm{Pa}·{\mathrm{cm}}^3$. Based on these figures, we have, by the order of magnitude,

$$\alpha =1.1\times {10}^4\mathrm{Pa}·{\mathrm{cm}}^{-4}$$

and

$$\beta =8\times {10}^7{\mathrm{cm}}^{-6}·{\mathrm{s}}^{-1}.$$

4. Conclusions

We applied the 1D-IRBF numerical method to solve the model of the flow between active walls adapted for a branching channel. We used homogeneous boundary conditions at the edges and continuity conditions at the branching (contact) point. We obtained and analyzed solutions in the form of auto-pulses penetrating through the branching point from the thick channel into the thin channel. The numerical experiments showed that the pulse decays while moving through the thin channel unless the excitation coefficient in the thin channel is increased. Therefore, we conducted a series of experiments using larger values of the excitation coefficient $\alpha$ in the thin channel until the pulse “survived.” The numerical results indicated that the thinner the channel the larger the excitation coefficient $\alpha$ needs to be in order to guarantee the pulse propagation. Based on the numerical experiments, we evaluated the empirical parameters $\alpha$ and $\beta$ responsible for the active component of the wall dynamics. So far we explored the branching configuration with only two branching points as maximum. In the future, it would be interesting to study a system of multiple vessels with three or more branching points.