Numerical method for solving the inverse problem of non-stationary flow of viscoelastic fluid in the pipe

It is known that viscoelastic fluids possess the property of elastic recovery of its shape, characteristic of solids and characteristics of viscous flow typical for fluids. Such properties are shown by mixtures of polymers, dough, oil and petroleum products with a high content of resins, bitumen, etc. For viscoelastic fluids, two different rheological models were proposed that correspond to two different approaches to determining the joint action of elastic forces and viscosity of fluids [1–3]. Usually, mechanical models of viscoelastic fluids are represented by a combination of elastic and viscous elements (Hooke and Newton models). The rheological model proposed by Maxwell is represented as a sequential connection of elastic and viscous elements. According to Maxwell’s model, the strain rate of viscoelastic fluids consists of the elastic strain rate and the viscous strain rate. A rheological model of viscoelastic fluids, proposed by Kelvin and Voigt, is represented as a parallel connection of elastic and viscous elements. In this case, the total tangent stress is represented as a simple sum of the stress corresponding to the elastic deformation and the stress caused by the viscous resistance. In the Kelvin–Voigt rheological model, the ratio between total stress and strain is written as an ordinary differential equation with respect to strain

∂ε

■ =^dt ⁺ E^£,

where σ is the tangent stress, ε is the strain that occurs under the influence of stress, E is the modulus of elasticity, µ is the coefficient of dynamic viscosity.

Currently, viscoelastic fluids, including artificially created ones, are widely used in the aviation, food, oil, chemical industries and many other branches of mechanical engineering. In many technological processes in these industries, the flow of viscoelastic fluids is one of the most important elements. Therefore, modelling the flow of viscoelastic fluids in different media is of great practical importance. General principles of construction of mathematical models of viscoelastic fluids, issues of numerical simulation of the flow of viscoelastic fluids in various media are studied in [4–7].

1. Problem Statement

Let us consider a non-stationary axisymmetric flow of an incompressible viscoelastic fluid in a horizontally arranged cylindrical tube with a constant cross-section. The mathematical model of this flow is presented as follows [8]:

p duujrtt = 1 d ^(r^) ^ ^tti, ₀ < r < R, о < t,         (2)

∂t r ∂r                l

where u(r, t) is the rate of flow of a viscoelastic fluid directed parallel to the axis of the pipe, l is the pipe length, AP(t) is the pressure drop along the length of the pipe, p is the fluid density, R is the radius of the tube.

Let us assume that a viscoelastic fluid satisfies the Kelvin–Voigt rheological model (1) and there is no deformation in the fluid at the initial time. Then, for the given ratio

8e(r, t)    6u(r,t)

∂t ∂r the Kelvin–Voigt rheological model is written as

t

*,1) = / E J      dT.

Substituting the last stress representation in equation (2), we obtain the following integrodifferential equation with respect to the flow velocity of a viscoelastic fluid

t f ди(г,т) d \    AP(t)

∂r            lρ

0 <  r < R, 0 < t <  T.

Suppose that equation (3) is endowed with the initial condition

u | t =0 = о,

natural boundary value condition for r = 0

MM = n

dr

and the adhesion condition on the pipe wall

u ( R, t ) = 0 .

Obviously, setting the law of change of pressure drop AP(t) in time, solving direct problem (3) – (6) to find the velocity distribution for viscoelastic fluid flow over the cross section of the pipe and the volumetric flow rate through the pipe. Let us assume that in problem (3) - (6), the function AP(t) is unknown along with the function u(r, t) , and we need to determine AP(t) by the specified volume flow of the fluid through the pipe

R j 2nru(r,t)dr=q(i),                                  (7)

where q(t) is the volume flow of fluid through the pipe.

Therefore, the problem is to determine the functions u(r,t) and AP(t) that satisfy equation (3) and conditions (4) – (7). This problem belongs to the class of inverse problems related to the recovery of the right parts of integro-differential equations [8–14].

• 2.    Method to Solve Problem

As a result of the replacement

t

w

⁽ r, ^{t) =} j

u(r, т)dr,

equation (3) is written as d2w(r,t)    ц d / d2w(r, t)\   Ed / dw(r, t)\   AP(t)

∂t2      rρ ∂r ∂r∂t      rρ ∂r ∂rlρ

0 < r < R, 0 < t <  T.

For equation (8), we have the following initial and boundary value conditions:

w(r,0) = 0,        ,    = 0,

^r-t—fl. w(R,t)=0.(10)

In this case, additional condition (7) is converted to the form

R

/            = Q(,)'

t where Q(t) = J д(т)dT.

Let us construct a difference analog of problem (8) – (11). To this end, we introduce the uniform difference grid ш = {(tj, ri) : ri = iAr, tj = jAt, i = 0,1, 2,..., n, j = 0,1, 2,..., m} in the rectangular area {0 < r < R, 0 < t < T} with the increment Ar = R/n of the variable r and the increment At = T/m of the time t. In the inner nodes of the grid ш, we associate equation (8) with an implicit difference scheme j+1 j     j-1                         j+1     j+1            j+1     j+1

w i   ^- 2w i + w = ц          ^w i +i ^- w i   -      w i ^- w i - 1

At ²            r i pArAt    i ^+1/ 2      Ar          i ^1/2      Ar

U        wj+1 - wj       wj - w— riPAA     —at" - ri-1/2   a

E Г    w j +1 - w j w j - wpH   AP ^{j +1}

+r^Ar ri+1/2ri-1/2   Ar -  + — i = 1, 2, 3,...,n - 1, j = 1, 2, 3,...,m - 1, where wj « w(ri, tj), APj « AP(tj), ri±1/2 = ri + Ar/2.

Approximating conditions (9) – (10), we have w0 = 0,

j +1 w ₁

-

w i ¹ -'

At j+1 w0

w0         ---

— = 0, i = 0, n,

Ar

= 0, w n ⁺¹ = 0.

And the discrete analog of additional condition (11) is written as

n

^ 2nY i r i w ^{j +1} = Q j ⁺¹ , i =1

where Q ^j+1 ~ Q(t j ₊₁ ) , Y i are coefficients of the quadrature formula. The resulting system of difference equations is converted to

1                             ________________ aiwp1 - Ciwj  + biwp1 = fi -   APj+ , i = 1, n - 1, lρ wj+1 = w0+1, wn+1 = 0,

V 2nYiriwj+1 = Qj+1, i=1

j = 1, 2, ...m - 1, w0 = 0, wi1 = w0, i = 0,n, where ai =

Uri-1/2    , ripAr2At

Er i - 1 / 2 ь r i pAr ² , i

Ur i +1 / 2     , Er i +1 / 2             , , , 1

---л—i— +-- л—т, C i — a i + b i +— л—о , r i pAr ² At   r i pAr ²                 At²

f i =

µ ripArAt

^{w j} +1 - w j         w j - ^{w j}- 11    ^{2 w j} - w j ¹

^{r ,+1/2} ^r     r ' ^-1/2    Ar           At 2

Difference problem (12) – (16) is a system of linear algebraic equations in which, as unknowns, we use the approximate values of the desired functions w(r, t) and AP(t) in nodes of the difference grid, i.e. w j ⁺¹ , AP ^j+1 , i = 0, ...,n , j = 1, ...,m .

i

In order to divide difference problem (12) – (14) into mutually independent subproblems, each of which can be solved independently, the solution to this system for each fixed value j , j = 1, 2,..., m - 1 , is represented as [8,9]

w j ⁺¹ = 0 j ⁺¹ + AP j ⁺¹ ф j ⁺¹ , i = 0,1, 2,...,n,                    (17)

where 0 ^j+1 , ф ^j+1 are unknown variables yet.

Substituting the expression for w j ⁺¹ in each equation of system (12) - (14), we get

[ ai Л

- C i 9 j ⁺¹ + b 9 ; - f i ] + AP ^j+1 ^ф^

^- c i ^ j' ⁺¹ + b i ^ j' +i + ipj — 0,

9 j ⁺¹ + AP j ' Ф 1/ — 9 j ₀ ⁺¹ + AP ^{j +}V ₀ ⁺¹ ,

9 n_n ⁺¹ + AP ^{j +1} ф ?_т + ¹ — 0.

From the last relations we get the following difference problems for determining the auxiliary variables 9 ^j+1 , ф 0 ⁺¹

ai9j+ - a-' + bi9j+ - fi — 0, i — 1, 2,...,n - 1,

9j+1 — 90+1,

- ' — 0.

ai ФД - Ci^i+1 + ьф <    0, i — 1, 2, ...,n - 1, lP

Ф1+1 — Ф0+1,

■     — 0.

j — 1, 2, 3,...,m - 1.

For each fixed value j — 1, 2, ...,m - 1 , resulting difference problems (18) - (20) and (21) – (23) are given by a system of linear algebraic equations with a tridiagonal matrix and solutions to these systems can be found independently of AP ^j+1 by the Thomas method [9].

And substituting representation (17) in (15), we have nn

^ 2nY i r i 9 j ⁺¹ + AP j ⁺¹ ^ 2пY i r i ф j' ⁺¹ — Q ^j+1 .

i =1                           i =1

From here, we can determine the approximate value of the desired function AP(t) for

t     t j +1 , i- ^e -

_AP j +1 _—

Q j ⁺¹ - EL2nY i r i ■■ ' EL 2nY i r ■

Thus, for each fixed value j — 1, 2, 3,..., m - 1 , the computational algorithm for solving difference problem (12) - (16) by definition of w ^{j +1} , i — 0, n and AP j ⁺¹ is as follows:

• 1.    Solutions to two second–order linear difference problems (18) – (20) and (21) – (23) with respect to the auxiliary variables 9 ^j+1 , ф 0 ⁺¹ , i — 0,n , are determined;

• 2.    Formula (24) defines AP ^{j +1} ;

• 3.    The values of variables w ^{j +1} , i — 0,n , are calculated using formula (17).

It should be noted that from the ratio

u(r, t) —

dw(r, t) ∂t

using numerical differentiation procedures, it is possible to find the velocity distribution for the viscoelastic fluid flow over the pipe section in each time layer.

Table

Results of the numerical experiment

t , s	AP t ,	AP3 ,	AT3 ,	MPa
	MPa	MPa	δ =0,02	δ =0,05
200	2,175	2,175	2,189	2,210
400	6,209	6,209	6,211	6,214
600	5,569	5,569	5,602	5,651
800	2,005	2,005	2,040	2,092
1000	5,264	5,264	5,342	5,458
1200	6,433	6,433	6,458	6,496
1400	2,315	2,315	2,340	2,378
1600	4,172	4,172	4,228	4,311
1800	6,925	6,925	6,987	7,079
2000	3,045	3,045	3,078	3,127
2200	3,144	3,144	3,156	3,173
2400	6,952	6,952	7,074	7,257
2600	4,054	4,054	4,073	4,100
2800	2,376	2,376	2,388	2,407
3000	6,507	6,507	6,507	6,508
3200	5,149	5,149	5,160	5,176
3400	2,016	2,016	2,048	2,095
3600	5,676	5,676	5,734	5,820
3800	6,120	6,120	6,123	6,129
4000	2,134	2,134	2,141	2,142

3.    Results of Numerical Calculations

To find out the effectiveness of the proposed computational algorithm, numerical experiments were performed for model problems. Numerical experiments were carried out according to the following scheme:

• -    for the given function AP(t) , 0 <  t <  T , find a solution to problem (8) - (10), i.e. the function w(r, t) , 0 <  r <  R, 0 <  t <  T ;

R

• -    consider the found dependency Q(t) = J 2nrw(r,t)dr as accurate data for solving 0

the inverse A P ( t ) recovery problem.

The first series of calculations was performed using these undisturbed data. The second series of calculations was performed by applying a function Q(t) to model the error of experimental data

Q(t) = Q(t) + 5n(t)Q(t), where n(t) is a random process simulated using the random number generator; 5 is the error level. For perturbation of input data, we consider the error level to be 5 = 0, 02; 0, 05.

Numerical calculations were performed using a spacetime difference grid with increments Ar = 0, 03m , At = 0, 005; 1; 10s . The results of the numerical experiment performed for the case ^ = 0,06Pa • s ; p = 900kq/m ³ ; R = 0,6m ; AP(t) = 4,5 — 2, 5sin10tMPa ; E = 200Pa ; At = 10s ; L = 10000m using undisturbed and disturbed input data are presented in Table; where t is time, AP t are the exact values of the function

AP(t) , AP are the calculated values of AP(t) for undisturbed data, AP are the calculated values of AP(t) for disturbed data.

Results of the numerical experiment show that with undisturbed input data, the desired function AP(t) is restored exactly for all calculated grids in time (Table, column 3). And when using perturbed input data, where the error has a fluctuating character, the desired function AP(t) is recovered with an error. At the same time, the use of fairly small time (At <  0.005s) steps gives the opposite effect compared to the numerical solution of direct boundary value problems, i.e. with a decrease in the time step, the error in restoring the function AP(t) increases. However, for the case of perturbed input data, it is not possible to theoretically determine the range of the time step at which the solution to the inverse problem is stable. Therefore, for perturbed input data, the time step was determined by numerical experimentation. Thus, when using At = 10s in calculations, the maximum relative error of restoring the values of the desired function AP(t) did not exceed 1, 76% at the error level 5 = 0, 02 and 4,42% at 5 = 0, 05 .

Analysis of the results of the numerical experiment shows that the proposed computational algorithm provides stability of the solution to errors in the input data.

Conclusion

The problem of determining the pressure drop in a non-stationary flow of a viscoelastic fluid in a cylindrical pipe is considered based on information about the change in time of the volume flow of the fluid through the pipe. To solve this problem, we propose a method based on converting an integro-differential equation into a third-order differential equation, discretizing the resulting problem, and using a special representation to separate the desired variables. The proposed method allows us to determine the pressure drop along the length of the pipe in each time layer.

Acknowledgements. The authors received financial support from Taif University Researches Supporting Project number TURSP-2020/031, Taif University, Taif, Saudi Arabia. All authors contributed equally to this work.