Аналитическое решение задачи о малых вынужденных колебаниях идеальной жидкости

2017 Математика и механика № 48

УДК 532.5.031

Б01 10.17223/19988621/48/7

А.В. Мерзляков, З.О. Матыева

АНАЛИТИЧЕСКОЕ РЕШЕНИЕ ЗАДАЧИ О МАЛЫХ ВЫНУЖДЕННЫХ КОЛЕБАНИЯХ ИДЕАЛЬНОЙ ЖИДКОСТИ

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

Ключевые слова: идеальная жидкость, потенциал скорости, уравнение Лапласа, метод разделения переменных.

В настоящее время в практических целях возникает задача о движении жидкости в сосуде (например, колебание топлива в баке автомобиля или ракеты, движение жидкости в автомобильной или железнодорожной цистерне и т.д.). Во многих случаях жидкость можно считать идеальной. Решение задачи о движении идеальной жидкости в сосуде при наличии свободной поверхности является важной технической задачей.

Задача о малых колебаниях идеальной жидкости в ограниченном объеме со свободной поверхностью под действием силы тяжести многократно рассматривалась в различных работах. Основные результаты аналитических решений этой задачи приведены в работах [1-8]. Есть и иностранные работы по этой теме, например [9-11]. Следует, однако, отметить, что во всех перечисленных работах рассматриваются только свободные колебания жидкости; кроме того, практически во всех перечисленных работах ограничиваются только определением частот колебаний и не приводят динамику изменения формы свободной поверхности. Появление современных компьютеров позволило снять этот недостаток. Сейчас можно аналитически решать задачи на определение формы свободной поверхности и пользоваться вычислительными средствами для реализации полученных результатов. Также кроме свободных колебаний оказалось возможным рассматривать вынужденные колебания идеальной жидкости.

В данной работе рассматривается один из возможных случаев движения идеальной жидкости - колебательное движение под действием переменного наружного давления и силы тяжести. Задача решается аналитически в линейном приближении.

Математическая постановка задачи

Рассмотрим происходящее под действием сил тяжести и переменного наружного давления волновое движение однородной несжимаемой идеальной жидкости, ограниченной снизу и с боков некоторыми неподвижными поверхностями, а сверху свободной поверхностью. Целью работы является определение формы свободной поверхности жидкости в любой момент времени. Движение жидкости

начинается тогда, когда имеет место некоторое возмущение жидкости. Пусть возмущение жидкости обусловливается причинами, действующими исключительно на ее свободную поверхность. Движение идеальной несжимаемой жидкости в плоском случае описывается уравнениями Эйлера (см. [12]):

dvx

dt

= X -

1 dp

р dx'

dv

= 7 -1—, it р dy

V. + ^ = 0,

dx dy

где ух,и уу - проекции вектора скорости жидкости на оси декартовой системы координат; X и У - проекции внешних сил на оси Ох, Оу соответственно; Р - давление, р - плотность, t - время.

В силу потенциальности движения идеальной жидкости можно ввести потенциал скорости жидкости и, для которого V = Уи . Тогда вместо уравнений Эйлера для описания движения жидкости можно использовать уравнение Лапласа для потенциала и:

d 2u

d 2u

Ди =-+-= 0.

дх2 дУ2

Граничными условиями для него являются:

- на твердой стенке - условие непротекания:

— = 0 (п - направление нормали к границе);

дп

- на свободной поверхности - интеграл Коши - Лагранжа:

ди V2 + уУ Р „

— +-- + — + П = [ (t),

дt 2 р

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

Рассмотрим движение идеальной жидкости в прямоугольном сосуде шириной а и высотой И

Движение жидкости будет определяться путем решения уравнение Лапласа (1) для потенциала скорости в области, занятой жидкостью внутри сосуда. Граничными условиями для него являются:

- на твердых границах Г и ГЗ -

ди =

условия непротекания — = 0 ;

дх

- на твердой границе Г2 - условие

ди =

непротекания — = 0 .

дУ

S

Oi

Гз

O

Г2

Рис. 1. Расчетная область задачи Fig. 1. Computational domain of the problem

На свободной поверхности в качестве граничного условия ставится интеграл Коши - Лагранжа со следующими дополнениями. Единственной массовой силой, действующей на жидкость, является сила тяжести. Если ввести величину г - отклонение точки свободной поверхности от равновесного положения (соответствующая координата представлена на рис. 1), то потенциал силы тяжести, действующей на единицу массы, определяется формулой П = gz.

При решении поставленной задачи предполагается, что отклонение свободной поверхности жидкости от положения равновесия настолько мало, что область, занятая жидкостью, сохраняет прямоугольную форму. Из-за этого в формуле можно пренебречь квадратами скоростей. Считается, что давление над поверхностью жидкости Р является некоторой функцией времени t и координаты х. Произвольную функцию / (0 можно сделать равной 0. Поэтому граничное условие на свободной поверхности примет следующий вид:

ди Р

— + + gz = 0.

дt р

В качестве начальных условий задачи задаются начальная форма свободной поверхности и начальная скорость точек свободной поверхности. Начальная форма свободной поверхности в рассматриваемой задаче считается равновесной, т.е.

г(0,х) = 0.

Начальная скорость точек свободной поверхности жидкости считается равной нулю, т.е.

ди ди 0

ду у=Ь,1=0 дг г=0Д=0

Метод решения

Уравнение Лапласа (1) решается методом разделения переменных (см. [13]), согласно которому потенциал вектора скорости представляется в виде следующего произведения:

и = Т ^) X (х)У (у),

где Т(^ - функция, зависящая только от времени, X(x) - функция, зависящая только от координаты х, У(у) - функция, зависящая только от координаты у.

После подстановки этого произведения в уравнение (1) имеем

ТХ "У + ТХУ" = 0.

Деление этого уравнения на произведение ТХУ приводит его к следующему виду:

XI = -Ц =-Х 2

X у ' где X2 - некоторое положительное число.

Из данного равенства получаются дифференциальные уравнения для определения функций Х и У.

Уравнение для X - задача Штурма - Лиувилля - сводится к следующему:

X " + Х2 X = 0.

Решение данного уравнения выглядит так:

X = A cos (Xx) + B sin(Xx).

Постоянные А и В определяются из граничных условий, которые заданы на границах Ti и Г3. Согласно им,

ди dx

x=0

= X '(0) = 0; dx

dx

= X'(a) = 0 .

Из первого условия следует В = 0; из второго условия

X. = ™ ,

где п - целое число. Поэтому выражение для функции Х - собственная функция задачи Штурма - Лиувилля - выглядит так:

Xn = c0s |™x

Уравнение для У сводится к следующему:

У "-Х2У = 0 Решение данного уравнения имеет вид

У = А сИ(Ху) + Вsh(Ху). Согласно граничному условию на границе Г2,

du ¥

= Y '(0) = 0.

y=0

Из этого условия следует В = 0. Поэтому выражение для функции У

Уп = сИ (Хпу) = сИ [Ну) .

Таким образом, для произвольно заданного натурального номера п решение уравнения Лапласа будет выглядеть так:

а общее решение - так:

, s Л (nn Л fnn Л

un=тп(t )ch i -^y Icos \~x I'

= ZTn (t)ch(ППУ 1 cos(nnx|.

n=0

Для поиска функций Г(0 необходимо использовать граничное условие на свободной поверхности, которое нужно продифференцировать по времени и учесть,

дя ди

что — = уу = —. Тогда выражение (2) после небольшого преобразования примет

вид

д 2u

It2

ди

~ду

y =h

i дР р dt

= Q (x, t).

a

Подстановка выражения для потенциала в это равенство и разложение в ряд Фурье правой части по собственным функциям задачи Штурма - Лиувилля приводит к следующему выражению:

^ rr,„ , (nn , A (nn Л ^ „ nn . inn , Л (nn Л ^ ^ , ч (nn Л X Tn ch I — h I cos I — x\ + g X Tn— sh ( — h I cos I — x I = X Qn (t )cos I — x I.

n=0 I a J I a J n=0 a { a J { a J n=0 { a J

После перенесения всех слагаемых в одну сторону, группировки и вынесения общих множителей имеем

X{t; ch {™h gTn nn sh {™h J-Qn (t )}cos (nnx | = 0.

n=01 { a J a { a J J { a J

Данная сумма будет равна нулю, когда каждая фигурная скобка будет равна нулю. Это приводит к следующей системе уравнений для определения функций T :

± n

__ , (nn , Л ^ nn , (nn , Л ^ / Ч

Tn ch I —h I + gTn— sh I —h j -- Qn (t ) = 0 .

{ a J a { a J

Небольшое преобразование и введение новых обозначений приводят систему к виду

n = 0, T'' = Q0 (t);

n > 0, ТЩ + v2nTn = Fn (t),

2 nn . f nn , Л „ , ч Qn (t) где = g—thI h I; Fn (t) = *n}>

0 V 0 J Ch\~a~hI

ch I — f

V a

Процедура решения уравнений проводится согласно методике, приведенной в работе [13]. Эта методика основана на том, что решение уравнения ищется в виде суммы общего решения однородного уравнения и частного решения неоднородного. Общее решение однородного уравнения соответствует свободным колебаниям жидкости и рассматриваться не будет, так как начальное состояние жидкости - равновесие, поэтому свободные колебания отсутствуют. Частное решение, соответствующее вынужденным колебаниям, определяется по формуле

t

Tn

0

где fn(t) - решение соответствующего однородного уравнения при следующих начальных условиях:

fn (0) = 0; f'(0) = 1.

Для n = 0 выражение (4) приобретает вид

t

t0 (t) = ¡(t-t)q (T)dt ,

0

а для случая n > 0 выражение

Tn (t ) = — j sin (con (t-T))Fn (t) d T .

Вычисление полученных интегралов и подстановка их в формулу для потенциала (3) дает окончательный вид решения уравнения Лапласа. Опираясь на это решение, можно получить выражение для формы свободной поверхности жидкости через отклонение z от состояния равновесия. Это выражение вытекает из формулы для граничного условия (2):

z = -1 Г* + PZÜ1.

g U р ) W

Результаты проверки метода

1. Вынужденные колебания идеальной жидкости под действием гармонически меняющегося давления

В данной части работы рассматривается вынужденное колебание идеальной жидкости под действием переменного наружного давления, которое описывается выражением

P

= P0 sinj sin(ю0/) ,

р

где Р0 - отношение амплитудного значения наружного давления к плотности жидкости, ю0 - частота колебаний наружного давления. Как указано выше,

ч 1 dP( x, t) . (nx Л , ,

Q (Лt) =----— = -Pc sinI — IcosW).

р dt v a )

Разложение в ряд Фурье по собственным функциям задачи Штурма - Лиувил-ля дает выражения для коэффициентов Qn:

Qn (t) = 02 f P0ro0 sin(nX| cos (<V))|cos ¡V^fx) dx .

Окончательно формулы для коэффициентов выглядят так:

a) n = 0: Q0 (t) = - 2^^ ) .

n

b) n = 1: Qi (t) = 0;

(-1)n +1

c) n > 1: Qn (t) = 2Pomo / 2 , ,cos Kt) = Gn cos Kt),

n

(-1)n +1 где Gn = 2Po Ю 2 .

n

Согласно изложенной выше методике, получаются выражения для Tn:

г 2P

a) n = 0: To (t) = f (t - т) Qo (т) dт =--(1 - cos (coot));

о ПЮ0

b) n = 1: T (t) = 0;

c) n > 1 Tn (t) = -L j sin (On (t - T)) Qnjn\ d T =

Юп 0 ch[^h]

ch I — h V о

cos

-n --0

Окончательное выражение для потенциала скорости жидкости

2 P G

u =--— (1 - cos (-01))+ X

П-0

f nnh Л

cos (-0Q-c°s (-nt) ch f ™ , I cos f ™ x I.

n=2 ch I -

V a J

-n --

0

V a

Подстановка в (5) дает выражение для определения отклонения свободной поверхности г от состояния равновесия:

1

. fn Л • / . ч 2 P0 . / ч -n sin (-nt) — -0 sin (-01) fnn

P0sinV"X Jsin(-01)--0sin(-0t) + X G

n=2

-n —-0

По приведенной формуле был проведен расчет формы свободной поверхности идеальной жидкости с помощью программы Mathcad. Размеры сосуда: ширина а = 1 м, глубина h = 1 м. Плотность жидкости было взята 1000 кг/м3. Амплитудное значение давления 100 Па. Период колебаний вынуждающего давления составляет Т0 = 2.1 с.

На рис. 2 приведена форма свободной поверхности жидкости в разные моменты времени, соответствующие началу колебаний, четверти периода, трем четвертям периода и первому периоду колебаний вынуждающего давления. Величина z указана в метрах.

0.01

_z

♦ ♦♦♦ z

---z

z

У'З-------------------

-0.01 -

z

0

0 0.n 0.4 0.6 0.8 x

Рис. 2. Форма свободной поверхности жидкости для непрерывно меняющегося внешнего давления

Fig. 2. Shape of the free surface of liquid in the case of external pressure varying continuously

2. Вынужденные колебания идеальной жидкости под действием импульса давления на свободной поверхности

Вынужденные колебания идеальной жидкости происходят под действием импульса наружного давления на свободной поверхности, который описывается следующим выражением:

P (x, t) fn( x - b )Л fnt Л

v ' = P0 sin [—^—^ J sin [—J при b < x < c, t < T,

P

—1—1 = 0 - в остальных случаях.

P

Здесь b и c - границы зоны действия импульса давления, Т - его длительность.

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

\ P0n . f п( x - b )Л fnt Л , ^

Q (x, t) =—TT—sin I—i—^ j cos I — J при b < x < c, t < T,

Q (x, t ) = 0 - в остальных случаях.

Так же, как и ранее, необходимо эту функцию разложить в ряд Фурье по собственным функциям задачи Штурма - Лиувилля. При разложении необходимо иметь в виду, что во время действия импульса интеграл достаточно вычислять только по интервалу его действия, а после окончания действия импульса функция Q = 0. Поэтому коэффициенты разложения принимают вид

T

ft (, 2 J (-T S11, j cos (t )j cos ) л= - p cos (

при t < T, Qn (t) = 0 - в остальных случаях. Для экономии места здесь введено обозначение вычисленного интеграла:

qn =

. í 1\ fnn(c-b)Л fnn(c + b)Л 4a (c - b )cos I—^- j cos I—^- j

n-

При реализации формулы (4) необходимо иметь в виду, что частные решения отличаются в зависимости от того, как соотносится временной промежуток интегрирования t и время действия импульса T. Поэтому при t < T:

T0 (t ) = f (t-t)Q0 (т)d т = - ^^ - cos (T JJ ,

i * Q (t) nP0qnT\cos(?j-cos(c°nt)j

Tn (t )= — f sin (con (t-т)) Qn(\ d т =--1 1 T J ,-^;

«f ch f=* J ()2-n2 ) [a n

при t > T:

T0 (t ) = j(t -t)Q (T)d T = - 4P°T(<2 - b) r. an

Tn

ch !™h

Таким образом, общий вид потенциала следующий: t < Т:

+nP0Tj

n=1 ()2 -П2 )chí^^h

l[1 - í nt

cos I T )J+

ínt 1

V T) ü ch í nn 1

I — y I

™h I V a )

nn

a

t > Т:

M =- 4P0T( - b) +nP0T jj qn(cos К(T - t)) + cos Qy» ch fnny! cos (™xj.

an

((cOnT )2-n2 )ch (™h

Подстановка в (5) дает выражение для определения отклонения свободной поверхности г от состояния равновесия:

t < Т:

(

P0 • z = —- sin

1 sin г— 1-2Í£zbl sin Int J +

c - b ) V T ) an V T

ш qn í— sin í —|-ron sin

+nT j V T , V T ) n i n }) cos (™x ¿1 ()) -n2) 1 a

t > Т:

z =

nP0T j qnCn (sinК (T - tsin(cnt)) cos ínnx

g £

По приведенным формулам был проведен расчет формы свободной поверхности идеальной жидкости с помощью программы Mathcad. Плотность жидкости было взята 1000 кг/м3. Размеры сосуда: ширина - 1 м, глубина - 0.5 м. Максимальное значение давления в импульсе 100 Па. Длительность импульса составляет 0.1 с. Границы действия импульса 0.44 и 0.56 м от левой границы сосуда.

На рис. 3 приведена форма свободной поверхности жидкости в разные моменты времени: а - начало импульса, б - окончание действия импульса, в - формирование волны на поверхности, г - достижение волной границы сосуда. Расчеты показывают четкий момент отражения волн от стенок и последующее наложение колебаний с образованием сложной волновой поверхности.

u

2

n=1

z

0.01

0

-0.01

-0.02 z

0.01

0

-0.01

-0.02

0 0.2 0.4 0.6 0.8 х 0 0.2 0.4 0.6 0.8 х

Рис. 3. Форма свободной поверхности жидкости для импульсного внешнего давления Fig. 3. Shape of the free surface of liquid in the case of impulsive external pressure

Таким образом, предложенная методика дает возможность построить форму свободной поверхности идеальной жидкости при ее вынужденных движениях в ограниченном сосуде под действием переменного давления на поверхности жидкости.

ЛИТЕРАТУРА

1. Богоряд И.Б., Дружинин И.А., Дружинина Г.З., Либин Э.Е. Введение в динамику сосудов с жидкостью. Томск. Изд-во Том. ун-та, 1977. 143 с.

2. Моисеев Н.Н. О колебаниях тяжелой идеальной несжимаемой жидкости в сосуде // ДАН СССР. 1952. Т. 85. Вып. 5. С. 963-965.

3. Моисеев Н.Н., Черноусько Ф.Л. Задачи колебаний жидкости, подверженной силам поверхностного натяжения // Журнал вычислительной математики и математической физики. 1965. Т. 5. № 6. С. 1071-1095.

4. Петров А.А. Приближенный метод расчета собственных колебаний жидкости в сосудах произвольной формы и потенциалов Жуковского для этих сосудов // Журнал вычислительной математики и математической физики. 1963. Т. 3. № 5. С. 958-964.

5. Петров А.А., Попов Ю.П., Пухначев Ю.В. Вычисление собственных колебаний жидкости в неподвижных сосудах вариационным методом // Журнал вычислительной математики и математической физики. 1964. Т. 4. № 5. С. 880-895.

6. Балабух Л.И., Молчанов А.Г. Осесимметричные колебания сферической оболочки, частично заполненной жидкостью // Изв. АН СССР. Механика твердого тела. 1967. № 5. С. 22-26.

7. Пожалостин А.А. Свободные колебания жидкости в жестком круговом цилиндрическом сосуде // Изв. вузов. Авиационная техника. 1963. № 3. С. 25-32.

8. Сретенский Л.Н. Теория волновых движений жидкости. М.: Наука, 1977. 816 с.

9. Linton C.M., McIver P. Embedded trapped modes in water waves and acoustics // Wave Motion. 2007. V. 45. - No. 7-8. P. 940-951.

10. Motygin O.V. On trapping of surface water waves by cylindrical bodies in a channel // Wave Motion. 2008. V. 45. P. 940-951

11. Huang D., Guo W. and Li X. An analytical solution of fluid-structure coupling oscillation in one-dimensional ideal condition under small disturbance // J. Sound and Vibration. 2002. V. 255. No. 3. P. 610-614

12. Кочин Н.Е., Кибель И.А., Розе Н.В. Теоретическая гидромеханика. Ч. 1. М.: ГИФМЛ, 1963. 584 с.

13. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1977. 742 с.

Статья поступила 23.01.2017 г.

Merzlyakov A.V., Matyeva Z.O. ANALYTICAL SOLUTION OF THE PROBLEM OF SMALL FORCED OSCILLATIONS OF THE IDEAL FLUID. Tomsk State University Journal of Mathematics and Mechanics. 48. pp. 70-81

DOI 10.17223/19988621/48/7

The study aims to determine the shape of the free surface of an ideal liquid vibrating under variable external pressure in a rectangular vessel. The two-dimensional problem has been solved analytically. The motion of the ideal fluid has been simulated by the solution of Laplace's equation for the fluid velocity potential in the flow region. On the solid boundary, the impermeability conditions have been assigned; on the free surface, the Cauchy-Lagrange integral. The equation was solved using the variable separation method. The free surface condition has been transformed in the case of small oscillations and varying pressure on the free surface and it has been used to determine the velocity potential as a function of time. The implementation of this condition required solving the system of linear second-order differential equations. The same condition gives the formula for determining the shape of the free surface of liquid at any time instant in the form of deviation of free surface points from the equilibrium position. The method described in this paper has been applied in the cases with the external pressure varying harmonically and acting on the restricted part of the free surface during a limited time. The obtained results have been compared with the currently available solutions of similar problems.

Keywords: ideal fluid, velocity potential, Laplace's equation, method of separation of variables.

MERZLYAKOVAleksandr Vladimirovich (Candidate of Physics and Mathematics, Tomsk State University, Tomsk, Russian Federation) E-mail: amerz@mail.ru

MATYEVA Zarina Olegovna (Tomsk State University, Tomsk, Russian Federation) E-mail: amerz@sibmail.com

REFERENCES

1. Bogoryad I.B., Druzhinin I.A., Druzhinina G.Z., Libin E.E. (1977) Vvedenie v dinamiku sosudov s zhidkost'yu [Introduction to dynamics of vessels with a liquid]. Tomsk: TSU Publ. P.143.

2. Moiseev N.N. (1952) O kolebaniyakh tyazheloy ideal'noy neszhimaemoy zhidkosti v sosude [On oscillations of a heavy ideal incompressible fluid in a vessel]. DAN SSSR. 5(85). pp. 963-965.

3. Moiseev N.N., Chernous'ko F.L. (1965) Problems on the oscillation of a fluid subject to surface tension forces. USSR Comput. Math. & Math. Phys. 6(5). pp. 128-160. DOI: 10.1016/0041-5553(65)90102-3.

4. Petrov A.A. (1963) An approximate method for the calculation of characteristic oscillations of a liquid in vessels of arbitrary shape and the Zhukovskii potentials for these vessels. USSR Comput. Math. & Math. Phys. 5(3). pp. 1307-1316. DOI: 10.1016/0041-5553(63)90116-2.

5. Petrov A.A., Popov Yu.P. Pukhnachev Yu.V. (1964) Evaluation of the proper vibrations of a fluid in a fixed vessel by the variational method. USSR Comput. Math. & Math. Phys. 5(4). pp. 880-895. DOI: 10.1016/0041-5553(64)90144-2.

6. Balabukh L.I., Molchanov A.G. (1967) Osesimmetrichnye kolebaniya sfericheskoy obolochki, chastichno zapolnennoy zhidkost'yu [Axisymmetric oscillations of a spherical shell partially filled with a liquid]. Izvestiya AN SSSR. Mekhanika tverdogo tela. 5. pp. 22-26.

7. Pozhalostin A.A. (1963) Svobodnye kolebaniya zhidkosti v zhestkom krugovom tsilindricheskom sosude [Free oscillations of a liquid in a rigid circular cylindrical vessel]. Izvestiya visshikh uchebnykh zavedeniy. Aviatsionnaya tekhnika. 3. pp. 25-32.

8. Sretenskiy L.N. (1977) Teoriya volnovykh dvizheniy zhidkosti [Theory of wave motions in a fluid]. Moscow: Nauka.

9. Linton C.M., McIver P. (2007) Embedded trapped modes in water waves and acoustics. Wave Motion. 45. pp. 16-29.

10. Motygin O.V. (2007) On trapping of surface water waves by cylindrical bodies in a channel. Wave Motion. 45(7-8). pp. 940-951. DOI: 10.1016/j.wavemoti.2007.04.009.

11. Huang D., Guo W., Li X. (2002) An analytical solution of fluid-structure coupling oscillation in one-dimensional ideal condition under small disturbance. Journal of Sound and Vibration. 255(3). pp. 610-614. DOI:10.1006/jsvi.2002.5193.

12. Kochin N.E., Kibel' I.A., Roze N.V. (1963) Teoreticheskaya gidromekhanika. Chast' 1 [Theoretical Fluid Mechanics. Part 1]. Moscow: State Edition of Physical and Mathematical Literature.

13. Tikhonov A.N., Samarskiy A.A. (1977) Uravneniya matematicheskoy fiziki [Equations of mathematical physics]. Moscow: Nauka.