CC BY
48
2014 Математика и механика № 5(31)
УДК 532.516.5
В.А. Якутенок, М.А. Пономарева, А.Е. Кузнецова
МОДЕЛИРОВАНИЕ ОСЕСИММЕТРИЧНЫХ ТЕЧЕНИЙ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ НЕПРЯМЫМ МЕТОДОМ ГРАНИЧНЫХ ЭЛЕМЕНТОВ1
Приводится вывод выражений для компонент фундаментальных тензоров скоростей и усилий для уравнений Стокса в цилиндрической системе координат для осесимметричного случая. На основе полученных выражений выписаны граничные интегральные уравнения в соответствии с представлениями теории потенциала. Предложен простой вариант получения приближенного решения, являющийся реализацией основных положений непрямого метода граничных элементов. Достоверность результатов подтверждается решением тестовой задачи. Разработанный метод решения краевых задач для уравнений Стокса может использоваться для моделирования ползущих течений вязкой жидкости со свободной поверхностью.
Ключевые слова: вязкая жидкость, метод граничных элементов, фундаментальные решения, осесимметричное течение.
Метод граничных элементов (МГЭ) является эффективным средством решения различных задач механики сплошных сред, в частности задач динамики жидкости. Существуют различные его варианты, причем главная разница в используемых подходах обусловлена выбором прямой или непрямой формулировки метода [1, 2]. Соотношения непрямого МГЭ для плоских и пространственных течений вязкой жидкости и примеры расчетов приведены в [3, 4]. Данный подход позволил провести ряд исследований практически важных двумерных течений, характерной особенностью которых является наличие свободной поверхности [5-8]. Возможен учет нелинейного поведения жидкости [9].
В некоторых случаях течение имеет осевую симметрию. Например, при заполнении цилиндрических пресс-форм или растекании капли жидкости по горизонтальной поверхности. Тогда становится возможным существенно упростить трехмерный анализ, сведя его практически к решению задачи в одной плоскости, проходящей через ось симметрии. Для этого необходимо использовать вместо декартовой цилиндрическую систему координат и соответствующие ей выражения для фундаментальных сингулярных тензоров скоростей и усилий. Последний необходим при моделировании течений со свободной поверхностью. Этому вопросу посвящено значительное количество работ, в основном касающихся теории упругости, например [10, 11]. В [12] приводятся соотношения для компонент фундаментального тензора скоростей (в том числе проинтегрированных по угловой координате) для уравнений Стокса при использовании прямого МГЭ. В итоге получаются довольно громоздкие формулы, содержащие различные комбинации эллиптических интегралов первого и второго рода. Как следует из работ по теории упругости, для усилий получаются еще более сложные выражения [13].
1 Исследование выполнено при финансовой поддержке гранта Президента РФ (МК-3687.2014.1) и РФФИ в рамках научного проекта № 14-08-31579 мол_а.
В настоящей работе приводится полный вывод выражений для компонент фундаментальных тензоров скоростей и усилий в цилиндрической системе координат, необходимых для решения осесимметричных задач динамики вязкой жидкости, и формулируется гранично-интегральная постановка задачи в непрямом представлении. Проверка полученных соотношений проведена с помощью решения смешанной краевой задачи (с граничными условиями первого и второго рода), имеющей аналитическое решение.
Фундаментальные сингулярные решения уравнений Стокса в осесимметричном случае
Уравнения Стокса и неразрывности для случая единичной сосредоточенной силы, действующей в точке 4 в направлении к декартовой системы координат можно записать в безразмерной форме в виде
м.
дх}
= д1к5(хк -%к), (1)
= 0, (2)
duk
дхг
где <jk = -pk5j + 2ék - компоненты тензора напряжений, pk - давление, 5k - сим-
вол Кронекера, ék =1
дх: дхг
i г
- компоненты тензора скоростей деформаций,
5(хк - - дельта-функция Дирака, хк, ^ (к = 1, 2, 3) - декартовые координаты точек х и 4.
Все переменные являются безразмерными, причем используются масштабы: длины - Я, скорости - и, давления - ди/Я, где д - коэффициент динамической вязкости, а Я, и - некоторые характерные величины задачи. При указанном выборе масштаба давления уравнения (1), (2) не содержат безразмерных комплексов.
Решение задачи (1), (2) можно получить с помощью преобразования Фурье, что дает [14]
ч- 1 (дгк + У>Ук 1 „к,% ч_ Ук
uk (х,5) = + ^f- I, pk (X, 5) = , (3)
р р3 ) 4пр3
где Уг = х, - р = [(хг - ^(х, - £;)]1/2.
При исследовании осесимметричных течений естественно использовать цилиндрическую систему координат r, 9, r (рис. 1).
Расположим точку наблюдения х в плоскости х10х3, в которой и будет определяться решение. Необходимо иметь выражения для компонент фундаментального тензора скоростей:
< (х, 5), < (х, 5), иГ (х, 5), u\(х, 5).
Здесь, по-прежнему, нижний индекс обозначает компоненту вектора скорости в точке х, соответствующую направлению действия силы в точке 5, которое обозначается верхним индексом. Пусть в точке 5 действует единичная сосредоточенная сила er(5)=cos 9^ e1+ sin 9^ e2 (рис. 2).
0 rx
Рис. 1. Системы координат
Рис. 2. Орты систем координат
Очевидно, что er(x) = e1, ez(x) = e3, так как 9x = 0. Тогда в точке x, согласно принципу суперпозиции, будем иметь
urr (x, 4) = cos 9¡ • u1 (x, 4) + sin 9¡ • «2 (x, 4), (4)
urz (x, 4) = cos 9¡ • «3 (x, 4) + sin 9¡ • «32 (x, 4). (5)
Аналогично для силы ez(4):
«Z (x, 4) = «i3 (x, 4), «z (x, 4) = «3 (x, 4). (6)
Таким образом, все необходимые компоненты выражаются через решение (3) для общего трехмерного случая в декартовой системе координат. Все остальные компоненты фундаментального тензора не нужны. Например, хотя
«е (x, 4) = cos 9¡ • «2 (x, 4) + sin 9¡ • «2 (x, 4) ^ 0, но при дальнейшем интегрировании по углу в силу осесимметричности получим ноль. Из рис. 1 можно усмотреть, что
р¡Т+д2=VcZX-zj^'+^+rF-^^Coser,
>"1 = rx - r¡j cos9?, >-2 = - r5 sin 94, >-3 = zx - z¡j. Тогда, используя формулы (3) - (6), получим выражения для компонент фундаментального тензора скоростей в цилиндрической системе координат:
< (x, 4) = -
1
cos 9¡ cos 9¡ (rx + r¡) - rxr¡¡ (cos 9¡ +1)
8n
«1 (x, 4) = -«Z (x, 4) = -
р р
1 (Zx - )(rx cos 9¡- r¡ )
8п р3 '
1 (Zx - )(rx - r¡ cos 9¡ )
(7)
8n
р
1 + (Zx - Z¡ )
2 Л
< (X, =
При решении задач со смешанными краевыми условиями (например, при наличии свободной поверхности) на части границы области задаются компоненты вектора усилий:
(8)
где п - компоненты внешней нормали к границе в точке х. Выражения (8), (3) позволяют получить формулу для компонент фундаментального тензора усилий для уравнений Стокса в декартовой системе координат:
^ (х, 4) = -
4пр5
• У,У,Укп,(х).
(9)
Поступая совершенно аналогично предыдущему, получим компоненты этого тензора в цилиндрической системе координат:
3
1ГГ (Л^= —— (гх -г сое0,)( сое0, -г)[з (гх -г% сое0,) + пг ( -г,)], 4пр
3
£(х,= —— ( -1,) соБ-пс)\_пг ( -г% сов) + п ( -4пр
3
К (Л = —— (( - 2 с) - гссо5 ес) [п ( - гс со!5 ес)+ п ( - гс) 4пр
3 2
К (Л = —— (( - 2С ) [пг (Гх - ГС ^ ес) + П2 (х - 2С ) .
4пр
Соотношения (7), (9) позволяют перейти к гранично-интегральной постановке краевой задачи (1), (2) в осесимметричном случае.
Гранично-интегральная формулировка задачи и метод решения
Течение рассматривается в осесимметричной области О с поверхностью S, получаемой вращением образующей Г вокруг оси 2 (рис. 3).
Рис. 3. Область течения
Само течение также осесимметричное, т.е. скорость движения жидкости зависит только от г, 2 и не зависит от угловой координаты 9. В соответствии со стандартными представлениями теории потенциала будем считать, что по поверхности £ области О распределены фиктивные силы с плотностью на единицу площади £ равной ф(4). Таким образом, используется представление, аналогичное тому, что приводит к понятию потенциала простого слоя. Причем вектор ф(4) в силу
осесимметричности имеет вид ф(%) = ф,£г(%) + фгег(%). Тогда, используя принцип суперпозиции и фундаментальные решения, соответствующие единичным сосредоточенным силам ег(%), ег(%), можно записать для точек Хо, принадлежащих образующей Г:
иг (Хо) = { [фг (5К (Хо, %) + фг (%К (Хо, %)] йГ(%), (10)
г
и2 (Хо) = { [фг (% к (Хо, %) + ф г (%К (Хо , %) ] й Г(%), (11)
г
(г (Хо) = { [фг (%)?; (Хо, %) + фг (%)?; (Хо, %)] й Г(%), (12)
г
(Хо) = { [фг (%)?; (Хо , %) + ф г (%)?/ (Хо , %) ] й Г(%), (13)
где
2п 2п К(Хо, %) = г%\ игг (Хо,%)й, йгг (Хо,%) = г%\ ыгг (Хо,%)й,
о о
2п 2п
К(Хо, %) = г%\ К(Хо, %)й, К(Хо, %) = г%\ и1(Хо, %)й,
о о
2п 2п
(14)
Гг (Хо, %) = ^ (Хо, %)й9^, Гг (Хо, %) = ?г (Хо, %)й9^,
о о
2п 2п
II (Хо, %) = г: (Хо, %)й9^, ?/ (Хо, %) = £ (Хо, %)й9^,
о о
а также принято во внимание, что элемент поверхности равен й5(%) = г^ йГ(%). В уравнениях (1о) - (13) подразумевается, что точка % принадлежит образующей Г.
Краевые условия в общем случае заключаются в том, что на части образующей Г задано значение вектора скорости (твердая стенка либо входная или выходная граница), а на другой части вектор усилий (свободная поверхность или граница раздела). Таким образом, левая часть в уравнениях (1о) - (13) известна и задача заключается в определении неизвестной функции плотности ф(%). После этого по формулам, аналогичным (1о) - (13), можно найти значения вектора скорости и других характеристик течения в любой внутренней точке Х.
Основные уравнения (1о) - (13) можно записать в компактной форме:
й (Хо) = { й/ (Хо, п)ф; (П)й Г( п), (15)
г
^(Хо) = {I/ (Хо,п)ф] (П)йГ(п), (16)
г
где индексы /, ] = 1, 2 соответствуют координатам Пъ"П2, т.е., фактически г и г (рис. 4).
2п
Тогда й/ (Хо, п) = П | К (Хо, п)й9^, (17)
о
2п
I/ (Хо,п) = П1 | Ч (Хо,п)й9л, (18)
о
причем точки п^Г.
2,П 2
Уравнения (15), (16) практически тождественны уравнениям для плоского случая, поэтому для их решения можно использовать метод, приведенный в [3]. Главное отличие заключается в том, что соответствующие интегралы крайне затруднительно вычислить аналитически, поэтому используются стандартные численные квадратуры. Отметим, что подынтегральное выражение в (15) имеет слабую особенность, эквивалентную 1п р. Интеграл (16) существует в смысле главного значения по Коши, поэтому при совпадении точек х0 и ц следует использовать анализ, проведенный в [14]. Это означает, что диагональные элементы матрицы, получаемые после
дискретизации задачи (рис. 4) равны 1ф/(х0).
2
Вычислительные свойства получаемой в итоге системы линейных алгебраических уравнений аналогичны плоскому случаю. Поэтому при получении решения были использованы технологии параллельных вычислений хорошо зарекомендовавшие себя ранее [15].
г ,Л1
Рис. 4. Гранично-элементная дискретизация
Решение тестовой задачи
Краевые условия тестовой задачи представлены на рис. 5, причем в качестве масштаба длины Я выбран радиус цилиндрической трубы, а масштаба скорости и - средняя скорость течения. Очевидно, что решение уравнений Стокса при таких условиях соответствует течению Пуазейля. Целью является проверка полученных выше соотношений НМГЭ, включая формулы для фундаментальных сингулярных решений (7), (9) путем сравнения с известным аналитическим решением. Отметим, что представленная краевая задача имеет все основные типы граничных условий, которые встречаются в практических приложениях, а именно: входная граница, твердая стенка и выходная граница с заданными значениями компонент вектора усилий (аналог свободной поверхности, либо границы раздела). Одновременно проверяется возможность использования простых численных квадратур (типа формул Гаусса) для вычисления интегралов со слабыми особенностями.
Н = -4П1, Ь = 0
и 1 = 0, и2 = 2(1-п12) Рис. 5. Тестовая краевая задача
П1
Результаты вычислений представлены на рис. 6 и 7. Аппроксимационная сходимость оценивалась в норме Ь2, т.е. отклонение профиля скорости, рассчитанного на выходной границе от аналитического, вычислялось по формуле
Е = N£>/(иР - и* )2 , л /=1
где игР - рассчитанные значения аксиальной компоненты вектора скорости на выходной границе, и" - аналитические значения, / - номер узла (середины) /-го
граничного элемента, N - количество элементов на выходной границе. Общее число граничных элементов равно 3N. Полученные результаты подтверждают работоспособность изложенного в работе варианта НМГЭ.
Рис. 6. Профили аксиальной скорости Рис. 7. Аппроксимационная сходимость
на выходной границе для различных N в норме Ь2
Заключение
Получены выражения для компонент фундаментальных тензоров скоростей и усилий, необходимых для моделирования ползущих течений вязкой жидкости в осесимметричном случае. Представлен упрощенный вариант непрямого метода граничных элементов, основанный на численном вычислении всех необходимых интегралов, включая интегралы со слабыми особенностями. Достоверность приведенных в работе соотношений подтверждена решением тестовой краевой задачи.
ЛИТЕРАТУРА
1. Бенерджи П., Баттерфилд Р. Методы граничных элементов в прикладных науках: пер. с англ. М.: Мир, 1984.
2. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов: пер. с англ. М.: Мир, 1987.
3. Якутенок В.А. Численное моделирование медленных течений вязкой жидкости со свободной поверхностью методом граничных элементов // Математическое моделирование. 1992. Т. 4. № 10. С. 62-70.
4. Якутенок В.А. Численное решение трехмерных задач о ползущем течении вязкой жидкости со свободной поверхностью методом граничных элементов // Математическое моделирование. 1999. Т. 11. № 10. С. 92-99.
5. Штоколова М.Н., Шрагер Г.Р., Якутенок В.А. Формирование свободной поверхности объема вязкой жидкости внутри вращающегося горизонтального цилиндра // Изв. РАН. Механика жидкости и газа. 2009. № 2. С. 179-185.
6. Пономарева М.А., Шрагер Г.Р., Якутенок В.А. Устойчивость плоской струи высоковязкой жидкости, натекающей на горизонтальную твердую плоскость // Изв. РАН. Механика жидкости и газа. 2011. № 1. С. 53-61.
7. Новошинцев А.В., Шрагер Г.Р., Якутенок В.А. и др. Численное моделирование истечения вязкой жидкости из объемного смесителя // Теоретические основы химической технологии. 2006. Т.40. № 6. С. 668-674.
8. Пономарева М.А., Шрагер Г.Р., Якутенок В.А. Использование уравнения Дюпре - Юнга для решения задачи о растекании жидкости при ограниченном смачивании // Вестник Томского государственного университета. Математика и механика. 2008. № 1. С. 90-96.
9. Якутенок В.А., Штоколова М.Н. Численное моделирование плоских течений неньютоновской жидкости со свободной поверхностью непрямым методом граничных элементов // Вычислительные технологии. 2006. Т. 11. № 5. С. 106-118.
10. Yun B.I., Ang W.T. A dual-reciprocity boundary element method for axisymmetric ther-moelastostatic analysis of nonhomogeneous materials // Engineering Analysis with Boundary Elements. 2012. № 36. С. 1776-1786.
11. Karageorghis A., Fairweather G. The method of fundamental solutions for axisymmetric elastidty problems // Computational Mechanics. 2000. № 25. С. 524-532.
12. Pozrikidis C. Boundary integral and singularity methods for linearized viscous flow. Cambridge University Press, 1992.
13. ParkK.H., Benerjee P.K. A new BEM formulation for transient axisymmetric poroelasticity via particular integrals // International Journal of Solids and Structures. 2007. № 44. С. 7276-7290.
14. Ладыженская О.А. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970.
15. Пономарева М.А. Анализ эффективности использования средств распределенных вычислений для систем с общей памятью при моделировании течений вязкой жидкости методом граничных элементов // Седьмая Сибирская конференция по параллельным и высокопроизводительным вычислениям / под ред. проф. А.В. Старченко. Томск: Изд-во Том. ун-та, 2014. С. 137-144.
Статья поступила 01.07.2014 г.
Yakutenok V.A., Ponomareva M.A., Kuznetsova A.E. MODELING OF AXISYMMETRIC VISCOUS FLOWS OF INCOMPRESSIBLE FLUID BY THE BOUNDARY ELEMENT METHOD
In this paper, the derivation of the fundamental velocity and traction tensor components is presented for Stokes equations in the cylindrical coordinate system for axisymmetric viscous flows. The derivation is based on a three-dimensional singular solution in Cartesian coordinates. Using the presented formulas, the boundary integral equations are constructed in accordance with the potential theory conception. A simple numerical solution algorithm that is an implementation of principles of the indirect boundary element method is proposed. Reliability of the results is verified by solving a test problem the role of which is played by the problem about a flow in a cylindrical tube (the Poiseuille flow) with mixed boundary conditions. The developed method of solving boundary problems for Stokes equations can be used for modeling creeping flows of a viscous fluid with a free surface.
Keywords: viscous fluid, boundary element method, fundamental solutions, axisymmetric flow
YAKUTENOK Vladimir Albertovich (Doctor of Physics and Mathematics, Prof., Tomsk State University, Tomsk, Russian Federation) E-mail: [email protected]
PONOMAREVA Maria Andreevna (Candidate of Physics and Mathematics, Assoc. Prof., Tomsk State University, Tomsk, Russian Federation) E-mail: [email protected]
KUZNETSOVA Alina Evgenievna (Tomsk State University, Tomsk, Russian Federation) E-mail: [email protected]
REFERENCES
1. Banerjee P.K., Butterfield R. Boundary Element Methods in Engineering Science. McGraf-Hill book company, UK, 1981.
2. Brebbia C.A. Telles J.C.F., Wrobel L.C. Boundary Element Techniques. Theory and applications in engineering. Berlin, Sprinqer-Verlag, 1984.
3. Yakutenok V.A. Chislennoe modelirovanie medlennykh techeniy vyazkoy zhidkosti so svobodnoy poverkhnost'yu metodom granichnykh elementov. Matematicheskoe modelirovanie, 1992, vol. 4, no. 10, pp. 62-70. (in Russian)
4. Yakutenok V.A. Chislennoe reshenie trekhmernykh zadach o polzushchem techenii vyazkoy zhidkosti so svobodnoy poverkhnost'yu metodom granichnykh elementov. Matematicheskoe modelirovanie, 1999, vol. 11, no. 10, pp. 92-99. (in Russian)
5. Shrager G.R., Shtokolova M.N., Yakutenok V.A. Formation of the free surface of a viscous fluid volume inside a rotating horizontal cylinder. Fluid Dynamics, 2009, vol. 44(2). pp. 322327.
6. Ponomareva M.A., Shrager G.R., Yakutenok V.A. Stability of a Plane Jet of a Highly Viscous Fluid Impinging on a Horizontal Solid Wall. Fluid Dynamics, 2011, vol. 46(1), pp. 44-50.
7. Novoshintsev A.V., Shrager G.R., Yakutenok V.A. et al. Numerical modeling of the outflow of a viscous liquid from a bulk mixer. Theoretical Foundations of Chemical Engineering, 2006, vol. 40(6). pp. 626-632.
8. Ponomareva M.A., Shrager G.R., Yakutenok V.A. Ispol'zovanie uravneniya Dyupre - Yunga dlya resheniya zadachi o rastekanii zhidkosti pri ogranichennom smachivanii. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika, 2008, no. 1, pp. 90-96. (in Russian)
9. Yakutenok V.A., Shtokolova M.N. Chislennoe modelirovanie ploskikh techeniy nen'yutonovskoy zhidkosti so svobodnoy poverkhnost'yu nepryamym metodom granichnykh elementov. Vychislitel'nye tekhnologii, 2006, vol. 11, no. 5, pp. 106-118. (in Russian)
10. Yun B.I., Ang W.T. A dual-reciprocity boundary element method for axisymmetric thermoelastostatic analysis of nonhomogeneous materials. Engineering Analysis with Boundary Elements, 2012, no. 36, pp. 1776-1786.
11. Karageorghis A., Fairweather G. The method of fundamental solutions for axisymmetric elastisity problems. Computational Mechanics, 2000, no. 25, pp. 524-532.
12. Pozrikidis C. Boundary integral and singularity methods for linearized viscous flow. Cambridge University Press, 1992.
13. Park K.H., Benerjee P.K. A new BEM formulation for transient axisymmetric poroelasticity via particular integrals. International Journal of Solids and Structures, 2007, no. 44, pp. 7276-7290.
14. Ladyzhenskaya O.A. The Mathematical Theory of Viscous Incompressible Flow. 2nd ed. New York, Gordon & Breach, 1969.
15. Ponomareva M.A. Analiz effektivnosti ispol'zovaniya sredstv raspredelennykh vychisleniy dlya sistem s obshchey pamyat'yu pri modelirovanii techeniy vyazkoy zhidkosti metodom granichnykh elementov. Sed'maya Sibirskaya konferentsiya po parallel'nym i vysokoproizvoditel'nym vychisleniy am. Tomsk, Izd-vo Tom. un-ta, 2014, pp. 137-144. (in Russian)
