Next Article in Journal
An Inertial Forward–Backward Splitting Method for Solving Modified Variational Inclusion Problems and Its Application
Previous Article in Journal
Distributed Finite-Time Cooperative Economic Dispatch Strategy for Smart Grid under DOS Attack
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Poiseuille-Type Approximations for Axisymmetric Flow in a Thin Tube with Thin Stiff Elastic Wall

by
Kristina Kaulakytė
1,
Nikolajus Kozulinas
1,
Grigory Panasenko
1,2,*,
Konstantinas Pileckas
1 and
Vytenis Šumskas
1
1
Institute of Applied Mathematics, Faculty of Mathematics and Informatics, Vilnius University, Naugarduko Str., 24, 03225 Vilnius, Lithuania
2
Institute Camille Jordan UMR CNRS 5208, University Jean Monnet, 23, Rue Dr Paul Michelon, 42023 Saint-Etienne, France
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(9), 2106; https://doi.org/10.3390/math11092106
Submission received: 23 February 2023 / Revised: 17 March 2023 / Accepted: 12 April 2023 / Published: 28 April 2023
(This article belongs to the Section Difference and Differential Equations)

Abstract

:
An asymptotic ansatz for the solution of the axisymmetric problem of interaction between a thin cylindrical elastic tube and a viscous fluid filling the thin interior of the elastic tube was recently introduced and justified by G. Panasenko and R. Stavre. The thickness of the elastic medium ( ε ) and that of the fluid domain ( ε 1 ) are small parameters with ε < < ε 1 < < 1 , while the scale of the longitudinal characteristic size is of order one. At the same time, the magnitude of the stiffness and density of the elastic tube may be large or finite parameters with respect to the viscosity and density of the fluid when the characteristic time is of order one. This ansatz can be considered as a Poiseuille-type solution for the fluid–structure interaction problem. Its substitution to the Stokes fluid–elastic wall coupled problem generates a one-dimensional model. We present a numerical experiment comparing this model with the solution of the full-dimensional fluid–structure interaction problem.

1. Introduction

The interaction between a fluid and an elastic shell arises in many real-life applications in engineering and biology. The corresponding mathematical model is represented by the fluid–structure interaction problem, coupling two media with different, often highly contrasting physical characteristics. These problems are of great interest for the mathematical community due to their connection with blood flow models, for example, papers [1,2,3] studying the existence of a solution for fluid–structure interaction problem, and papers [4,5,6] considering viscous flows in elastic tubes. The wall is modeled either by a two-dimensional shell or by a full-dimensional elasticity or visco-elasticity equation; the small parameter is the thickness of a wall with respect to its diameter or the diameter of the pipe compared to its length, and the stiffness of the wall is considered finite. The stiffness of the wall as a large parameter was considered in [7,8,9], where the fluid–plate and fluid–cylindric wall interactions were studied, and the wall was modeled by a special boundary condition. In [10,11,12] the wall was modeled by the elasticity equation in two or three dimensions as a thin plate or thin cylindric layer, but the fluid domain was not thin. The asymptotic behavior of the solution when both the thickness of the wall and the diameter of the pipe are small parameters, and the wall has high stiffness, was first considered in [13]. Note that this setting is crucial for understanding an analog of a Poiseuille flow in a tube with an elastic wall. Paper [13] was devoted to this setting. Clearly, different orders of stiffness of the wall lead to different reduced models, but the ansatz of the asymptotic solution does not change for the cases considered in the cited paper. This property allows us to introduce this ansatz as an approximate Poiseuille-type flow in a stiff elastic tube. Plugging this ansatz in to the variational formulation of the full-dimensional problem, one obtains a one-dimensional model of the flow, taking into consideration the elasticity of the wall. The passage of the fluid through such a pipe is governed by the given inflow flux and an outflow flux, which can now be different functions of time. Recall that, for the rigid tube, they should be the same, due to the fluid incompressibility condition. The main result is a numerical comparison of the full-dimensional model of the flow and the one-dimensional approximation for a test problem. The results show good accuracy of the one-dimensional model.
Alternative approaches can be found in [14] (see also [15,16,17]), where the one-dimensional model is derived by engineering methods from the mass conservation law and the impulse-momentum theorem. It is presented by a hyperbolic-type system of two equations relating the cross-section area A of the vessel and the average velocity u. This model is applied in the case when the viscous term is dominated by the convective term. The viscosity is introduced by an empirical friction. In contrast to this model, we derive asymptotically the one-dimensional model from the Stokes equation coupled with the elasticity equation for the wall, and it is applied in the case that the viscous term is dominant in the Stokes equation. Thus, the model presented below is derived using more rigorous mathematical tools. However, its limits of application correspond to the modest Reynolds numbers. The goal of the present study is to evaluate numerically, using computational tests, the difference between the solution to the asymptotically derived one-dimensional model and the solution to the three-dimensional Navier–Stokes equation set in a thin tube with elastic wall (see Figure 1). The main conclusion of the paper is that, for modest Reynolds numbers, the asymptotically derived one-dimensional model is in good agreement with the solution to the three-dimensional Navier–Stokes equation. This allows this 1D model to be used for networks of blood vessels.
Note that, although the main application of this model is related to hemodynamics, it can be applied to viscous flows in industrial installations, for example, for the thin tube structures in [18].

2. Formulation of the Three-Dimensional Problem and Its Asymptotic One-Dimensional Approximation

Let us recall the formulation of the problem and the main result of [13].
Consider a fluid occupying the interior of a thin cylindrical elastic tube (as shown in Figure 1) with the radius denoted by ε 1 . More precisely, ref. [13] studies a periodic direction problem in longitudinal space, where ε 1 is defined as the ratio between the radius of the fluid domain and the length of the period, taken to be equal to one. The second small parameter is the elastic tube thickness ε with
ε < < ε 1 < < 1 .
The fluid and the elastic parts of the domain are defined, respectively, by
C f = { ( x 1 , x 2 , x 3 ) R 3 / x 1 2 + x 2 2 < ε 1 2 , x 3 R } , C ε e = { ( x 1 , x 2 , x 3 ) R 3 / ε 1 2 < x 1 2 + x 2 2 < ( ε 1 + ε ) 2 , x 3 R } .
Supposing that the coupled problem is axisymmetric and denoting r = x 1 2 + x 2 2 , let us associate the following layers in the plane ( x 3 , r ) to C f and C ε e :
L f = { ( x 3 , r ) R 2 / x 3 R , r ( 0 , ε 1 ) } , L ε e = { ( x 3 , r ) R 2 / x 3 R , r ( ε 1 , ε 1 + ε ) } .
The boundaries of the fluid and elastic layers L f and L ε e are denoted by
F 0 = { ( x 3 , 0 ) / x 3 R } , F ε 1 = { ( x 3 , ε 1 ) / x 3 R } , F ε 1 + ε = { ( x 3 , ε 1 + ε ) / x 3 R } ,
with the second one representing the interface surface.
The characteristics of the elastic medium are represented by the variable density ρ ˜ e = ρ ˜ e ( ξ ) depending on the radial variable, and by the matrix-valued coefficients A ˜ i j = A ˜ i j ( ξ ) , i , j { 1 , 2 , 3 } , which depend on the Young’s modulus E ˜ = E ˜ ( ξ ) and on the Poisson’s ratio ν ^ = ν ^ ( ξ ) , with ξ = r ε 1 ε , r ( ε 1 , ε 1 + ε ) . Because the density and the Young’s modulus characterizing a solid phase may be very different (high contrast) from one material to the other, it is interesting to analyze the changes in the mathematical model with respect to these physical characteristics. For this purpose, we introduce two additional parameters related to ε :
ω ρ = ε α , ω E = ε β , with   α , β R .
The density of the elastic tube is supposed to be of order ω ρ and its Young’s modulus to be of order ω E , while the characteristic time, the dynamic viscosity and the density of the fluid are supposed to be scaled so that they are of order one. So, ρ ˜ e = ω E ρ e and E ˜ = ω E E , with ρ e and E of order one. With respect to the physical characteristics of the elastic medium, we consider the following assumptions imposed for physical reasons:
ω ρ ω E ; ε 1 ω E .
With respect to the magnitude of ω E , the following different cases are identified:
(1)
ε 1 ω E < < ε 1 ε 1 2 ;
(2)
ω E = ε 1 ε 1 2 ;
(3)
ε 1 ε 1 2 < < ω E < < ε 1 ε 1 3 ;
(4)
ω E = ε 1 ε 1 3 ;
(5)
ω E > > ε 1 ε 1 3 .
In most practical applications, the relation between the density of an elastic medium and its Young’s modulus is ω ρ < < ω E . However, in order to simplify the illustration, we start with the relation
ω ρ = ω E ;
and then we pass to the case
ω ρ < < ω E ,
and we present the results if there are major differences with respect to the case (7).
Some practical examples from biology and engineering are given in the introductory part of [13].
Let us now give the formulation of the problem. Consider the stiffness matrices A ˜ i j having the elements A ˜ i j = ( a ˜ i j k l ) 1 k , l 3 , which satisfy the relations
a ˜ i j k l = ω E a i j k l , a i j k l = E 2 ( 1 + ν ^ ) 2 ν ^ 1 2 ν ^ δ i k δ j l + δ i j δ k l + δ i l δ j k , with the properties:
(i)
a i j k l ( ξ ) = a k j i l ( ξ ) = a j i l k ( ξ ) , ( ) i , j , k , l { 1 , 2 , 3 } , ( ) ξ [ 0 , 1 ] ,
(ii)
( ) κ > 0 independent of ε , ε 1 such that i , j , k , l = 1 3 a i j k l ( ξ ) η j l η i k κ j , l = 1 3 ( η j l ) 2 , ( ) ξ [ 0 , 1 ] , ( ) η = ( η j l ) 1 j , l 3 , with η j l = η l j .
The functions ρ e , ν ^ , E are supposed to be piecewise-smooth, and ρ e satisfies the following condition:
( ) m ρ , M ρ   independent of   ε , ε 1 , 0 < m ρ < M ρ ,   such that   m ρ ρ e ( ξ ) M ρ , ( ) ξ [ 0 , 1 ] .
So, the wall can be a natural laminate, as is the case with blood vessels. The characteristics of the viscous fluid, independent of ε , ε 1 , are the positive constants ρ f and ν representing its density and its dynamic viscosity, respectively. They are supposed to be of order one.
In addition to the data ρ e , E , ν ^ (for the elastic medium) and ρ f , ν (for the viscous fluid), we also consider as data of the problem the forces g and f, which represent an external action on the elastic medium and on the fluid, respectively. In the case of blood vessels, it corresponds to the action of the surrounding tissues; see [15].
The mathematical model of the periodic, axisymmetric, time-dependent interaction between an incompressible, Stokes fluid and a cylindrical, stratified, elastic tube uses the following notation and expressions for: the linear elasticity operator, the divergence operator (for a vector-valued function and a symmetric tensor-valued function), the gradient operator (for a vector-valued function) and the velocity strain tensor with respect to the cylindrical coordinates:
L u · e 3 = x 3 ( λ + 2 μ ) u 3 x 3 + λ u r r + 1 r u r + r μ u 3 r + u r x 3 + μ r u 3 r + u r x 3 , L u · e r = x 3 μ u 3 r + u r x 3 + r λ u 3 x 3 + 1 r u r + ( λ + 2 μ ) u r r + 2 μ r u r r 1 r u r , div c u = u 3 x 3 + u r r + 1 r u r , div c S = S 33 x 3 + 1 r r ( r S r 3 ) e 3 + S 3 r x 3 + 1 r r ( r S r r ) S θ θ r e r , c u = u 3 x 3 0 u 3 r 0 1 r u r 0 u r x 3 0 u r r . D c ( u ) = 1 2 c u + ( c u ) T .
Here ( x 3 , r , θ ) are the cylindrical coordinates; e 3 and e r are the unit vectors of O x 3 and O r axes, respectively; u represents a vectorial function; S represents a tensorial axisymmetric function; and λ = E ν ^ ( 1 + ν ^ ) ( 1 2 ν ^ ) , μ = E 2 ( 1 + ν ^ ) . In a two-dimensional setting, taking into account only the “ x 3 ” and “r” components, we will use the notation
D c ( u ) = u 3 x 3 1 2 u 3 r + u r x 3 1 2 u 3 r + u r x 3 u r r .
The mathematical model describing the fluid–elastic structure interaction in a thin-walled elastic tube with periodicity condition along the tube and homogeneous initial conditions is given by
ω ρ ρ e 2 u t 2 ω E L u = ε 1 g in L ε e × ( 0 , T ) , ρ f v t 2 ν d i v c D c ( v ) + p = f div c v = 0 in L f × ( 0 , T ) , v r = 0 on F 0 × ( 0 , T ) , u 3 r + u r x 3 = 0 λ ( 1 ) u 3 x 3 + ( λ ( 1 ) + 2 μ ( 1 ) ) u r r + λ ( 1 ) ε 1 + ε u r = 0 on F ε 1 + ε × ( 0 , T ) , v = u t ν v 3 r + v r x 3 = ω E μ ( 0 ) u 3 r + u r x 3 p + 2 ν v r r = ω E λ ( 0 ) u 3 x 3 + ( λ ( 0 ) + 2 μ ( 0 ) ) u r r + λ ( 0 ) ε 1 u r on F ε 1 × ( 0 , T ) , u , v , p 1 periodic in x 3 , u ( 0 ) = u t ( 0 ) = 0 in L ε e , v ( 0 ) = 0 in L f .
The positive constant T appearing in this system gives the time interval of the problem, and it can be arbitrarily chosen, independently of ε , ε 1 .
The main result of [13] is the construction of a complete asymptotic expansion of the solution in all the cases cited above, the evaluation of the norms of the error for the truncated expansions, and the computation of the leading term of the expansion, which is some analogue of the Poiseuille flow for tubes with elastic walls. Classical Poiseuille flow is defined as a flow with a given pressure drop at the ends of a tube. Formally, the problem (11) does not include such a case, because of the periodicity condition along the tube. However, such a pressure drop can be created in half of a tube. Then, the leading term of the asymptotic expansion of the solution, independently of the considered cases, can be interpreted as a Poiseuille-type flow; the velocity of the fluid and the displacements in the walls are expressed via three basic functions depending on the normal variable and time only. One of them is the scaled average velocity Q (the averaged velocity is equal to 2 ε 1 2 Q ), the second is the longitudinal displacement of the wall-fluid interface w 3 , and the third is the leading term for the pressure, q:
v 3 ( x 3 , r , t ) = 4 ε 1 2 1 r 2 ε 1 2 Q ( x 3 , t ) + w 3 t ( x 3 , t ) + ε 1 2 4 1 r 2 ε 1 2 ρ f ν 2 w 3 t 2 ( x 3 , t ) + 3 w 3 t x 3 2 ( x 3 , t ) , v r ( x 3 , r , t ) = ε 1 3 r ε 1 ( 2 r 2 ε 1 2 ) Q x 3 ( x 3 , t ) ε 1 r 2 ε 1 2 w 3 t x 3 ( x 3 , t ) ε 1 3 16 r ε 1 2 r 2 ε 1 2 ρ f ν 3 w 3 t 2 x 3 ( x 3 , t ) + 4 w 3 t x 3 3 ( x 3 , t ) , p ( x 3 , r , t ) = q ( x 3 , t ) , u 3 ( x 3 , r , t ) = w 3 ( x 3 , t ) + ε r ε 1 ε ε 1 3 0 t 2 Q x 3 2 ( x 3 , θ ) d θ + ε 1 2 2 w 3 x 3 2 ( x 3 , t ) ν ω E 1 ε ε 1 0 r ε 1 ε 1 τ μ ( τ ) d τ 8 Q ( x 3 , t ) ρ f 2 ν 2 w 3 t 2 ( x 3 , t ) + 3 w 3 t x 3 2 ( x 3 , t ) , u r ( x 3 , r , t ) = ε 1 3 1 ε 0 r ε 1 ε 1 ε 1 + ε τ λ ( τ ) λ ( τ ) + 2 μ ( τ ) d τ 0 t Q x 3 ( x 3 , θ ) d θ ε 1 2 1 ε 0 r ε 1 ε 1 ε 1 + ε τ λ ( τ ) λ ( τ ) + 2 μ ( τ ) d τ + ε 0 r ε 1 ε λ ( τ ) λ ( τ ) + 2 μ ( τ ) d τ w 3 x 3 ( x 3 , t ) + ω E 1 ε 0 r ε 1 ε 1 τ λ ( τ ) + 2 μ ( τ ) d τ 2 ν ε 1 2 Q x 3 ( x 3 , t ) ν 2 w 3 t x 3 ( x 3 , t ) q ( x 3 , t ) .
Here, for the leading terms, we keep the same notation as for the exact solution.
Note that the leading term for pressure, q, is related to the scaled average velocity Q by
q x 3 ( x 3 , t ) + 16 ν Q ( x 3 , t ) = f 3 ,
so that there are only two independent basic functions of the leading term of the ansatz, and the radial displacement of the wall-fluid interface, w r , can be approximately calculated as
w r ( x 3 , t ) = ε 1 3 0 t Q x 3 ( x 3 , θ ) d θ ε 1 2 w 3 x 3 ( x 3 , t ) ,
and so,
w r t ( x 3 , t ) = ε 1 3 Q x 3 ( x 3 , t ) ε 1 2 2 w 3 t x 3 ( x 3 , t ) .
If we need a continuous approximation of the velocity at the interface, then we have to add the third order terms in the approximation of u r :
ε 1 3 16 ρ f ν 2 w 3 t x 3 ( x 3 , t ) + 3 w 3 x 3 3 ( x 3 , t ) .
Note that the classical notion of the Poiseuille flow in an axisymmetric thin cylinder corresponds to the parabolic shape of the normal velocity of the flow and vanishing radial velocity. It was first described by J.L.M. Poiseuille and J. Boussinesq [19,20], and it takes place only for the stationary Stokes or Navier–Stokes equations. In the time-dependent case with no slip condition at the boundary (absolutely rigid wall), its analogue is the so-called Womersley flow [21,22], when the normal velocity is just independent of the normal coordinate, and the radial velocity still vanishes. However, in the case of a flow in a tube with an elastic wall, the normal velocity depends on the normal coordinate because the cross-section changes over time. Therefore, a Poiseuille-type flow is quite different from a flow with a rigid wall. The first term in formula (12) 1 is the “trace” of the rigid wall case. If the longitudinal velocity of the wall vanishes, then this term coincides with the result of [23]; the relation (13) coincides with the earlier known relations for the pressure slope and the average velocity (Darcy law), and the relations for the cross-section area and the average velocity (mass conservation law).
Relations (12) and (13) with two free basic functions Q and w 3 can be used as an ansatz in the frame of the method of asymptotic partial decomposition of the domain for the fully three-dimensional problem for the Stokes or Navier–Stokes equations in a network of vessels with fluid–structure interaction at the interface with elastic walls. Namely, the weak formulation of the full problem can be restricted to a subspace of functions having the form (12) and (13) within all cylindrical parts of a network at some distance from the bifurcations, similarly to the approach proposed in simpler cases in [10,24]. For the Navier–Stokes equations, formal analysis shows that the same ansatz can be used, but relation (13) should be modified in order to take into account non-linear terms. This approach allows us to reduce the number of dimensions to one in all cylindrical parts of the network and reduces the computational costs and computer memory required.
Comparing the ansatz (12) to the one-dimensional model [14], one can check that the leading term of the velocity satisfies the first equation of the system of equations
A t + F x 3 = 0 ,
where A ( x 3 , t ) is the area of the cross-section at the point x 3 , and F ( x 3 , t ) is the flux at the same point. However, the second equation of [14] differs from the equation below. This is explained by a different approach to the derivation of the one-dimensional model and the absence of the convective term in the formulation of the fluid–structure interaction problem.
Note that in models [15,16,17], the first equation of the system is the same.

3. The Variational Framework of the Problem

The first part of this section is devoted to the presentation of the functional spaces we are dealing with. Because our functions are periodic with respect to x 3 , we introduce the 3 D periodicity domains for the fluid and for the elastic tube, respectively:
Ω f = { ( x 1 , x 2 , x 3 ) R 3 / x 1 2 + x 2 2 < ε 1 2 , x 3 ( 0 , 1 ) } , Ω ε e = { ( x 1 , x 2 , x 3 ) R 3 / ε 1 2 < x 1 2 + x 2 2 < ( ε 1 + ε ) 2 , x 3 ( 0 , 1 ) }
and the corresponding periodicity plane domains
D f = { ( x 3 , r ) R 2 / x 3 ( 0 , 1 ) , r ( 0 , ε 1 ) } , D ε e = { ( x 3 , r ) R 2 / x 3 ( 0 , 1 ) , r ( ε 1 , ε 1 + ε ) } .
The boundaries corresponding to D f and D ε e are denoted by
Γ 0 = { ( x 3 , 0 ) / x 3 ( 0 , 1 ) } , Γ ε 1 = { ( x 3 , ε 1 ) / x 3 ( 0 , 1 ) } , Γ ε 1 + ε = { ( x 3 , ε 1 + ε ) / x 3 ( 0 , 1 ) } .
Let us introduce the following weighted spaces for the fluid part of the domain:
( L r 2 ( D f ) ) 2 = { ψ : D f R 2 / D f r ψ 2 ( x 3 , r ) d x 3 d r < } , ( H r 1 ( D f ) ) 2 = { ψ ( L r 2 ( D f ) ) 2 / D f r | c ψ | 2 ( x 3 , r ) d x 3 d r < } , ( H 0 , r 1 ( D f ) ) 2 = { ψ ( H r 1 ( D f ) ) 2 / r ψ = 0   on   Γ ε 1 } , ( H r 2 ( D f ) ) 2 = { ψ ( H r 1 ( D f ) ) 2 / D f r | c 2 ψ | 2 ( x 3 , r ) d x 3 d r < } ,
where
| c 2 ψ | 2 = 2 ψ 3 x 3 2 2 + 2 ψ 3 r 2 2 + 2 2 ψ 3 x 3 r 2 + 2 ψ r x 3 2 2 + 2 ψ r r 2 2 + 2 2 ψ r x 3 r 2 + 1 r 2 ψ 3 r 2 + 2 ψ r x 3 2 + 3 ψ r r 1 r ψ r 2 .
For periodic in x 3 functions, we introduce the periodic weighted spaces: H r , p e r l ( D f ) is the space of functions defined on L f , 1 periodic in x 3 , with the restriction to any arbitrary rectangle ( a , b ) × ( 0 , 1 ) in H r l ( ( a , b ) × ( 0 , 1 ) ) , ( ) l N and L r , p e r 2 ( D f ) = { h : L f R / h / D f L r 2 ( D f ) , h / D i f = h / D f , ( ) i Z } , where D i f = { ( x 3 , r ) R 2 / x 3 ( i , i + 1 ) , r ( 0 , 1 ) } , i Z ( D 0 f = D f ) . The space of the traces of functions from H r , p e r l ( D f ) is denoted by H r , p e r l 1 / 2 ( D f ) and, if we are interested only of the properties on a subset Γ D f , we shall write H r , p e r l 1 / 2 ( Γ ) , for l integer, l 1 . The space H r , p e r 1 / 2 ( Γ ) represents the dual of the subspace of H r , p e r 1 / 2 ( D f ) containing functions that vanish on D f Γ .
For the elastic part of the domain, the classical Sobolev spaces are used: H l ( D ε e ) , see [25,26], (as in D ε e we have r > ε 1 , H r l ( D ε e ) = H l ( D ε e ) ). Additionally, we introduce the periodic spaces H p e r l ( D ε e ) , L p e r 2 ( D ε e ) defined in a similar way as for the fluid domain and the classical trace spaces. However, we write H r l ( D ε e ) , etc., to show that the integrands in the norm contain the factor r.
The spaces that contain some of the properties given by (11) for the displacement and for the velocity, as well as the expected regularity of these functions, are defined by
U = φ ( H r , p e r 1 ( D ε e ) ) 2 / 0 1 φ r ( x 3 , 1 ) d x 3 = 0 , , V = ψ ( H r , p e r 1 ( D f ) ) 2 / div c ψ = 0 , ψ r = 0 on Γ 0 , H U = φ H 1 ( 0 , T ; U ) / 2 φ t 2 L 2 ( 0 , T ; U ) , H V = ψ L 2 ( 0 , T ; V ) / ψ t L 2 ( 0 , T ; V ) .
In order to present the space containing the coupling condition, let us set for any φ defined on D ε e and belonging to a certain functional space H,
V φ = ψ V / ψ = φ on Γ ε 1
and then
S H = ( φ , ψ ) / φ H , ψ V φ .
Note that the space H should have enough regularity to give sense to the trace of φ on Γ ε 1 .
Remark 1. 
In the axisymmetric problem, associate to any arbitrary function ψ : L f ( o r L ε e ) R 2 the function Ψ : C f ( o r C ε e ) R 3 such that
Ψ ( x 1 , x 2 , x 3 ) i s   t h e   3 D   c a r t e s i a n   f o r m   o f   ψ ( x 3 , r ) f o r ( x 1 , x 2 , x 3 ) C f ( o r C ε e ) ,
with x 1 2 + x 2 2 = r . The function Ψ is expressed with respect to the Cartesian coordinates and unit vectors and standard calculations give ψ ( H r l ( D f ( o r D ε e ) ) ) 2 if and only if Ψ ( H l ( Ω f ( o r Ω ε e ) ) ) 3 , ( ) l { 0 , 1 , 2 } .
In [13], it is assumed that
ρ e , a i j k l L ( 0 , 1 ) , i , j , k , l { 1 , 2 , 3 } , g H 1 ( 0 , T ; ( L r , p e r 2 ( D ε e ) ) 2 ) , f H 1 ( 0 , T ; ( L r , p e r 2 ( D f ) ) 2 ) .
In the framework presented above, the variational formulation of system (11), studied in [12], is given by
Find ( u , v ) H U × H V , such that ω ρ d d t D ε e r ρ e u ( t ) t · φ + ω E a L ( u ( t ) , φ ) + ρ f d d t D f r v ( t ) · ψ + 2 V D f r D c ( v ( t ) ) : D c ( ψ ) = ε 1 D ε e r g ( t ) · φ + D f r f ( t ) · ψ ( φ , ψ ) S U , a . e . in ( 0 , T ) , u t = v in L 2 ( 0 , T ; ( H p e r 1 / 2 ( Γ 1 ) ) 2 ) , u ( 0 ) = u t ( 0 ) = 0 in ( L r , p e r 2 ( D ε e ) ) 2 , v ( 0 ) = 0 in ( L r , p e r 2 ( D f ) ) 2 ,
where a L , defined by
a L ( u , φ ) = D ε e r μ 2 u 3 x 3 φ 3 x 3 + u r r φ r r + u 3 r + u r x 3 φ 3 r + φ r x 3 + 2 u r r φ r r + λ div c u div c φ ,
is the bilinear form associated with the elasticity operator L. The results concerning the existence, the uniqueness and the regularity of the unknown functions are obtained in [12].

4. Modified Variational Formulation for the Numerical Setup

We will modify the boundary conditions at the ends of the tube. Instead of the periodic solution with respect to the variable x 3 , we introduce some given inflow and outflow, supposing that the elastic tube is clamped at the ends of the tube. Namely, the periodic boundary conditions (11) 10 are replaced by a given velocity v r = 1 4 ν ( ε 1 2 r 2 ) g i n ( t ) , v 3 = 0 for x 3 = 0 and v r = 1 4 ν ( ε 1 2 r 2 ) g o u t ( t ) , v 3 = 0 for x 3 = 1 , and by vanishing displacements u = 0 for x3 = 0 and x3 = 1. Here gin and gout are given differentiable functions vanishing at t = 0. On the right-hand side, g and f are taken to be equal to zero. So, the problem is
ω ρ ρ e 2 u t 2 ω E L u = 0 in L ε e × ( 0 , T ) , ρ f v t 2 ν d i v c D c ( v ) + p = 0 div c v = 0 in L f × ( 0 , T ) , v r = 0 on F 0 × ( 0 , T ) , u 3 r + u r x 3 = 0 λ ( 1 ) u 3 x 3 + ( λ ( 1 ) + 2 μ ( 1 ) ) u r r + λ ( 1 ) ε 1 + ε u r = 0 on F ε 1 + ε × ( 0 , T ) , v = u t ν v 3 r + v r x 3 = ω E μ ( 0 ) u 3 r + u r x 3 p + 2 ν v r r = ω E λ ( 0 ) u 3 x 3 + ( λ ( 0 ) + 2 μ ( 0 ) ) u r r + λ ( 0 ) ε 1 u r on F ε 1 × ( 0 , T ) , v r = 1 4 ν ( ε 1 2 r 2 ) g i n ( t ) , v 3 = 0 , u = 0 for x 3 = 0 , v r = 1 4 ν ( ε 1 2 r 2 ) g o u t ( t ) , v 3 = 0 , u = 0 for x 3 = 1 , u ( 0 ) = u t ( 0 ) = 0 in L ε e , v ( 0 ) = 0 in L f .
Respectively, the variational formulation is changed. Let us define the spaces
U ˜ = φ ( H r 1 ( D ε e ) ) 2 , V ˜ = ψ ( H r 1 ( D f ) ) 2 / div c ψ = 0 , ψ r = 0 on Γ 0 ,
and
H ˜ U ˜ = φ H 1 ( 0 , T ; U ˜ ) / 2 φ t 2 L 2 ( 0 , T ; U ˜ ) ,
H ˜ V ˜ = ψ L 2 ( 0 , T ; V ˜ ) / ψ t L 2 ( 0 , T ; V ˜ ) .
Find ( u , v ) H ˜ U ˜ × H ˜ V ˜ , such that v 3 = 1 4 ν ( ε 1 2 r 2 ) g i n ( t ) , v r = 0 , u = 0 for x 3 = 0 , v 3 = 1 4 ν ( ε 1 2 r 2 ) g o u t ( t ) , v r = 0 , u = 0 for x 3 = 1 , and ω ρ d d t D ε e r ρ e u ( t ) t · φ + ω E a L ( u ( t ) , φ ) + ρ f d d t D f r v ( t ) · ψ + 2 ν D f r D c ( v ( t ) ) : D c ( ψ ) = 0 a . e . in ( 0 , T ) ( φ , ψ ) S U ˜ , such that φ ε | x 3 = 0 ; 1 = 0 , ψ ε | x 3 = 0 ; 1 = 0 , u t = v in L 2 ( 0 , T ; ( H 1 / 2 ( Γ 1 ) ) 2 ) , u ( 0 ) = u t ( 0 ) = 0 in ( L r 2 ( D ε e ) ) 2 , v ( 0 ) = 0 in ( L r 2 ( D f ) ) 2 .
Multiplying (28) 4 by functions of the base η k in H 0 2 ( 0 , T ) with some coefficients c k , and summing up and varying c k , we obtain the equivalent integral equation.
Find ( u , v ) H ˜ U ˜ × H ˜ V ˜ , such that v 3 = 1 4 ν ( ε 1 2 r 2 ) g i n ( t ) , v r = 0 , u = 0 for x 3 = 0 , v 3 = 1 4 ν ( ε 1 2 r 2 ) g o u t ( t ) , v r = 0 , u = 0 for x 3 = 1 , and ω ρ 0 T D ε e r ρ e 2 u ( t ) t 2 · φ d t + 0 T ω E a L ( u ( t ) , φ ) d t + 0 T ρ f D f r t v ( t ) · ψ d t + 2 0 T ν D f r D c ( v ( t ) ) : D c ( ψ ) d t = 0 ( φ , ψ ) H ˜ U ˜ × H ˜ V ˜ , such that φ | x 3 = 0 ; 1 = 0 , ψ | x 3 = 0 ; 1 = 0 , and φ t = ψ in L 2 ( 0 , T ; ( H 1 / 2 ( Γ 1 ) ) 2 ) , φ ( 0 ) = φ ( T ) = φ t ( 0 ) = φ t ( T ) = 0 in ( L r 2 ( D ε e ) ) 2 , ψ ( 0 ) = ψ ( T ) = 0 in ( L r 2 ( D f ) ) 2 .
Here, 2 ε 2 Q is the average velocity, and 2 π ε 4 Q is the flux.
Equation (29) 4 can be rewritten in the following form:
ω ρ 0 T D ε e r ρ e ( 2 u 3 t 2 φ 3 + 2 u r t 2 φ r ) + ω E 0 T D ε e ( 2 μ r ( u 3 x 3 φ 3 x 3 + u r r φ r r ) + μ r ( u 3 r + u r x 3 ) ( φ 3 r + φ r x 3 ) + 2 μ r u r r φ r r + λ r ( u 3 x 3 + u r r + u r r ) ( φ 3 x 3 + φ r r + φ r r ) ) + ρ f 0 T D f r ( v 3 t ψ 3 + v r t ψ r ) + 2 ν 0 T D f r ( v 3 x 3 ψ 3 x 3 + 1 2 ( v 3 r + v r x 3 ) ( ψ 3 r + ψ r x 3 ) + v r r ψ r r ) = 0 .
Finally, we consider a projection of this variational formulation to the “ansatz space” of test functions having the form (12) with q given by (13), with arbitrary free basic functions Q and w 3 from C 0 ( [ 0 , T ] × [ 0 , 1 ] ) . The solution is sought in the closure with respect to the norm H ˜ U ˜ × H ˜ V ˜ of the space of functions having the form (12) with smooth base functions.
In the following, we consider a simplified ansatz, assuming that w 3 = 0 and taking only a part of the terms in (12). This assumption is often satisfied in the blood vessels, such that the vessel does not move very much in the direction of blood flow. Namely, we consider the following space of test functions:
v 3 ( x 3 , r , t ) = 4 ε 1 2 1 r 2 ε 1 2 Q ( x 3 , t ) , v r ( x 3 , r , t ) = ε 1 3 r ε 1 ( 2 r 2 ε 1 2 i g ) Q x 3 ( x 3 , t ) , u 3 ( x 3 , r , t ) = ε r ε 1 ε ε 1 3 0 t 2 Q x 3 2 ( x 3 , θ ) d θ 8 ν ω E 1 ε ε 1 0 r ε 1 ε 1 τ μ ( τ ) d τ Q ( x 3 , t ) , u r ( x 3 , r , t ) = ε 1 3 1 ε 0 r ε 1 ε 1 ε 1 + ε τ λ ( τ ) λ ( τ ) + 2 μ ( τ ) d τ 0 t Q x 3 ( x 3 , θ ) d θ + ω E 1 ε 0 r ε 1 ε 1 τ λ ( τ ) + 2 μ ( τ ) d τ 2 ν ε 1 2 Q x 3 ( x 3 , t ) + 16 ν 0 x 3 Q ( s , t ) d s ,
where Q H 2 ( 0 , T ; H 4 ( 0 , 1 ) ) vanishes with its first derivatives in x 3 for x 3 = 0 and x 3 = 1 , and vanishes with its first derivatives in t for t = 0 and t = T . Furthermore, we will replace the notation of the basic function Q for the space of test functions with R. Additionally, we assume that functions g i n and g o u t satisfy the conditions g i n T = g o u t T = 0 and 0 t g i n ( τ ) d τ T = 0 t g o u t ( τ ) d τ T = 0 , where · T = 0 T · d τ , and we seek a solution with Q satisfying the relations Q T = 0 t Q ( τ ) d τ T = 0 . The test functions are also considered, with R satisfying the same relations: R T = 0 t R ( τ ) d τ T = 0 .
After substituting (31) in (30), we obtain the following differential equation:
C 1 ˜ 4 Q ( x 3 , t ) x 3 4 + C 2 ˜ 2 Q ( x 3 , t ) t 2 + C 3 ˜ 2 Q ( x 3 , t ) x 3 2 + C 4 ˜ 4 Q ( x 3 , t ) x 3 2 t 2 + C 5 ˜ 0 x 3 0 θ 2 Q ( s , t ) t 2 d s d θ + C 6 ˜ 0 t 0 τ 6 Q ( x 3 , θ ) x 3 6 d θ d τ + C 7 ˜ 0 t 0 τ 2 Q ( x 3 , θ ) x 3 2 d θ d τ + C 8 ˜ Q ( x 3 , t ) + C 9 ˜ 0 x 3 0 θ Q ( s , t ) d s d θ + C ˜ 10 0 t 0 τ 4 Q ( x 3 , θ ) x 3 4 d θ d τ + C ˜ 11 Q ( x 3 , t ) t + C ˜ 12 3 Q ( x 3 , t ) x 3 2 t = 0 ,
where
C 1 ˜ = ω ρ ρ e ε 1 6 ε 3 ( 3 ε + 4 ε 1 ) 12 + 4 ω E μ ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε 1 4 ε ( 11 ε + 16 ε 1 ) 120 + 11 24 ν ε 1 8 ,
C 2 ˜ = 64 ω ρ ρ e ( ν ω E 1 ε ε 1 μ ) 2 ε ( 11 ε + 16 ε 1 ) 120 64 ω ρ ρ e ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε 1 2 ε ( 11 ε + 16 ε 1 ) 120 ,
C 3 ˜ = 1 4 ω ρ ρ e ε 1 6 ( ( ( 2 ε 2 + 4 ε 1 ε + 2 ε 1 2 ) ln 2 ( ε + ε 1 ε 1 ) + ( 2 ε 2 4 ε 1 ε 2 ε 1 2 ) ln ( ε + ε 1 ε 1 ) + ε 2 + 2 ε 1 ε ) ( λ λ + 2 μ ) 2 + ( ( 4 ε 2 8 ε 1 ε 4 ε 1 2 ) ln ( ε + ε 1 ε 1 ) + 2 ε 2 + 4 ε 1 ε ) λ λ + 2 μ + 2 ε 2 + 4 ε 1 ε ) 64 ω E ( 2 μ + λ ) ( ν ω E 1 ε ε 1 μ ) 2 ε ( 11 ε + 16 ε 1 ) 120 4 ω E ( 2 μ + λ ) ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε 1 4 ε + 4 ε 1 12 ε + 64 ω E μ ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε 1 2 ε ( 11 ε + 16 ε 1 ) 120 + ω E ( 2 μ + λ ) ( ω E 1 λ + 2 μ ) 2 ν 2 ε 1 4 ( ( 4 ε 1 ε + 2 ε 1 2 ) ln ( ε 1 + ε ) 3 ε 2 + ( 4 ε 1 ln ( ε 1 ) 2 ε 1 ) ε 2 ε 1 2 ln ( ε 1 ) ) + 16 ω E λ ν ω E 1 ε ε 1 μ ω E 1 ε λ + 2 μ ν ε 1 2 ( 4 ε 15 + 7 ε + 15 ε 1 60 ) ω E λ ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε 1 4 ν ε 1 6 ,
C 4 ˜ = 4 ω ρ ρ e ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε 1 4 ε ( 11 ε + 16 ε 1 ) 120 , C 5 ˜ = 256 ω ρ ρ e ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε ( 11 ε + 16 ε 1 ) 120 ,
C 6 ˜ = ω E ( 2 μ + λ ) ε 1 6 ε 3 ( 3 ε + 4 ε 1 ) 12 ,
C 7 ˜ = ω E ( 2 μ + λ ) ( λ ε 1 λ + 2 μ ) 2 ln ( ε + ε 1 ε 1 ) + ω E ( 2 μ + λ ) ε 1 6 ( 1 3 ln 3 ( ε + ε 1 ε 1 ) ( λ λ + 2 μ ) 2 ln 2 ( ε + ε 1 ε 1 ) λ λ + 2 μ + ln ( ε + ε 1 ε 1 ) ) ω E λ 2 ε 1 6 λ + 2 μ ( ln ( ε + ε 1 ε 1 ) 1 2 ln 2 ( ε + ε 1 ε 1 ) λ λ + 2 μ ) ω E λ 2 ε 1 6 λ + 2 μ ( ln ( ε + ε 1 ε 1 ) 1 2 ln 2 ( ε + ε 1 ε 1 ) λ λ + 2 μ ) ,
C 8 ˜ = 64 ω E ( 2 μ + λ ) ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε 1 2 ε + 4 ε 1 12 ε + 64 ω E μ ( ν ω E 1 ε ε 1 μ ) 2 ε + 4 ε 1 12 ε + 256 ω E μ ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε ( 11 ε + 16 ε 1 ) 120 + 16 ω E ( 2 μ + λ ) ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε 1 2 ε 2 ( ( 4 ε 1 ε + 2 ε 1 2 ) ln ( ε 1 + ε ) 3 ε 2 + ( 4 ε 1 ln ( ε 1 ) 2 ε 1 ) ε 2 ε 1 2 ln ( ε 1 ) ) 16 ω E λ ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε 1 2 + 8 ω E λ ν ω E 1 ε ε 1 μ ω E 1 ε λ + 2 μ ( 4 ε 15 + 7 ε + 15 ε 1 60 ) + 16 ν ε 1 4 ,
C 9 ˜ = 64 ω E ( 2 μ + λ ) ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε + 4 ε 1 3 ε + 64 ω E ( 2 μ + λ ) ( ω E 1 ε λ + 2 μ ) 2 ν 2 ε 2 ( ( 4 ε 1 ε + 2 ε 1 2 ) ln ( ε 1 + ε ) 3 ε 2 + ( 4 ε 1 ln ( ε 1 ) 2 ε 1 ) ε 2 ε 1 2 ln ( ε 1 ) ) 64 ω E λ ( ω E 1 ε λ + 2 μ ) 2 ν 2 ,
C ˜ 10 = ω E μ ε 1 6 ( ε ε 1 + ε 2 2 ) 1 2 ω E μ ε 1 6 ( ( ( 2 ε 2 + 4 ε 1 ε + 2 ε 1 2 ) ln ( ε 1 + ε ε 1 ) ε 2 2 ε ε 1 ) λ λ + 2 μ 2 ε 2 4 ε ε 1 ) 1 4 ω E μ ε 1 6 ( ( ( 2 ε 2 + 4 ε 1 ε + 2 ε 1 2 ) ln 2 ( ε + ε 1 ε 1 ) + ( 2 ε 2 4 ε 1 ε 2 ε 1 2 ) ln ( ε + ε 1 ε 1 ) + ε 2 + 2 ε 1 ε ) ( λ λ + 2 μ ) 2 + ( ( 4 ε 2 8 ε 1 ε 4 ε 1 2 ) ln ( ε + ε 1 ε 1 ) + 2 ε 2 + 4 ε 1 ε ) λ λ + 2 μ + 2 ε 2 + 4 ε 1 ε ) + ω E λ ε 1 6 ε 2 2 ( λ + 2 μ ) + 1 2 ω E λ ε 1 6 ( ( ( 2 ε 2 2 ε 1 2 ) ln ( ε 1 + ε ε 1 ) ε 2 2 ε 1 ε ) λ λ + 2 μ 2 ε 2 ) + ω E ε 1 6 λ 2 λ + 2 μ ε 2 2 ,
C ˜ 11 = 8 3 ρ f ε 1 6 , C ˜ 12 = 11 24 ρ f ε 1 8 .
Next, we will consider a shorter approximation for the solution:
v 3 ( x 3 , r , t ) = 4 ε 1 2 1 r 2 ε 1 2 Q ( x 3 , t ) , v r ( x 3 , r , t ) = ε 1 3 r ε 1 ( 2 r 2 ε 1 2 ) Q x 3 ( x 3 , t ) , u 3 ( x 3 , r , t ) = 8 ν ω E 1 ε ε 1 0 r ε 1 ε 1 τ μ ( τ ) d τ Q ( x 3 , t ) , u r ( x 3 , r , t ) = 2 ω E 1 ε ν 0 r ε 1 ε 1 τ λ ( τ ) + 2 μ ( τ ) d τ ε 1 2 Q x 3 ( x 3 , t ) .
We assume that μ and λ are constants; therefore, from (33), we have the following expressions:
v 3 ( x 3 , r , t ) = 4 ε 1 2 1 r 2 ε 1 2 Q ( x 3 , t ) , v r ( x 3 , r , t ) = ε 1 3 r ε 1 ( 2 r 2 ε 1 2 ) Q x 3 ( x 3 , t ) , u 3 ( x 3 , r , t ) = 8 ν ω E 1 ε ε 1 μ ( r ε 1 ε ( r ε 1 ) 2 2 ε 2 ) Q ( x 3 , t ) , u r ( x 3 , r , t ) = 2 ω E 1 ε ν λ + 2 μ ( r ε 1 ε ( r ε 1 ) 2 2 ε 2 ) ε 1 2 Q ( x 3 , t ) .
As the boundary conditions for the solution, we use the value of 2 ε 1 2 Q equal to the average velocity at the inlet and outlet, and the x 3 derivative is equal to zero.
Substituting (34) into the following integral identity:
ω ρ d d t D ε e r ρ e ( u 3 t φ 3 + u r t φ r ) + ω E D ε e ( 2 μ r ( u 3 x 3 φ 3 x 3 + u r r φ r r ) + μ r ( u 3 r + u r x 3 ) ( φ 3 r + φ r x 3 ) + 2 μ r u r r φ r r + λ r ( u 3 x 3 + u r r + u r r ) ( φ 3 x 3 + φ r r + φ r r ) ) + ρ f d d t D f r ( v 3 ψ 3 + v r ψ r ) + 2 ν D f r ( v 3 x 3 ψ 3 x 3 + 1 2 ( v 3 r + v r x 3 ) ( ψ 3 r + ψ r x 3 ) + v r r ψ r r ) = 0 ,
we obtain the fourth order differential equation describing the elastic wall Poiseuille equation (EWPE):
C 1 2 Q ( x 3 , t ) t 2 + C 2 4 Q ( x 3 , t ) x 3 2 t 2 + C 3 Q ( x 3 , t ) t + C 4 3 Q ( x 3 , t ) x 3 2 t + C 5 4 Q ( x 3 , t ) x 3 4 + C 6 2 Q ( x 3 , t ) x 3 2 + C 7 Q ( x 3 , t ) = 0 ,
where
C 1 = 8 ω ρ ρ e ν 2 ω E 2 ε 1 2 ε 3 ( 16 ε 1 + 11 ε ) 15 μ , C 2 = ω ρ ρ e ν 2 ω E 2 ε 1 4 ε 3 ( 16 ε 1 + 11 ε ) 30 ( λ + 2 μ ) 2 ,
C 3 = 8 ρ f ε 1 6 3 , C 4 = 11 ρ f ε 1 8 24 , C 5 = ω E μ ν 2 ω E 2 ε 1 4 ε 3 ( 16 ε 1 + 11 ε ) 30 ( λ + 2 μ ) 2 + 2 ν 11 ε 1 8 48 ,
C 6 = ω E ( 8 ( λ + 2 μ ) ν 2 ω E 2 ε 1 2 ε 3 ( 16 ε 1 + 11 ε ) 15 μ 2 + 2 λ ν 2 ω E 2 ε 1 4 ε ( 4 ε 1 + ε ) 3 ( λ + 2 μ ) 2 + ( λ 2 μ ) ν 2 ω E 2 ε 1 4 ( ( 4 ε 1 ε + 2 ε 1 2 ) ln ( ε 1 + ε ) 3 ε 2 + ( 4 ε 1 ln ( ε 1 ) 2 ε 1 ) ε 2 ε 1 2 ln ( ε 1 ) ) ε ( λ + 2 μ ) 2 + 4 λ ν 2 ω E 2 ε 1 3 ε 2 ( 15 ε 1 + 7 ε ) 15 μ ( λ + 2 μ ) + 32 λ ν 2 ω E 2 ε 1 3 ε 3 15 μ ( λ + 2 μ ) + λ ν 2 ω E 2 ε 1 4 ε 2 ( λ + 2 μ ) 2 + 2 λ ν 2 ω E 2 ε 1 3 ε 3 ( 16 ε 1 + 11 ε ) 15 μ ( λ + 2 μ ) ) + ω E 2 μ ν 2 ω E 2 ε 1 3 ε 2 ( 15 ε 1 + 7 ε ) 15 μ ( λ + 2 μ ) + ω E 2 μ ν 2 ω E 2 ε 1 3 ε 2 ( 15 ε 1 + 7 ε ) 15 μ 2 ν 19 ε 1 6 6 + 4 ν 4 ε 1 6 3 ,
C 7 = ω E 16 μ ν 2 ω E 2 ε 1 2 ε ( 4 ε 1 + ε ) 3 μ 2 + 2 ν 8 ε 1 4 .

5. Numerical Scheme for the Poiseuille Equation

In this section, we numerically solve the following boundary value problem for the elastic wall Poiseuille equation (EWPE) (36):
C 1 2 Q ( x 3 , t ) t 2 + C 2 4 Q ( x 3 , t ) x 3 2 t 2 + C 3 Q ( x 3 , t ) t + C 4 3 Q ( x 3 , t ) x 3 2 t + C 5 4 Q ( x 3 , t ) x 3 4 + C 6 2 Q ( x 3 , t ) x 3 2 + C 7 Q ( x 3 , t ) = 0 , Q ( 0 , t ) = g i n ( t ) , Q ( 1 , t ) = g o u t ( t ) , Q ( 0 , t ) x 3 = Q ( 1 , t ) x 3 = 0 , Q ( x 3 , 0 ) = Q ( x 3 , 0 ) t = 0 .
We introduce uniform meshes of points in space and time with step sizes equal to h and τ , respectively:
x 3 m : = m h , m = 0 , , M , x 3 M = 1 , t n : = n τ , n = 0 , , N , t N = T .
Then, the numerical solution to (37) at grid point ( x 3 m , t n ) is denoted by U m n .
The derivatives of even order in (37) are discretized using central differences, while the derivatives of the first order are discretized using forward differences. The resulting scheme is implicit. Applied at each inner point of the spatial mesh, it forms a tridiagonal matrix problem, which is effectively solved using the Thomas algorithm.
Using the standard method of Taylor expansions, it is easy to check that, with sufficient smoothness conditions, the order of accuracy of this scheme is O ( h 2 + τ ) . Note that a central-difference discretization of the first-order derivatives would result in higher accuracy in time; however, the stability of such a scheme would be hardly attainable in practical computations, where the parabolic terms of (37) are dominant.
Using the classical stability analysis, it is easy to show that, for the practical problems that we are interested in, the numerical stability Courant–Friedrichs–Lewy (CFL) condition τ c h 2 holds, with some constant c that depends on the parameters of the model.
We seek the scaled average velocity Q in m/s inside a small elastic tube or capillary with the Young modulus E = 10 6 Pa, the Poisson ratio ν ^ = 0.35 , the length of the tube 10 3 m, the radius 5 · 10 5 m, the thickness ε = 2 · 10 5 m and the dynamic viscosity of blood ν = 4 · 10 3 Pa·s. The following set of constants C i is obtained:
C 1 = 9.2287 · 10 38 , C 2 = 2.0734 · 10 42 , C 3 = 4.1667 · 10 2 , C 4 = 1.7904 · 10 1 , C 5 = 7.1615 · 10 1 , C 6 = 6.2500 · 10 2 , C 7 = 4.000 · 10 4 .
To test the accuracy of the constructed solver, let us define the error in the maximum norm by e, and the experimental convergence rates in space and time by e h and e τ :
e ( h , τ ) = max m U m N U ( x 3 m , T ) , e h ( h ) = log 2 e ( 2 h , τ ) e ( h , τ ) , e τ ( τ ) = log 2 e ( h , 2 τ ) e ( h , τ ) .
Here, U ( x 3 m , T ) is a benchmark solution, which is computed using very small step sizes.
The results of experimental convergence tests are presented in Table 1 below. Here, the functions g i n ( t ) = 0.005 sin ( 2 t ) , g o u t = 0.001 e t were used, and the problem was solved for T = 1 . We can see from the results of Table 1 that the experimental convergence rate agrees well with the theoretical estimate.

6. Numerical Comparison of the Approximate 1D Model and the Full-Dimensional Problem

In the present section, we compare the numerical solutions of two problems: the full-dimensional Navier–Stokes equations with the fluid–structure interaction (FSI) conditions for the tubes with an elastic wall, and the EWPE approximation.
  • Solution in a single vessel
We consider two types of inflow/outflow boundary conditions: first, g i n ( t ) = g o u t ( t ) ; second, with g o u t ( t ) = g i n ( t + δ ) ; and compare the solutions in the middle point of the tube (Figure 2 and Figure 3).
For δ = 0 and δ = 0.1 , the difference between the solutions is less than 10 per cent. When δ increases ( δ = 0.2 , δ = 0.25 ), the difference is more significant. This can be explained by important variations of the cross-section (more than two times). In this case, the assumption of small displacements of the wall is no more realistic, and so the model must take into consideration the difference between Eulerian and Lagrangian approaches.
  • Solution in a Y-shaped network of vessels
Next, we compare the solutions acquired from (37) and the full Navier–Stokes in a Y-shaped network (see Figure 4), where the inflow is on the left and two outflows are on the right. The corresponding boundary conditions are g i n = 0.005 sin ( 2 t ) and g o u t = 0.0025 sin ( 2 t + 0.02 ) for each of the two outlets to the right, while the Young modulus is E = 10 6 Pa for the tube-shaped blood vessel wall and E = 10 10 Pa for the sphere-shaped connection of three tubes, the Poisson ratio ν ^ = 0.35 , the length of each of the three tubes 10 2 m, the radius 10 3 m, the wall thickness ε = 10 4 m, and the dynamic viscosity of blood ν = 4 · 10 3 Pa·s. The large Young’s modulus of the wall in the bifurcation zone simulates the no-slip boundary condition there. It corresponds to the experimental observation that, near the bifurcations, the vessels become much stiffer than the “tube” walls.
The EWPE approximation uses the following junction conditions at the bifurcation node: (1) the flux 4 π ε 1 4 Q coming from the left tube is equal to the sum of the fluxes going out into the left tubes (this condition corresponds to a stiffer bifurcation zone of vessels); (2) the first derivatives Q x 3 = 0 for all tubes at the bifurcation node; and (3) the second derivatives 2 Q x 3 2 = 0 at the bifurcation node for the tubes on the right.
Figure 5 and Figure 6 display the graphs of average velocity in the midpoints of the left and right vessels, respectively, up to T = 10.

7. Conclusions

A 1D elastic wall quasi-Poiseuille approximation for the fluid–structure interaction problem in a cylindrical tube is introduced and tested numerically by comparison with the full-dimensional Navier–Stokes–elastic shell FSI problem. One-tube and tube network geometries are considered. The inflows and outflows with and without flux conservation are tested. The results show a good approximation for the fluid–structure interaction model in the case of small displacements of the wall and modest Reynolds numbers.

Author Contributions

Conceptualization: G.P.; Formal analysis: K.K., N.K., G.P., K.P., V.Š.; Software, numerical computations N.K., V.Š. All authors have read and agreed to the published version of the manuscript.

Funding

The work has received funding from the European Social Fund (project No 09.3.3-LMT-K-712-17-0003) under a grant agreement with the Research Council of Lithuania (LMTLT).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Desjardins, B.; Esteban, M.J.; Grandmont, C.; le Talec, P. Weak solutions for a fluid-structure interaction model. Rev. Mat. Comput. 2001, 14, 523–538. [Google Scholar] [CrossRef] [Green Version]
  2. Gilbert, R.P.; Mikelić, A. Homogenizing the acoustic properties of the seabed: Part I. Nonlin. Anal. Theory Methods Appl. 2000, 40, 185–212. [Google Scholar] [CrossRef] [Green Version]
  3. Grandmont, C.; Maday, Y. Existence for an unsteady fluid-structure interaction problem. Math. Model. Numer. Anal. 2000, 34, 609–636. [Google Scholar] [CrossRef]
  4. Čanić, S.; Mikelić, A. Effective equations modeling the flow of a viscous incompressible fluid through a long elastic tube arising in the study of blood flowthrough small arteries. SIAM J. Appl. Dyn. Syst. 2003, 2, 431–463. [Google Scholar] [CrossRef] [Green Version]
  5. Čanić, S.; Tambača, J.; Guidoboni, G.; Mikelić, A.; Hartley, C.J.; Rosenstrauch, D. Modeling viscoelastic behavior of arterial walls and their interaction with pulsate blood flow. SIAM J. Appl. Math. 2006, 67, 164–193. [Google Scholar] [CrossRef] [Green Version]
  6. Mikelić, A.; Guidoboni, G.; Čanić, S. Fluid-structure interaction in a pre-stressed tube with thick elastic walls. I. The stationary Stokes problem. Netw. Heterog. Media 2007, 2, 397–423. [Google Scholar] [CrossRef]
  7. Panasenko, G.P.; Stavre, R. Asymptotic analysis of a periodic flow in a thin channel with visco-elastic wall. J. Math. Pures Appl. 2006, 85, 558–579. [Google Scholar] [CrossRef] [Green Version]
  8. Panasenko, G.P.; Sirakov, Y.; Stavre, R. Asymptotic and numerical modeling of a flow in a thin channel with visco-elastic wall. Int. J. Multiscale Comput. Eng. 2007, 5, 473–482. [Google Scholar]
  9. Panasenko, G.P.; Stavre, R. Asymptotic analysis of a non-periodic flow with variable viscosity in a thin elastic channel. Netw. Heterog. Media 2010, 5, 783–812. [Google Scholar] [CrossRef]
  10. Malakhova-Ziablova, I.; Panasenko, G.; Stavre, R. Asymptotic analysis of a thin rigid stratified elastic plate—Viscous fluid interaction problem. Appl. Anal. 2016, 95, 1467–1506. [Google Scholar] [CrossRef]
  11. Panasenko, G.P.; Stavre, R. Viscous fluid-thin elastic plate interaction: Asymptotic analysis with respect to the rigidity and density of the plate. Appl. Math. Optimiz. 2020, 81, 141–194. [Google Scholar] [CrossRef]
  12. Panasenko, G.P.; Stavre, R. Viscous fluid-thin cylindrical elastic body interaction: Asymptotic analysis on contrasting properties. Appl. Anal. 2019, 98, 162–216. [Google Scholar] [CrossRef]
  13. Panasenko, G.P.; Stavre, R. Three dimensional asymptotic analysis of an axisymmetric flow in a thin tube with thin stiff elastic wall. J. Math. Fluid Mech. 2020, 22, 35. [Google Scholar] [CrossRef]
  14. Formaggia, L.; Lamponi, D.; Quarteroni, A. One-dimensional models for blood flow in arteries. J. Eng. Math. 2003, 47, 251–276. [Google Scholar] [CrossRef]
  15. Gamilov, T.M. Mathematical Modeling of Blood Flow under Mechanical Action at Vessels. Ph.D. Thesis, Moscow Physics and Technology Institute, Moscow, Russia, 2017. [Google Scholar]
  16. Koshelev, V.B.; Mukhin, S.I.; Sosnin, N.V.; Favorsky, A.P. Mathematical Models for Quasi-One-Dimensional Hemodynamics; MAKS Press: Moscow, Russia, 2010. [Google Scholar]
  17. Simakov, S.S.; Kholodov, A.S.; Evdokimova, A.V. Methods for computing global blood flow in the human’s body by heterogeneous numerical models. In Medicine in the Mirror of Informatics; Nauka: Moscow, Russia, 2008; pp. 124–170. [Google Scholar]
  18. Long, G.; Liu, Y.; Xu, W.; Zhou, P.; Zhou, J.; Xu, G.; Xiao, B. Analysis of crack problems in multilayered elastic medium by a consequtive stiffness method. Mathematics 2022, 10, 4403. [Google Scholar] [CrossRef]
  19. Boussinesq, J. Essaii sur la Théorie des Eaux Courrantes. In Mémoires Présentés par Divers Savants de l’Académie des Sciences de l’Institut de France; Impr. Nationale: Paris, France, 1877. [Google Scholar]
  20. Poiseuille, J.L.M. Recherches sur les Causes du Mouvement du Sang dans les Vaisseaux Capillaires; Impr. Royale: Paris, France, 1839. [Google Scholar]
  21. Galdi, G.P.; Pileckas, K.; Silvestre, A.L. On the unsteady Poiseuille flow in a pipe. Z. Angew. Math. Phys. 2007, 58, 994–1007. [Google Scholar] [CrossRef]
  22. Womersley, J.R. Method of calculation of velocity, rate of the flow and viscous drag in arteries when the pressure gradient is known. J. Physiol. 1955, 127, 553–563. [Google Scholar] [CrossRef] [PubMed]
  23. Berntsson, F.; Ghosh, A.; Kozlov, V.A.; Nazarov, S.A. A one-dimensional model of blood flow through a curvilinear artery. Appl. Math. Model. 2018, 63, 633–643. [Google Scholar] [CrossRef]
  24. Panasenko, G.P.; Stavre, R. Asymptotic analysis of a viscous fluid-thin plate interaction: Periodic flow. Math. Mod. Meth. Appl. Sci. 2014, 24, 1781–1822. [Google Scholar] [CrossRef]
  25. Evans, L.C. Partial Differential Equations. In Graduate Studies in Mathematics, 2nd ed.; American Mathematical Society: Providence, Rhode Island, 2010; Volume 19. [Google Scholar]
  26. Temam, R. Navier-Stokes Equations; North-Holland Publishing Company: Amsterdam, The Netherlands; New York, NY, USA, 1977. [Google Scholar]
Figure 1. A thin tube with radius ε 1 , covered by an elastic wall of thickness ε .
Figure 1. A thin tube with radius ε 1 , covered by an elastic wall of thickness ε .
Mathematics 11 02106 g001
Figure 2. (a) presents the average velocity computed by EWPE (the fourth order PDE (37)) over time, (b) depicts the average solution (velocity) of the full Navier–Stokes model. In both cases, g i n = 0.005 sin ( 2 t ) , g o u t = 0.005 sin ( 2 t + 0.02 ) and the solution was taken from the middle cross-section of the tube.
Figure 2. (a) presents the average velocity computed by EWPE (the fourth order PDE (37)) over time, (b) depicts the average solution (velocity) of the full Navier–Stokes model. In both cases, g i n = 0.005 sin ( 2 t ) , g o u t = 0.005 sin ( 2 t + 0.02 ) and the solution was taken from the middle cross-section of the tube.
Mathematics 11 02106 g002
Figure 3. (a) presents the average velocity according to the EWPE approximation over time, (b) depicts the average velocity of the full Navier–Stokes model. In both cases, g i n = 0.005 sin ( 2 t ) , g o u t = 0.005 sin ( 2 t + 0.1 ) and the solution was taken from the middle cross-section of the tube. Note that here we deliberately simulate a case of large displacement of the wall, in order to explore the limitations of the ansatz.
Figure 3. (a) presents the average velocity according to the EWPE approximation over time, (b) depicts the average velocity of the full Navier–Stokes model. In both cases, g i n = 0.005 sin ( 2 t ) , g o u t = 0.005 sin ( 2 t + 0.1 ) and the solution was taken from the middle cross-section of the tube. Note that here we deliberately simulate a case of large displacement of the wall, in order to explore the limitations of the ansatz.
Mathematics 11 02106 g003
Figure 4. Geometry of the Y-shaped bifurcation model.
Figure 4. Geometry of the Y-shaped bifurcation model.
Mathematics 11 02106 g004
Figure 5. (a) presents the average velocity given by the EWPE approximation (37) over time, (b) depicts the average velocity given by the full Navier–Stokes model. In both cases, g i n = 0.005 sin ( 2 t ) , g o u t = 0.0025 sin ( 2 t + 0.02 ) , and the solution was taken from the middle cross-section of the left tube.
Figure 5. (a) presents the average velocity given by the EWPE approximation (37) over time, (b) depicts the average velocity given by the full Navier–Stokes model. In both cases, g i n = 0.005 sin ( 2 t ) , g o u t = 0.0025 sin ( 2 t + 0.02 ) , and the solution was taken from the middle cross-section of the left tube.
Mathematics 11 02106 g005
Figure 6. (a) presents the average velocity given by the EWPE approximation (37) over time, (b) depicts the average velocity given by the full Navier–Stokes model. In both cases, g i n = 0.005 sin ( 2 t ) , g o u t = 0.0025 sin ( 2 t + 0.02 ) , and the solution was taken from the middle cross-section of either of the right tubes.
Figure 6. (a) presents the average velocity given by the EWPE approximation (37) over time, (b) depicts the average velocity given by the full Navier–Stokes model. In both cases, g i n = 0.005 sin ( 2 t ) , g o u t = 0.0025 sin ( 2 t + 0.02 ) , and the solution was taken from the middle cross-section of either of the right tubes.
Mathematics 11 02106 g006
Table 1. Errors and experimental convergence rates for sequences of step sizes in space (left) and time (right). The results confirm experimentally that the order of accuracy of the numerical scheme is O ( h 2 + τ ) .
Table 1. Errors and experimental convergence rates for sequences of step sizes in space (left) and time (right). The results confirm experimentally that the order of accuracy of the numerical scheme is O ( h 2 + τ ) .
h e ( h , τ ) e h ( h ) τ e ( h , τ ) e τ ( τ )
1 · 10 4 3.1982 · 10 6 2.0019 5 · 10 5 1.1683 · 10 10 1.0203
5 · 10 5 7.9850 · 10 7 1.9957 2.5 · 10 5 5.7601 · 10 11 1.0409
2.5 · 10 5 2.0022 · 10 7 1.9984 1.25 · 10 5 2.7996 · 10 11 1.0900
1.25 · 10 5 5.0111 · 10 8 2.0105 6.25 · 10 6 1.3152 · 10 11 0.9782
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Kaulakytė, K.; Kozulinas, N.; Panasenko, G.; Pileckas, K.; Šumskas, V. Poiseuille-Type Approximations for Axisymmetric Flow in a Thin Tube with Thin Stiff Elastic Wall. Mathematics 2023, 11, 2106. https://doi.org/10.3390/math11092106

AMA Style

Kaulakytė K, Kozulinas N, Panasenko G, Pileckas K, Šumskas V. Poiseuille-Type Approximations for Axisymmetric Flow in a Thin Tube with Thin Stiff Elastic Wall. Mathematics. 2023; 11(9):2106. https://doi.org/10.3390/math11092106

Chicago/Turabian Style

Kaulakytė, Kristina, Nikolajus Kozulinas, Grigory Panasenko, Konstantinas Pileckas, and Vytenis Šumskas. 2023. "Poiseuille-Type Approximations for Axisymmetric Flow in a Thin Tube with Thin Stiff Elastic Wall" Mathematics 11, no. 9: 2106. https://doi.org/10.3390/math11092106

APA Style

Kaulakytė, K., Kozulinas, N., Panasenko, G., Pileckas, K., & Šumskas, V. (2023). Poiseuille-Type Approximations for Axisymmetric Flow in a Thin Tube with Thin Stiff Elastic Wall. Mathematics, 11(9), 2106. https://doi.org/10.3390/math11092106

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop