CC BY
15
MSC 65M32, 76A10 DOI: 10.14529/mmp220408
NUMERICAL METHOD FOR SOLVING THE INVERSE PROBLEM OF NON-STATIONARY FLOW OF VISCOELASTIC FLUID IN THE PIPE
A.R. Aliev1'2, Kh.M. Gamzaev2, A.A. Darwish3, T.A. Nofal45
Azerbaijan State Oil and Industry University, Baku, Azerbaijan 2Institute of Mathematics and Mechanics of ANAS, Baku, Azerbaijan 3Helwan University, Cairo, Egypt 4Taif University, Taif, Saudi Arabia 5El-Minia University, Minia, Egypt
E-mail: [email protected], [email protected], [email protected], [email protected]
The process of unsteady flow of incompressible viscoelastic fluid in a cylindrical tube of constant cross-section is considered. To describe the rheological properties of a viscoelastic fluid, the Kelvin-Voigt model is used and the mathematical model of this process is presented as an integro-differential partial differential equation. Within the framework of this model, the problem is to determine the pressure drop along the length of the pipe, which ensures the passage of a given flow rate of viscoelastic fluid through the pipe. This problem belongs to the class of inverse problems related to the recovery of the right parts of integro-differential equations. By replacing variables, the integro-differential equation is transformed into a third-order partial differential equation. First, a discrete analog of the problem is constructed using finite-difference approximations. To solve the resulting difference problem, we propose a special representation that allows splitting the problems into two mutually independent second-order difference problems. As a result, an explicit formula is obtained for determining the approximate value of the pressure drop along the length of the pipeline for each discrete value of the time variable. Based on the proposed computational algorithm, numerical experiments were performed for model problems.
Keywords: viscoelastic fluid; Kelvin-Voigt model; integro-differential equation; pressure drop along the length of the pipe; inverse problem.
Introduction
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
<J = ^ + Ee, (1)
where a is the tangent stress, e 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-
du(r, t) д t
1 d . . .. AP(t) r dr l
0 < r < R, 0 < t < T,
(2)
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
de(r,t) du(r,t) dt dr
the Kelvin-Voigt rheological model is written as
du(r, t) dr
¿t.
Substituting the last stress representation in equation (2), we obtain the following integro-differential equation with respect to the flow velocity of a viscoelastic fluid
du(r, t) dt
ц d rp dr
du(r, t) dr
E д rp dr
du(r, t) dr
¿t +
0 <r <R, 0 <t < T. Suppose that equation (3) is endowed with the initial condition
u|t=o = 0,
natural boundary value condition for r = 0
du(0, t)
dr
0
and the adhesion condition on the pipe wall
u(R,t) = 0.
mt) lp
(3)
(4)
(5)
(6)
t
Вестник !Ю"УрГ"У. Серия «Математическое моделирование и программирование» (Вестник ЮУрГУ ММП). 2022. Т. 15, № 4. С. 90-98
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(t), (7)
0
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
w(r,t) = J u(r,T)dr, 0
equation (3) is written as
d2w(r,t) » d dt2 rpdr
d2w(r, t)
E d
+ —7T 1 r"
rp or
8w(r,t)\ AP (t)
dr
+
Ip
drdt
0 < r < R, 0 <t < T. For equation (8), we have the following initial and boundary value conditions:
(8)
w(r, 0) = 0
dw(0,t) _ dr
dw(r, 0) ' dt
0, w(R,t)
0,
0.
In this case, additional condition (7) is converted to the form
R
2nrw(r,t)dr = Q(t),
(9) (10)
(11)
where Q(t) = f q(r)dr. 0
Let us construct a difference analog of problem (8) - (11). To this end, we introduce the uniform difference grid
u
{(tj,n) : r = 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 u, we associate equation (8) with an implicit difference scheme
wj+1 — 2wj + wj 1
»
At2
ripArAt
ri+l/2~
j+1 wj+1 — w'
j+1
w
j+1
—w
j+1
Ar
— ri-
1/2-
i-1
Ar
t
ß
ri pArAt
ri+l/2-
wi+i — wi
, j
w — w
3
Ar
— Vi-
1/2-
i-1
Ar
E
ripAr
„3+1
Гг+1/2"
wi+i — w:
3+i
3+i
w — w,
3+i
Ar
—ri-
i-1/2"
i-i
Ar
+
+
AP 3+l
lp
i = 1, 2, 3, ...,n — 1, j = 1, 2, 3, ...,m — 1,
where wj « w(r,,tj), APj « AP(tj), ri±1/2 = r + Ar/2. Approximating conditions (9) - (10), we have
n wj — w®
w°i = 0, -L_—L = o, i = 0,n,
At
3+i 3+i
wl — W0
Ar
0, wj+l = 0.
And the discrete analog of additional condition (11) is written as
n
Y,2nllrlwW^1 = Qj+1,
i=l
where Qj+1 ~ Q(tj+1), y, are coefficients of the quadrature formula. The resulting system of difference equations is converted to
j+i
ai wi-l — ciwi
3+i
+ bi w:
3+l i+l
fi — —APj+1, г = l,n — 1, lp
w
j+l
w
w
3+l
3+l о ,
0,
y^ 2-KYiViw:
3+l
Q
3+l
i=l
j = 1, 2, ...m — 1,
w
0 = 0, wl
w0, i
0, n,
where a,
ßTi-i/2 EVj-1/2
л ^ i
ripAr2At ripAr2
ßVi+l/2 , Eri+l/2
-'- ~г -— C'
ripAr2At ripAr2 i
ai + bi +
At2
fi
ß
ripArAt
Гг+1/2"
wi+l — wi
Ar
— ri-
1/2-
33
wi — w
i- l
Ar
2wi3 — wi
3-l
At2
(12)
(13)
(14)
(15)
(16)
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. wj+1, APj+1, i = 0, ...,n, j = 1, ...,m.
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
3+l = 3+l + AP3+lф3+l, i = 0,1, 2,...,
n,
(17)
where 9j+1, +1 are unknown variables yet.
Вестник ЮЮУрГУ. Серия «Математическое моделирование и программирование» (Вестник ЮУрГУ ММП). 2022. Т. 15, № 4. С. 90-98
l
1
Substituting the expression for wj+1 in each equation of system (12) - (14), we get
h — ci ej+1 + — fi] + APj+1
a^l - Cl<pr + bittlï + Yp
j+1
j+1
0,
j1 + ap j+V1+1 = eo+1 + APj+Vo+\ en+1 + ap j+10n+1 = 0.
From the last relations we get the following difference problems for determining the auxiliary variables ej+1, 0j'+1
ajl — c% ej+1 + bi ej+11 — fi = 0, i = 1, 2,...,n — 1,
ej
j +1
e0+1,
en+1 = 0.
aicpl^l - Ci4>l+1 + bi^i+l + = i
1, 2, ...,n — 1,
j1,
= 0.
j = 1, 2, 3, ...,m — 1.
(18)
(19)
(20)
(21)
(22) (23)
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 APj+1 by the Thomas method [9].
And substituting representation (17) in (15), we have
£ 2nllrlej+1 + APj+1 £ 2nllrl^l+1 = Qj+1.
i=1
i=1
From here, we can determine the approximate value of the desired function AP(t) for t = tj+1, i.e.
APi+i = EIL^vr^r (24)
Thus, for each fixed value j = 1, 2, 3,..., m — 1, the computational algorithm for solving difference problem (12) - (16) by definition of wj+1 ,i = 0,n and APJ+1 is as follows:
1. Solutions to two second-order linear difference problems (18) - (20) and (21) - (23) with respect to the auxiliary variables 0{+1, (¡>{+l, i = 0,n, are determined;
2. Formula (24) defines APj+1;
3. The values of variables wj+1, i = 0,n, are calculated using formula (17). It should be noted that from the ratio
u(r, t)
dw(r, t) dt
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 А Р\ МРа АР, МРа АР, МРа
¿=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) = f 2nrw(r,t)dr as accurate data for solving
0
the inverse AP(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) + SV(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/m3; 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, AP1 are the exact values of the function
Вестник !Ю"УрГ"У. Серия «Математическое моделирование и программирование» (Вестник ЮУрГУ ММП). 2022. Т. 15, № 4. С. 90-98
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 8 = 0, 02 and 4, 42% at 8 = 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.
References
1. Wilkinson W.L. Non-Newtonian Fluids: Fluid Mechanics. Mixing and Heat Transfer. New York, Pergamon Press, 1960.
2. Astarita G., Marrucci G. Principles of Non-Newtonian Fluid Mechanics. London, New York, McGraw-Hill, 1974.
3. Joseph D.D. Fluid Dynamics of Viscoelastic Liquids. New York, Springer, 1990.
4. Huilgol R. R., Phan-Thien N. Fluid Mechanics of Viscoelasticity: General Principles, Constitutive Modelling, Analytical and Numerical Techniques. Amsterdam, Elsevier, 1997.
5. Crochet M.J., Davies A.R., Walters K. Numerical Simulation of Non-Newtonian Flow. Amsterdam, Elsevier, 2012.
6. Aristov S.N., Skul'skii O.I. Exact Solution of the Problem of Flow of a Polymer Solution in a Plane Channel. Journal of Engineering Physics and Thermophysics, 2003, vol. 76, no. 3, pp. 577-585. DOI: 10.1023/A:1024768930375
7. Carrozza M.A., Hulsen M.A., Hutter M., Anderson P. D. Viscoelastic Fluid Flow Simulation Using the Contravariant Deformation Formulation. Journal of Non-Newtonian Fluid Mechanics, 2019, vol. 270, pp. 23-35. DOI: 10.1016/j.jnnfm.2019.07.001
8. Gamzaev Kh.M. Numerical Method of Pipeline Hydraulics Identification at Turbulent Flow of Viscous Liquids. Pipeline Science and Technology, 2019, vol. 3, no. 2, pp. 118-124. DOI: 10.28999/2514-541X-2019-3-2-118-124
9. Samarskii A.A., Vabishchevich P.N. Numerical Methods for Solving Inverse Problems of Mathematical Physics. Berlin, De Gruyter, 2007.
10. Vabishchevich P.N., Vasil'ev V.I., Vasil'eva M.V. Computational Identification of the Right-Hand Side of a Parabolic Equation. Computational Mathematics and Mathematical Physics, 2015, vol. 55, no. 6, pp. 1015-1021. DOI: 10.1134/S0965542515030185
11. Deng Z.C., Qian K., Rao X.B., Yang L., Luo G.W. An Inverse Problem of Identifying the Source Coefficient in a Degenerate Heat Equation. Inverse Problems in Science and Engineering, 2015, vol. 23, no. 3, pp. 498-517. DOI: 10.1080/17415977.2014.922079
12. Borukhov V.T., Zayats G.M. Identification of a Time-Dependent Source Term in Nonlinear Hyperbolic or Parabolic Heat Equation. International Journal of Heat and Mass Transfer, 2015, vol. 91, pp. 1106-1113. DOI: 10.1016/j.ijheatmasstransfer.2015.07.06
13. Ashyralyev A., Erdogan A.S., Demirdag O. On the Determination of the Right-Hand Side in a Parabolic Equation. Applied Numerical Mathematics, 2012, vol. 62, no. 11, pp. 1672-1683. DOI: 10.1016/j.apnum.2012.05.008
14. Gamzaev Kh. I. The Problem of Identifying the Trajectory of a Mobile Point Source in the Convective Transport Equation. Bulletin of the South Ural State University. Mathematical Modelling, Programming and Computer Software, 2021, vol. 14, no. 2, pp. 78-84. DOI: 10.14529/mmp210208
Received October 23, 2021
УДК 532.546+519.6 DOI: 10.14529/mmp220408
ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ НЕСТАЦИОНАРНОГО ТЕЧЕНИЯ ВЯЗКОУПРУГОЙ ЖИДКОСТИ В ТРУБЕ
А.Р. Алиев1,2, Х.М. Гамзаев2, А.А. Дарвиш3, Т.А. Нофал4,5
1 Азербайджанский государственный университет нефти и промышленности, г. Баку, Азербайджан
2Институт математики и механики НАН Азербайджана, г. Баку, Азербайджан 3Хелуанский университет, г. Каир, Египет 4Таифский университет, г. Таиф, Саудовская Аравия 5Университет Эль-Миния, г. Миниа, Египет
Рассматривается процесс нестационарного течения несжимаемой вязкоупругой жидкости в цилиндрической трубе постоянного сечения. Для описания реологических свойств вязкоупругой жидкости используется модель Кельвина - Фойгта и математическая модель данного процесса представляется в виде интегро-дифференциального уравнения в частных производных. В рамках данной модели поставлена задача определения перепада давления по длине трубы, обеспечивающего пропуск заданного расхода вязкоупругой жидкости по трубе. Поставленная задача относится к классу обратных задач, связанных с восстановлением правых частей интегро-дифференциальных уравнений. Путем замены переменных интегро-дифференциальное уравнение преобразуется в дифференциальное уравнение третьего порядка в частных производных. Сначала строится дискретный аналог задачи с использованием конечно-разностных аппроксимаций. Для решения полученной разностной задачи предлагается специальное представление, позволяющее расщепить задачи на две взаимно независимых разностные задачи второго порядка. В результате получена явная формула для определения приближенного значения перепада давления по длине трубопровода при каждом дискретном значении временной переменной. На основе предложенного вычислительного алгоритма были проведены численные эксперименты для модельных задач.
Ключевые слова: вязкоупругая жидкость; модель Кельвина - Фойгта; интегро-дифференциальное уравнение; перепад давления по длине трубы; обратная задача.
Вестник !Ю"УрГ"У. Серия «Математическое моделирование и программирование» (Вестник ЮУрГУ ММП). 2022. Т. 15, № 4. С. 90-98
Араз Рафиг оглу Алиев, доктор физико-математических наук, профессор, заведующий кафедрой, кафедра «Общая и прикладная математика:», Азербайджанский государственный университет нефти и промышленности (г. Баку, Азербайджан); главный научный сотрудник, отдел «Функциональный анализ», Институт математики и Механики НАН Азербайджана, (г. Баку, Азербайджан), [email protected].
Ханлар Мехвали оглу Гамзаев, доктор технических наук, профессор, кафедра «Общая и прикладная математика», Азербайджанский государственный университет нефти и промышленности, (г. Баку, Азербайджан), [email protected].
Адель Абделфаттах Дарвиш, кандидат физико-математических наук, профессор, кафедра «Математика», Хелуанский университет (г. Каир, Египет), [email protected].
Тахер Абдехамед Нофал, кандидат физико-математических наук, профессор, кафедра «Математика», Таифский университет (г. Таиф, Саудовская Аравия); профессор, кафедра «Математика», Университет Эль-Миниа, (г. Миниа, Египет), [email protected].
Поступила в редакцию 23 октября 2021 г.
