Mathematical modeling of oscillatory processes with a free boundary Текст научной статьи по специальности «Математика»

UDC 519.6

Mathematical modeling of oscillatory processes with a free boundary* A.E. Chistyakov 1, E. A. Protsenko 2, E.F. Timofeeva 3**

1 Don State Technical University, Rostov-on-Don, Russian Federation

2 Rostov State Economic University, Rostov-on-Don, Russian Federation

3 North-Caucasus Federal University, Stavropol, Russian Federation

Introduction. The paper is devoted to research of the wave processes with free boundary based on the finite-difference method.

Materials and methods. A mathematical model describing the dynamics of distribution of wave fluctuation was proposed on the basis of heterogeneous wave equation with the appropriate initial and boundary conditions. Discretization of the model was conducted using the integro-interpolation method taking into account the partial "filling" of computational cells. The adaptive modified alternating triangular iterative method of variational type with the highest rate of convergence in the class of two-layer iterative methods for solving the developed difference equations. Results. The developed discrete mathematical model for numerical simulation of wave propagation. The results of numerical experiments were obtained. The developed numerical algorithms and their computer implementation were used to research the dynamics of distribution of wave processes in the presence of the free boundary.

Discussion and conclusions. The obtained results can be used for research of the dynamics of distribution of the wave processes with a free boundary and controlling in conducting experimental researches, evaluation and diagnosis, etc.

Keywords: wave oscillations, grid equations, adaptive modified alternating triangular iterative method.

Introduction. The improving of monitoring systems and increasing of detailing the experimental information lead to the need of the consideration of wave fields with more precision at present. Despite the large number of studies in this direction, the problem of mathematical modeling of the distribution of wave distribution is relevant because of the wide diversity of problems, the specifics of which should be considered in the development of methods and algorithms for constructing numerical solutions and their computer implementation. Methods of mathematical modeling using high-performance computers are the most effective approach for obtaining information about research processes. Theoretical and experimental researchers of these phenomena are limited to methodological and technical difficulties.

Traditional methods, used at modeling processes of wave distribution, are the following: asymptotic (radial) methods; integral methods (based on the Huygens' principle); direct numerical methods.

It should be noted that the characteristic feature of radial methods (method of zero radial approximation, matrix methods, and methods of the generalized ray) is to research a limited part of the wave field. The wave field is represented as sum of waves, which are distributed with the local

*

The paper was supported by RFFR project № 15-07-08626, № 15-01-08619.

** E-mail: cheese 05@mail.ru, eapros@rambler.ru, teflena@mail.ru

velocity along the ray trajectories imposed by law of Snellius. The amplitude of these waves is determined by the geometric divergence of rays in the path from source to receiver. In Haskell-Thomson matrix method and it's modifications, the problem is investigated in the frequency domain: when a set of frequencies solve the system of Helmholtz equations is solved with a further shift in the time domain at a certain set of frequencies.

We can calculate the total wave field using direct numerical methods for solving the wave equation. The modeling of total wave field is performed at solving differential equations of wave motion. The wave field is calculated at a set of nearby discrete grid nodes by approximation of derivatives by finite differences and recursive solution of differential equation.

The purpose of this paper is mathematical modeling and development of the software complex, implemented the model and intended to describe the distribution of wave fluctuations with free boundary based on the finite-difference method.

The following problems were solved in accordance with the purpose of this paper:

- the discrete model taking into account the fullness of computational cells was developed which guaranteed the implementation of the basic conservation laws at the discrete level;

- the dependence of the approximation error from step of temporary variable was researched;

- the conditions of stability of three-layer difference schemes were obtained;

- the optimal values for the weight parameter of three-layer difference schemes were calculated;

- the variant of adaptive modified alternating triangular iterative method of variational type (MATM), which has the best rate of convergence, was developed;

- the software complex for modeling the distribution of oscillatory processes with a free boundary was designed.

Materials and methods. Problem statement. We must solve the heterogeneous wave equation [1-3]: p"t = c2Ap + f with the following initial conditions:

p(x,y,0) = (p0{x,y), p't(x,y,0) = q\(x,y)

And boundary conditions:

- on the solid boundary: p(x, y, t) = 0, (x, y) e y,

- on the soft boundary: p'n(x, y, t) = 0, (x, y) e y,

- on the free boundary: p't = cp'n, (x, y) e y,

where n is the internal surface normal.

The computational domain is inscribed in a rectangle and covered by a uniform computational grid (o = atxaxxay : ot = {tn = nht ,0 < n < Nt, lt = htNt}, &x= (x, = ihx ,0 < i < N, h = \NX},

(°y= (yj = jhy,0 < j <N ly = hy (N _l)}, where l, j, n -are indexes by the time and spatial

directions Ox, Oy, accordingly; ht, hx, hy are steps by time and spatial coordinate directions; Nt, Nx, Ny is the number of nodes by the time variable and spatial coordinate directions; lt, lx, ly is the length of computational domain by the time variable and spatial coordinate directions.

Discrete model. The approximation of model equations was performed on rectangular grids using the modification of integro-interpolation method [4, 5] that takes into account "filling" of

computational cells. So, we can dynamically vary the discrete computational domain without transforming the grid, i.e. without additional computational costs [6, 7].

The approximation of operator of the second differential derivative has the form:

r

(q\j № ) x =(q v+1/2 ,j i++ i,J - (\j u- li] h2 -(qL-(q2)uj ach+Jx, i=12...,Nx -1, j=12...,Ny -1,

in the case of third order boundary condition u'n (x, y , t} = au + J.

"Filling" coefficients of control domain qm, m = 0,4 can be calculated by formulas:

in) in) -°',J +0i+1J +q+l,j+1 + 0,J_^+1,j + 0+1,j+1

(qm )i,J = s ' (qo )i J = 4 ' (q1 )i,J = 2 ,

o.+o,, , , 0.,,.,,+a/ x 0,J + o+1,j

j 2 v 'i, j 2 v 7'', j 2 o,, - the "filling" of cell (i, j) , i = 0,1,...,Nx -1, j = 0,1,...,Ny -1.

The location of the nodes relative to cells is shown in Fig. 1.

We used the integro-interpolation method for obtaining the discrete model. The difference scheme, approximated the wave equation, is in the form:

Pi j1 - 2Ph + Pu1 _ 2 P+W - 2Pu + P^j , 2 P,j+1 - 2P + Pj1

= ,+,J ^ ,I-,J + ,J+1 , 2,J , -1 + fn. n = 0.1..... Nt -1,

h2 h2 h2 ht hx hy

i = 1,2,...,Nx -1, j = 1,2,...,Ny -1,

where p t j = a1 pij1 +(1 -a1 -a2 ) pnUj +^2 p^j1 are weights of scheme.

Fig. 1. Location of nodes relative to cells

Discrete analog of wave equation with boundary conditions in the Dirichlet's form (p=0) is following:

PiJ -2Ph+PU ,n 2 ■ ( \Pi+1,i - 2Pij+Pi-1,

+qo,ij-j2—- = qoijfij +c min(qvj.q2,,j)-

2

h2 '.j h2 .'. J . J v 1 '. J .'. J/ K

2

i+1 - 2P,

+c" min (q3ij, -^ Fi,J-1, n = 0,1,..., N -1, i = 1,2,..., Nx -1, j = 1,2,..., Ny -1.

hy

Discrete analog of wave equation with boundary conditions in the Neumann's form (p=0) is following:

Pu - 2Ph + Pu _ 2 Pi+u - Pu 2 Pi,j - P.-U , 2 p,j+1 - pi,j

q0, i, j" h?-- = C ^^hT^ - C + c ^^--

t x x y

-c\ii,j Pi,J -2Pi-j~1 + , n = 0,1,...,Nt -1, i = 1,2,...,Nx -1, j = 1,2,...,Ny -1.

hy

The condition of free boundary for problem of distribution of wave fluctuations takes the form:

Pu - Pil.Ji+U - P,j

= с , i = 0, j = 1,2,...,Ny -1,

ht К

pn+1 - pn. p +1 - p

p,j = +1 p,j, j = о, . = 1,2,...,n -1.

ht hy

Discrete model research. We analyzed the wave equation c"tt = div(^grad (c)) with initial conditions c[=0 = c0, c;|(=0 = q.

Approximation of problem by spatial variables can be defined in the form: c"t=-Ac или

r

"t = Z (^т) , where r is a space dimension.

i=i x

The analytical solution in the basis, composed of eigenvectors, has the form:

c = ZaX'' AXi = AiXI ' a. (t) = ai,0 COS (№) + "7= Sln (^) •

i

ctt =

Symmetric scheme for the wave problem takes the form:

'— = -Xpa^1 - \ (1 - 2a) x1 - Xlaan-1, n > 1;

a"+1 - 2an + an-1

r2

,n+1 n

---^ = -AC+1 - ^ (1/2 -a)a,", n = 0.

T T

Numerical solution is corresponded to the following: a" = A cos(n^) + B sin(n^), where cos(p = k/2, k = 2- Att2 /(1 + A&T2).

The stability condition of difference scheme taking into account the estimation of the maximum condition number Amx <4rj/h2 has the form: z<h((1 -4a)rj Numerical solution is represented in the form:

-1/2

an =x 0COS((p/T)t) + ^ sin((p/r)t), p:

, v ; ^(1 + 1/ 4)r2) v ;

The approximation error of the fluctuation frequency is following:

\ Г /2 Л

v 1 + j

L Ar2/2 1

arceos

1 -

1 + Aar

y

Ar

1 + A<T2 +

f Ar2/2 1

x 1 + A.ar ,

V ï y

+

O (r )

= 1+O (r).

r^A r^A r^A 24 '

According to the obtained estimation, we concluded that the fluctuation frequencies are slightly differed from the actual values and depend on the time step related to the wave period and weight in computing the distribution of oscillatory processes on the basis of differential methods.

The dependence of approximation error of the fluctuation frequency yi r^fA -1 from t^JA is

shown in Fig. 2. This parameter describes the time step, referred to the wave period; the unit interval is corresponded to the 1/2n of wave period.

0.05

0.05

-0.1

CT=0 CT=1/1 2

0.2

0.4

0.6

0.8

Fig. 2. The dependence of approximation error of the fluctuation frequency from parameter describing

the time step, referred to the wave period

1 -

1+<y

We considered the function f (y, c) = (1 i y) arccos The dependence of the function f on y is shown in the Fig. 3.

0.003 0.001

-0.001

f y2/2 1 r-

-1, where y = t^JAí < 1.

G= 1/12

G=0.0869172

0.2 0.4 0.6 0.8 Fig. 3. The dependence of the function f on y

Optimal values from the point of view of preserving values of the frequency of distribution of oscillatory processes depended on the constraints to the y and deviation of the fluctuation frequency from the actual values max(f are given in Table 1.

Table 1

Deviation of the fluctuation frequency from the actual values depended on the weight of the scheme

Value of weight parameter a Constraints to the y Value off,%

1 0.08375 0.34667 0.0005198

2 0.084 0.43794 0.001329

3 0.0842042 0.5 0.002266

4 0.0845 0.57782 0.004061

5 0.085 0.68882 0.008268

6 0.08531575 0.75 0.01168

7 0.0855 0.78332 0.01394

8 0.086 0.86674 0.02107

9 0.0865 0.94204 0.02964

10 0.0869172 1 0.0379

The approximation error of the second derivative differential operator based on schemes of the second and fourth order accuracy. Difference schemes of the second and fourth order accuracy were obtained in paper [9].

The scheme of the second order of accuracy: ax (r) = 1 - 2 (l - cos (n / r)) / (n / r)2.

The scheme of the fourth order of accuracy:

a2( r) = 1 - (l5 -16 cos (n / r ) + cos (2n / r)) /6 (n / r )2.

The value r describes the number of nodes for the half of wave period (for a description of the object). Based on the obtained estimations we can calculate numerical values of winnings in computational time using schemes of high order of accuracy.

The dependence of approximation error of the second derivative operator for schemes of second and fourth order of accuracy are shown in Fig. 4.

0 1 0.01

xl0~3

1x10" 5 1x10" 6

1x10"7 3

1 10 100 1x10

Fig. 4. The dependence of approximation error of the second derivative operator: 1 - schemes of fourth order of accuracy, 2 - schemes of second order of accuracy

Modified alternating triangular iterative method for calculating grid equations with non-self-adjoint operator. Using schemes taking into account the "filling" of control domains requires the implementation of the modified variant of alternating triangular iterative method [5-8, 10]. This method has the quality evaluation of convergence rate and is effective for solving problems on grids taking into account the complex geometry of researching objects. At present the modified alternating triangular iterative method (MATM) is widely used at solving problems of aero-hydrodynamics and transport of bottom materials [11-16].

We use the implicit iterative process for the solution of grid equations [5, 7-9]:

m+1 _ m

B--— + Axm = f, B : H ^ H,

where A is the linear, positive definite operator; m is an iteration number; t>0 is an iterative parameter; B is some conversion operator. Note that the conversion of the operator B should be significantly easier than the direct conversion of the source operator A. We assume the additive decomposition of the operator A in the construction of operator B:

A =(A + A*)/2 = R + R, R = R*, A =(A_A*)/2.

The operator-stabilizer is defined as follows:

B = (D + oR )D_1(D + aR), D = D* > 0, o > 0, y e H, where D is some operator.

The value ro is minimal at c =

( Dwm, wm )

where wm is the correction vector.

V( D-%wm, Rwm )'

Iteration parameters for MATM of minimal corrections are defined by the formula:

(Awm, wm )

t -^, Bwm = Axm - f, m = 0,1,...

m+1 (B- Awm, Awm ), J

The algorithm of the modified alternating triangular iterative method of minimal corrections is in the form:

r'" = Ax'" -/, B{com)w"'=r"\ r5„,=

V(/) ]R2w"',R2W"')' (K ]Ay,A,y)(KW'",w"')'

(B-lAwm,Awm) f I s 2k ^ // xx (Awm,wm)

k = V ' Q = 1 Sm km / /l , k /l _ _ 2\\ = g > J

^ (B_1Awm , Awm ) , Q + km )J/(1 + km (1 )) , ^ Q ^^m , ^ ) ,

The estimation of convergence rate of this method is in the form:

* — 1 2 ||z"+1|| z"||, =v{JTTl +4k ) ,

( b- A®m, A®m)

where v is the condition number of the matrix Co at k = y-T < 1.

( B- A com, A com )

The additional apriory information about initial problem is the important at this approach. This information for MATM method is associated with the estimations ô and A:

D <1R1D1R2 A a 4

The estimation of the condition number is the following: v <

1 + # e_S

Software implementation. The developed software complex for calculating the distribution of oscillatory processes with a free boundary consist of the following units: initial data input unit; geometry calculation unit; unit for calculation coefficients of grid equations; unit for calculation functions of right parts of grid equations; transition unit for a coarser grid; unit for calculation grid equations using the modified alternating triangular method; calculation unit for computational window; unit of account of boundary conditions; output unit of computing pressure functions; output unit of spectrum; unit of calculating phase; unit of calculating the phase gradient; calculating orientation unit; output orientation unit.

Results. The developed numerical algorithms and software complex that implements them were used for research of the dynamics of wave processes distribution in the presence of a free boundary.

The dynamics of wave processes distribution in the presence of a free boundary is given in Fig. 5. The bottom and left boundary are free, the top and right boundary are solid.

Fig. 5. Distribution of wave processes with free boundary

Discussion and conclusions. The paper covered the research of wave fluctuations and development software complex designed for describing wave radiation processes with a free boundary. The proposed mathematical model is based on the heterogeneous wave equation with appropriate initial and boundary conditions. The grid method was used for solving problem. The discrete model was constructed using the integro-interpolation method taking into account the "filling" of computational cells which guaranteed the performing the basic conservation laws (for flux of the electric field and the circulation of magnetic field) on the discrete level. Optimal values of weight parameters were calculated.

The obtained grid equations were solved using the adaptive modified alternating triangular iterative method of variational type that has the quality evaluation of convergence rate in the class of two-layer iterative methods. The software complex for modeling distribution of electromagnetic waves in waveguides with the complex geometry was designed on the basis of developed parallel algorithms for adaptive MATM. The developed numerical algorithms and software complex implemented of them were used for researching the dynamics of wave processes distribution with the free boundary.

References

1. Markov, G.T., Chaplin, A.F. Excitation of electromagnetic waves. Moscow, Radio and communication, 1983, 296 p.

2. Vladimirov, V.S. Equations of mathematical physics. The textbook for phys. and mechanical-mathematical. spec. universities. 4-th ed. Moscow, Nauka, 1981, 512 p.

3. Sukhinov, A.I., Zuev, V.N., Semenisty, V.V. Equations of mathematical physics. Taganrog, TRTU, 2005.

4. Samarskii, A.A. The Theory of Difference Schemes. Moscow, Nauka, 1989.

5. Samarskii, A.A., Nikolaev, E.S. Numerical Methods for Grid Equations. Moscow, Nauka,

1978.

6. Sukhinov A.I., Chistyakov, A.E., Shishenya, A.V. Error estimate for diffusion equations solved by schemes with weights. Mathematical Models and Computer Simulations, 2014, vol. 6, iss. 3, pp. 324-331.

7. Sukhinov, A.I., Chistyakov, A.E., Fomenko, N.A. Method of construction difference scheme for problems of diffusion-convection reaction, takes into the degree filling of the control volume. Izvestiya SFedU, Engineering Sciences, 2013, no. 4 (141), pp. 87-98.

8. Konovalov, A.N. To the Theory of the Alternating Triangle Iteration Method. Sib. Math. J., 2002 no. 43, pp. 552-572.

9. Sukhinov, A.I., Chistyakov, A.E., Yakobovskii, M.V. Accuracy of the numerical solution of the equations of diffusion-convection using the difference schemes of second and fourth order approximation error. Bulletin of the South Ural State University, series "Computational Mathematics and Software Engineering", 2016, vol. 5, iss. 1, pp. 47-62.

10. Sukhinov, A.I., Chistyakov, A.E. Adaptive modified alternating triangular iterative method for solving grid equations with a non-self-adjoint operator. Mathematical Models and Computer Simulations, 2012, vol. 4, iss. 4, pp. 398-409.

11. Sukhinov, A.I., Chistyakov, A.E., Timofeeva, E.F., Shishenya, A.V. Mathematical model for calculating coastal wave processes operator. Mathematical Models and Computer Simulations,

2013, vol. 5, iss. 2, pp 122-129.

12. Sukhinov, A.I., Chistyakov, A.E., Protsenko, E.A. Two-Dimensional Mathematical Model of Sediment Transportation. Izv. Yuzhn. Federal. Univ., Tekh. Nauki, 2011, no. 8, pp. 32-44.

13. Sukhinov, A.I., Chistyakov, A.E., Protsenko, E.A. Mathematical modeling of sediment transport in coastal aquatic systems on a multiprocessor computer system. Computational methods and programming, 2014, vol. 15, pp. 610-620.

14. Sukhinov, A.I., Protsenko, E.A., Chistyakov, A.E., Shreter, S.A. Comparison of computational efficiencies of explicit and implicit schemes for the sediment transport problem in coastal zones. Numerical Methods and Programming, 2015, vol. 16, no. 3, pp. 328-338.

15. Nikitina, A.V., Semenyakina, A.A., Chistyakov, A.E., Protsenko, E.A., Yakovenko, I.V. Application schemes of higher order accuracy for solving biological kinetics on a multiprocessor computer system. Basic Research, 2015, no. 12-3, pp. 500-504.

16. Sukhinov, A.I., Chistyakov, A.E., Protsenko, E.A. Mathematical modeling of sediment transport in the coastal zone of shallow reservoirs. Mathematical Models and Computer Simulations,

2014, vol. 6, no. 4, pp. 351-363.

Authors:

Chistyakov Alexander Ev., Professor, Department "Software computer engineering and automated systems", Don State Technical University (1 Gagarin sq., Rostov-on-Don, Russian Federation), Doctor of Science in Physics and Maths.

Protsenko Elena An., Associate Professor, Mathematics Department, Taganrog Institute named after A. P. Chekhov (branch) of Rostov state University (RINE) (48 Initiative St., Taganrog, Russian Federation), Candidate of Science in Physics and Maths.

Timofeeva Elena F., Associate Professor, Department of Applied Mathematics and Computer Technology, North Caucasus State Technical University (2, Kulakov Pr., Stavropol, 355029, Russian Federation), Candidate of Science in Physics and Maths.

УДК 519.6

Математическое моделирование колебательных процессов со свободной

А

границей

А.Е. Чистяков , Е.А. Проценко , Е.Ф. Тимофеева

1 Донской государственный технический университет, г. Ростов-на-Дону, Российская Федерация

2 Ростовский государственный экономический университет, г. Ростов-на-Дону, Российская Федерация

3 Северо-Кавказский федеральный университет, г. Ставрополь, Российская Федерация

Введение. Статья посвящена исследованию процесса волновых колебаний и разработке комплекса программ, предназначенного для описания распространения волновых процессов. Целью настоящей работы является математическое моделирование и разработка реализующего модель, комплекса программ, предназначенного для описания распространения волновых колебаний со свободной границей на основе конечно-разностного метода.

Материалы и методы. Предложена математическая модель, описывающая динамику распространения волновых колебаний, в основе которой лежит неоднородное волновое уравнение с соответствующими начальными и граничными условиями. Дискретизация модели проведена интегро-интерполяционным методом, при этом реализован подход, учитывающий частичную «заполненность» расчетных ячеек. Для решения, полученных сеточных уравнений применен адаптивный модифицированный попеременно-треугольный итерационный метод вариационного типа, имеющий наиболее высокую скорость сходимости в классе двухслойных итерационных методов.

Результаты исследования. Разработаны математическая модель и программное обеспечение для численного моделирования распространения волн. Приведены результаты численных экспериментов. Разработанные численные алгоритмы и их компьютерная реализация использованы для исследования динамики распространения волновых процессов при наличии свободной границы.

Обсуждение и заключения. Полученные результаты могут быть использованы в процессе исследования динамики распространения волновых процессов при наличии свободной границы, контроля при проведении экспериментальных исследований и т.д.

Ключевые слова: волновые колебания, сеточные уравнения, адаптивный модифицированный попеременно-треугольный итерационный метод

Авторы:

Чистяков Александр Евгеньевич, профессор кафедры «Программное обеспечение вычислительной техники и автоматизированных систем» Донского государственного технического университета (РФ, 344000, г. Ростов-на-Дону, пл. Гагарина, 1), доктор физико-математических наук.

Проценко Елена Анатольевна, доцент кафедры «Математика» Таганрогского института имени А.П. Чехова (филиал) РГЭУ (РИНХ) (РФ, Ростовская область, г. Таганрог, ул. Инициативная, 48), кандидат физико-математических наук.

Тимофеева Елена Федоровна, доцент кафедры «Прикладная математика и компьютерные технологии» Северо-Кавказского государственного технического университета (РФ, 355029, г. Ставрополь, просп. Кулакова, 2), кандидат физико-математических наук.

*

Работа выполнена при финансовой поддержке РФФИ по проектам № 15-01-08619, № 15-07-08626. Авторы: e-mail: capros@ramblcr.ru. teflena@mail.ru. cheese 05@mail.ru