Possibility of existence an exact solution of Navier-Stokes equations on example of calculation the initial section in a circular tube



Cite item

Full Text

Abstract

The solution of the equations of motion of Newtonian fluids (Navier-Stokes equation) in vector form, containing a projection of speed and rotor speed, is considered. There was analyzed the axisymmetric laminar flow of a viscous incompressible fluid in the initial section of entrance to the circular tube, where the flow is described by only two velocity components.

Full Text

Движение частиц вязких сред описывается так называемыми уравнениями Навье-Стокса (далее - уравнения Стокса), которые выводятся согласно методу, впервые применённому Л.Эйлером, при котором рассматривается равновесие дифференциально малой частицы [1]. По существу уравнения Стокса представляют собой запись уравнения движения Ньютона при учёте внешних сил вязкости, которые пропорциональны градиентам скорости по граням частицы (для так называемых ньютоновских жидкостей). Уравнения Стокса представляют собой с математической точки зрения нелинейные дифференциальные уравнения в частных производных, которые не могут до настоящего времени быть решены в общем случае. Но данные уравнения могут быть рассмотрены с физической точки зрения как количественное описание механического движения сплошной среды, и с этой стороны можно указать на их неполноту, которая заключается в том, что они не описывают и не используют такую особенность движения как возможность движения вязких частиц по траекториям, имеющим кручение. Это игнорирование кручения является естественным следствием метода Эйлера. Как следует из [1, с. 419] уравнение Стокса при условии несжимаемости жидкой среды можно представить в следующем виде, в котором явно выделено выражение для полной энергии (суммы энергии давления, положения и кинетической энергии) жидкой частицы: . (1) Здесь вектор скорости имеет проекции - частное от давления в жидкости и её плотности, П - потенциал внешних сил, - кинематическая вязкость жидкости. Исследуем возможность упрощения вида уравнения (1), воспользовавшись тем, что для любого скаляра ротор его градиента равен нулю. Если для простоты рассматривать установившийся режим течения, то из (1) получим уравнение: . (2) Это векторное уравнение является векторной записью уравнения движения относительно скорости в частных производных третьего порядка. Существенно, что оно содержит только два члена. Один член квадратичный, а второй - первой степени от производных скорости, т.е. линейный. Поэтому из уравнения Стокса, записанного в виде (2) и не содержащего в такой форме давления и координат положения жидкой частицы, следует, что при условии сохранения картины течения (при условии сохранения линий токов и ламинарного режима течения и тем самым сохранения траекторий движения жидких частиц) при изменении скорости в m раз один член уравнения изменяется в m2 раз, а второй - только в m раз. Но при этом оба члена уравнения остаются равными. Такое положение возможно только, если каждый из членов уравнения тождественно равен нулю. В ином случае при изменении скорости течения согласно (2) должна изменяться и картина линий токов. Это утверждение можно проверить экспериментально, но желательно также иметь и теоретические значения скоростей течения в зависимости от величины расхода через трубу Q. В случае осесимметричного течения в круглой трубе вдоль координаты z проекция скорости т.е. . Поскольку одна из проекций скорости равна нулю, то запись уравнений движения для осесимметричного движения проще, чем пространственного. С учётом этого условие несжимаемости имеет вид [1, с. 324]: . (3) Выражение для ротора скорости V принимает вид: . (4) Здесь для краткости введено обозначение: . При таком обозначении векторные уравнения (2) сокращаются до одного (ненулевой проекции): . (5) Решение уравнения (5) будем искать в виде комбинации из степенного и тригонометрического рядов, представив аксиальную составляющую скорости в виде суммы степенного ряда и ряда Фурье по синусам, а именно в виде: . (6) Причём: Здесь - длина начального участка (участка Гагена, т.е. участка, на котором происходит перестроение начального профиля скоростей в пуазейлевый [2]). Ввиду громоздкости метода рядов степень полярного радиуса, как правило, выбираем не более 5, т.е. m ≤ 5. Начальное значение радиальной проекции скорости в силу её непрерывности считаем возможным представить в виде степенного ряда: . (7) Таким образом, возможное выражение для проекции скорости в развёрнутом виде, причём удовлетворяющее граничным условиям, может быть представлено как: (8) Следовательно, для определения радиальной проекции скорости (и тем самым и её осевой проекции) необходимо определить n значений постоянных коэффициентов mi, а также n2 коэффициентов bnm. Поскольку решение уравнений (1) является поиском решения задачи как математической, так и физической, то далее всюду будет иметь место обращение к физическому смыслу задачи. В этой связи отметим, что степень разложения в ряд должна начинаться с первой степени, так как при радиальная составляющая скорости по смыслу задачи равна нулю , но , и при этом на расстоянии x0 от входа в трубу радиальная составляющая становится равной нулю по всему сечению трубы: . Выражение (8) удовлетворяет перечисленным граничным условиям. Кроме того, ещё одним граничным условием является условие того, что на стенке трубы, т.е при r = r0 величина радиальной компоненты равняется нулю на всей длине начального участка, т.е. соответственно (8) получим n уравнений: (9) Аксиальную проекцию скорости определим из уравнения неразрывности (3) (при этом значение в силу её непрерывности и однозначности также считаем разлагаемым в степенной ряд) с учётом (6) и (7) следующим образом: (10) Здесь, как и в (7): Тогда: (11) Граничные условия для осевой проекции скорости состоят в том, что при входе в трубу распределение скоростей Uz (r; 0) имеет вид (10). По окончании начального участка длиной , когда радиальная проекция скорости полностью гасится для всех значений r, течение становится чисто Пуазейлевым и имеет место параболическое распределение скорости, т.е.: (12) Здесь - максимальная скорость течения на оси, т.е. - радиус трубы. Из (11) и (12) при получим: Но поскольку написанное выражение должно быть тождеством, то должны быть справедливы следующие соотношения между коэффициентами: (13) Далее, используя выражения для проекций скорости (8) и (11), вычислим проекцию ротора скорости согласно (4): (14) Или, группируя слагаемые по степеням полярного радиуса, получим вместо (14): (15) Полученное выражение громоздко, поэтому обозначим для краткости значения выражений в фигурных скобках при степенях полярного радиуса через А0(z), А1(z), А2(z) и т.д. Тогда вместо (15) получим более компактное выражение: (16) Теперь используем значения граничных условий. Из физического смысла задачи полагаем, что при входе в трубу, как и для всей массы жидкости вокруг трубы, значение вихря по всему входному сечению трубы равно нулю, т.е. . В конце начального участка для установившегося течения имеет место пуазейлевое распределение скоростей, в котором значение ротора согласно (11) равно: . Отсюда из (15) для указанных граничных условий получим соотношения: (17) Т.к. коэффициенты степенного ряда, сумма которого равна нулю, тоже должны быть равны нулю, то отсюда получим следующие соотношения: (18) Уравнения (18) связывают значения коэффициентов начального распределения скорости с коэффициентами, определяющими распределение скорости на начальном участке и т.д. Значение ротора скорости для пуазейлевого течения, т.е. в конце начального участка, известно и равно (2Ur/r02). Вычислим значение проекции ротора, даваемое рядом вида (15) в конце начального участка, где значение ротора становится равным пуазейлевому, т.е. при . (19) Из (19) с учётом следует: (20) Используя уравнение неразрывности (3), упростим выражение (5) с учётом (16), получим: (21) Отметим, что согласно (16): В соответствии с этим: . (22) Из этого выражения следует, что поскольку рассматриваются лишь конечные величины, то в (22) должно быть . Поэтому в соответствии с (15) получим: Отсюда: (23) Т.е. все коэффициенты при вторых степенях полярного радиуса в (8) равны нулю. Помимо этого отметим, что запись уравнения (5) движения с учётом того, что и (23) может быть представлена в виде: . (24) В этом уравнении все члены, кроме члена , имеют общий множитель r. Поэтому должно иметь место тождественное равенство . Отсюда в соответствии с (24) получим: Отсюда И, кроме того, используя (23) определим значения коэффициентов: (25) С учётом этого выражение для уравнения движения (5) можно записать в виде: (26) При указанном значении коэффициентов согласно (23) и (25) выражение для проекции скорости Ur становится гораздо менее громоздким, а именно, вместо (8) получим: (27) Таким образом, разложение радиальной проекции скорости содержит только нечётные степени полярного радиуса. Согласно (18) с учётом (23) получим . Кроме того: . Тогда для второй проекции скорости Uz согласно (11) получим: (28) Выражение (15) для проекции ротора скорости с учётом (23) и (25) также упрощается, а именно (поскольку ): Здесь: (29) Пользуясь формулами (27) и (28), с учётом того, что k1 = 0, получим следующие выражения для коэффициента при первой степени r при разложении в степенной ряд слагаемых уравнения (26): (30) Теперь отметим, что с точностью до первой степени разложения в степенной ряд по степеням радиуса в соответствии с (16) и с учётом (23) оказывается верным соотношение: (31) С учётом этого уравнение движения (26) принимает с точностью до первой степени полярного радиуса более простой вид: Или, после сокращения на r , следующее выражение: (32) Подставляя в (32) слагаемые, выраженные согласно (30) и с учётом (25), получим с точностью до первого порядка разложения в степенной ряд по радиусу следующее равенство: (33) Данное выражение должно быть тождеством, т.е. должно быть справедливым при всех значениях переменной z. Левая часть выражения (33) состоит из произведений суммы многочлена и тригонометрического ряда. Обозначим многочлены левого выражения через и , а суммы тригонометрических выражений через и . Соответственно, в левой части через и . Тогда (33) схематически можно записать как: (34) или: . Слагаемые и есть многочлены, а слагаемые и - это суммы, соответственно, косинусов и синусов кратных углов. Слагаемое представляет собой произведение тригонометрических функций и . Но произведение синусов и косинусов согласно известным тригонометрическим формулам может быть представлено в виде суммы синусов кратных углов и поэтому можно считать, что в левой части (34) имеется только сумма синусов кратных углов, а в правой части - сумма косинусов кратных углов. Тождественное равенство синусов и косинусов в выражении (33) возможно лишь в случае, если постоянные коэффициенты при этих функциях равны нулю. Поэтому тождественное равенство (34) может выполняться лишь в следующих рассматриваемых далее случаях. Во-первых, если все слагаемые равны нулю: . Это приводит к уравнениям , откуда: . Условие приводит к системе уравнений: (35) Из уравнения , получим: (36) Поскольку написанные равенства должны удовлетворяться тождественно, то . Но тогда . Однако тогда задача теряет физический смысл. Поэтому не может быть равен нулю и рассмотренный случай не может иметь место в действительности. Помимо рассмотренного случая тождественное равенство (34) может формально иметь место в случае, если , но при этом и - любые. Т.е.: . Отсюда . Из условия должны выполняться уравнения (35). Условие требует чтобы: Отсюда . Тогда: (37) Наконец из условия В3(z) = 0 получим систему уравнений (36). Используем теперь систему уравнений (9), которая после сокращения на r0 и с учётом (23) и (25) имеет вид однородной системы уравнений: (38) Подставляя в (38) значения и т.д., из (35) получим: (39) Выражая из этой системы значения и подставляя их в (36) с использованием (35), получим из (36) следующую систему уравнений: (40) Каждое уравнение системы состоит из двух сомножителей, хотя бы один из которых равен нулю. Поскольку значения корней уравнений в квадратных скобках оказываются для всех уравнений только мнимыми, то решением системы (40) может быть только случай: . Но в таком случае и остальные коэффициенты тригонометрического ряда, согласно (35) и (39), также равны нулю. Из (37) получим: , так что с учётом того, что m3 = 0 получим, что . Значение проекций скоростей тогда для рассматриваемого случая вместо (27) и (28) принимают вид: (41) При этих значениях проекций скорости уравнение неразрывности (3) выполняется. Выражение для ротора скорости в соответствии с (4) приобретает вид: . (42) Подставляя в уравнение движения (5) значения проекций скоростей (41) получим, что должно быть . Поэтому и . Следовательно, рассмотренный вариант не представляет интереса. Исследуем теперь возможность . Т.е. все коэффициенты тригонометрических рядов равны нулю. Следует отметить, что этот случай не учитывает изменение скорости по радиусу трубы, и поэтому он не удовлетворяет физическому смыслу задачи. Но формально, согласно (34), получим: . Или в развёрнутом виде: (43) Поскольку сумма коэффициентов при z3 в написанном тождестве должна быть равна нулю, то . Поэтому либо , либо . Аналогично коэффициент при также должен быть равен нулю, откуда . Отсюда . Сумма коэффициентов при также должна быть равна нулю . Сравнивая это выражение с предыдущим, получим: . Поэтому либо , либо . Наконец, сумма постоянных коэффициентов, не зависящих от переменной равенства (43), приводит к соотношению или . Таким образом, в рассматриваемом случае коэффициенты степенных рядов (как и тригонометрических) оказываются либо равны нулю, либо они не определены. Ещё один формальный случай выполнения равенства (33) требует условий: (44) Отсюда следует, что , , Исследуем эти условия. Условие означает справедливость системы уравнений (35). Условие эквивалентно в соответствии с (33) справедливости тождества: Из условия суммы коэффициентов при в написанном тождестве равной нулю получим: . Из условия суммы коэффициентов в этом равенстве при равной нулю получим то же самое соотношение. Равенство свободных членов приводит к уравнению: . Условие в развёрнутом виде приводит с учётом неизменной системы (38) к системе уравнений, аналогичных системе (36), а именно: (45) В этих уравнениях постоянная одна и та же, поэтому приходится считать, что коэффициенты ряда Фурье с индексом выше третьего равны нулю, т.е. что . Из (45) получим, исключая из двух уравнений и учитывая системы (35) и (38), одно уравнение: . Это алгебраическое уравнение не имеет действительных корней. Поэтому приходится считать, что все коэффициенты ряда Фурье равны нулю: . Кроме того, оказывается, что и все коэффициенты степенных рядов и равны нулю. Таким образом, согласно изложенному применение универсальных точных математических методов (степенных и тригонометрических рядов) не позволило добиться решения системы уравнений математической физики - (3) и (5). Это может быть объяснено либо тем, что эти уравнения противоречивы или не точны, или тем, что при решении поставленной задачи выбраны противоречивые граничные условия. Между тем, указанные уравнения при отмеченных граничных условиях допускают существование различных приближённых решений. Отметим, что первым на существование начального участка обратил внимание Гаген в своих опытах [2], а теоретически первым длину начального участка приближенно вычислил Буссинеск в 1891 г. [3]. Он полагал, что формирование ламинарного потока в трубе радиуса заканчивается, если скорость в конце начального участка на оси трубы достигает 0,99 от значения максимальной скорости U пуазейлевого течения. Другими словами он считал, что экспериментальная погрешность имеет порядок одного процента от максимальной скорости потока. При расчёте он использовал упрощённые уравнения Навье-Стокса и допущение о равномерном профиле поля скоростей при входе в трубу. Согласно Буссинеску длина начального участка равна: , где число Рейнольдса (46) Помимо Буссинеска длину начального участка в 1931 г. вычислил Шиллер Л. [2], который использовал для этого приближённую теорию ламинарного пограничного слоя Прандтля. Шиллер получил для длины начального участка выражение: . (47) Наконец, Слёзкин Н.А. [4] в 1955 г., используя метод операционного исчисления и также допуская погрешность в 0,01 от максимальной скорости течения по оси трубы, получил следующее значение для длины начального участка: . (48) Таким образом, известно несколько разных приближённых теоретических формул, согласно которым результаты различаются более, чем в два раза. Как показано выше, попытка автора использовать универсальные точные аналитические методы (методы степенных и гармонических рядов) для получения точного решения (1) не позволила получить точное решение. Из приближённых формул (46) - (48), полученных совершенно разными методами, следует, что длина начального участка пропорциональна скорости течения в трубе. Это означает, что линии тока на начальном участке меняются при изменении скорости и возможно, что переход от ламинарного течения в турбулентное в опытах Рейнольдса обусловлен тем, что толщина пристеночного слоя растёт с ростом средней скорости течения и достигает оси трубы, вследствие чего и происходит переход от ламинарного течения в турбулентное (как иногда говорят, наступает потеря устойчивости течения), рисунок 1. Возможно также, что толщина пристеночного слоя имеет волнистую, колеблющуюся форму (рисунок 1, справа, граница пристеночного слоя 3). а) возможный характер изменения скорости и образования пристеночного слоя на протяжении начального участка б) начало турбулизации потока при приближении пристеночных слоём толщиной (r. z) друг к другу Рисунок 1 На рисунке 1а показан характер перестроения скоростей на начальном участке при входе в круглую трубу радиусом при ламинарном режиме течения: r - расстояние от оси трубы 0z, I - сечение трубы с максимальным сужением потока, II - сечение трубы с участком замедления скорости течения, III - сечение, в котором профиль скоростей становится квадратичным, - угол наклона траектории вблизи кромки трубы при входе, - толщина пристеночного слоя. На рисунке 1б: 1, 2 и 3 - различные возможные формы линии тока, касающейся входной острой кромки трубы. Xmax - наибольшая возможная длина начального участка, в котором сохраняется ламинарный режим течения. Это предположение основано на таком опытном факте, что граница ламинарного потока и неподвижной жидкости всегда приобретает волнистую форму [5], как это показано на рисунке 2 для течения в прямоугольной полости (каверне). Амплитуда волн пропорциональна скорости потока и при достаточно высокой скорости вершины волн размываются в потоке, а при дальнейшем увеличении скорости вся граница приобретает беспорядочный, турбулентный характер движения. а) б) Рисунок 2 На рисунке 2а характер линий токов (зарисовка с видеоклипа) при обтекании прозрачной жидкости (воды) со скоростью U прямоугольной полости шириной 5 сантиметров, заполненной подкрашенной жидкостью (тоже вода) [5]. Уже при самых малых скоростях течения граница между неподвижной подкрашенной жидкостью и потоком с неподкрашенной жидкостью начинает колебаться с образованием волн (как на поверхности воды волн от ветра), амплитуда которых растёт с ростом скорости потока. Амплитуда волн нарастает по ходу течения. Перемешивания жидкости в полости (каверне) практически не наблюдается. По мере роста скорости течения основного потока и роста амплитуды волн начинает происходить перемешивание подкрашенной жидкости полости и прозрачной жидкости основного потока вплоть до полного вымывания подкрашенной жидкости из полости. На рисунке 2б схема переходного от ламинарного к турбулентному режиму течения в опыте Рейнольдса с подкрашиванием центральной струйки втекающей в трубу воды: 1 - сосуд с подкрашенной водой, 2 - резервуар с отстоявшейся водой, 3 - стеклянная труба с плавным входом, 4 - центральная подкрашенная струйка в переходном от ламинарного к турбулентному режиму течения. Колебания центральной струйки происходят как в вертикальной, так и в горизонтальной плоскостях. Кроме того, в опытах уже самого Рейнольдса было установлено, что переход от слоистого, ламинарного течения в беспорядочное, турбулентное при возрастании скорости происходит не внезапно, а через промежуточный смешанный режим, при котором наблюдается попеременное чередование ламинарного и турбулентного течения (рисунок 2). Отметим также, что переход от ламинарного течения в турбулентное очень сильно зависит от времени предварительного отстоя жидкости в основной ёмкости, из которой и происходит истечение в круглую трубу. Поэтому критическое число Рейнольдса в трубах, как свидетельствуют опыты на воде, может меняться от 2300 до почти 50 000. Соответственно, и длина начального участка может меняться при одной и той же скорости течения в десятки раз. Фактически используется нижняя граница начального участка и не ясно, какова его верхняя граница. Отметим также, что при переходе от ламинарного режима течения к турбулентному в опыте Рейнольдса наблюдается синусообразная форма первоначально центральной струйки тока, и поскольку изогнутость траектории течения должна наблюдаться не только в вертикальной, но и в горизонтальной плоскостях, то, возможно, реальная траектория имеет спиралевидную форму. При теоретических расчётах все указанные особенности течения на начальном участке не учитываются, а имеющийся экспериментальный материал не указывает на то, как эти особенности учесть. Помимо изложенного имеются и другие причины для сомнений в том, что уравнения Навье-Стокса достаточно адекватно описывают поведение вязкой жидкости. Это, например, то обстоятельство, что при их выводе не учитывается возможность кручения линий тока. И эта особенность вообще характерна для всей современной теории поля [2]. В результате гидромеханика не может объяснить такие наблюдаемые явления, как сравнительно давно известные вихри Гёртлера (рисунок 3) или винтовое движение шаров в жидкости, наблюдаемое в природе в виде винтообразного подъёма воздушных пузырьков в воде [6, 7]. а) б) Рисунок 3 На рисунке 3а - схема вращения шариков, соединённых в виде «гантели» 6 с угловой скоростью ω1 в воде, буксируемой с постоянной скоростью V нитью 3, прикреплённой к поводку 5. Зарисовка с видеоклипа [6]. На рисунке 3б - вихри Гёртлера при обтекании цилиндра равномерным потоком воздуха. Зарисовка с фотографии. Характер линий тока обнаруживается по отклонению шелковинок при продувке цилиндра в аэродинамической трубе при (опыты ЦАГИ) [7]. Положение не меняет и использование для решения уравнений Навье-Стокса численных методов. В настоящее время достаточно широко известны и применяются многоцелевые пакеты программ типа NASTRAN, ANSIS, KATIA и др. В этих программах для динамического анализа механических систем используется метод конечных элементов. В них, однако, не учитывается возможность движения частиц по траекториям, имеющим кручение. Логично пытаться применить эти пакеты математического программного обеспечения для расчёта начального участка течения в круглой трубе и сравнить с результатами, полученными классической математикой, а также с опытными данными. В любом случае целесообразно было бы сравнить результаты компьютерного счёта при использовании разных пакетов программ, и, как максимум, получить опытные данные о поведении жидкостей (и газов) на входном участке круглой трубы при входе среды из неограниченного пространства, в частности, определить характер траекторий течения на входном участке. То же самое можно отнести и к пространственным движениям жидкости при обтекании цилиндра и шара [2, 6, 7]. По мнению автора уже известные опытные данные позволяют сделать вывод о том, что вращательное движение есть изначальное свойство сплошных сред и оно всегда проявляется при достаточно больших скоростях движения. Можно отметить, не упомянув о других подобных примерах, резкое отличие картины течения для так называемой «тестовой» задачи по течению жидкости в каверне в случае использования численных методов и результатов опытов. Принципиальное различие меду этими двумя случаями сохраняется вне зависимости от быстродействия и объёма памяти используемых компьютеров. Ситуацию, сложившуюся к настоящему времени, возможно, удастся изменить уточнением записи уравнений движения, изменив их современный вид (в форме Навье-Стокса). Это можно сделать только на основании гораздо более обширного экспериментального материала, чем тот, который послужил для вывода уравнения движения вязких жидкостей более 150 лет назад сначала Навье (1821 год), а потом Стоксу (1845 год) [1]. По мнению автора именно опыты по наблюдению за потоком на начальном участке входа в цилиндрическую трубу способны дать соответствующую «базу» данных. Можно отметить, что результаты предварительных опытов по наблюдению характера движения на входном участке круглой трубы (опыты на воздухе, диаметр трубы 4,0 см, ) показывают, что траектории движения на этом участке имеют спиральный характер, в опытах Рейнольдса можно различить свидетельство наличия на начальном участке вращательного движения частиц жидкости, поскольку колебания подкрашенной центральной струйки воды в его опытах не должны иметь предпочтения к вертикальной плоскости, а должны происходить во всех иных плоскостях [8, 9].
×

About the authors

V. G. Vyskrebtsov

Moscow State University of Mechanical Engineering (MAMI)

Email: tm@mami.ru
Ph.D.

References

  1. Лоцянский Л.Г. Механика жидкости и газа. Изд-во «Наука», М, 1973.
  2. Пискунов Н.С. Дифференциальное и интегральное исчисления для ВТУЗов. Том 2, М, Издательство «Наука», М, 1976.
  3. Шиллер Л. Движение жидкости в трубах. М - Л, 1936.
  4. Тарг С.М. Основные задачи теории ламинарных течений, Гостехиздат, 1951
  5. Слёзкин Н.А. Динамика вязкой сжимаемой жидкости. Изд-во технико-теоретической литературы, М, 1955.
  6. Выскребцов В.Г. Установившееся течение вязкой жидкости при входе в трубу. Известия МГТУ «МАМИ», М, №1(15), серия 3. Естественные науки, 2013. Том 3.
  7. Короткин А.И. О трёхмерном характере обтекания кругового цилиндра. «Учёные записки ЦАГИ», М, 1973, том VI, №4.
  8. Ван-Дайк М. Альбом течений жидкости и газа. М, Наука, 1986.
  9. Тарг С.М. Основные задачи теории ламинарных течений. М, Гостехиздат, 1951.

Supplementary files

Supplementary Files
Action
1. JATS XML

Copyright (c) 2015 Vyskrebtsov V.G.

Creative Commons License
This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License.

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies