Solving axisymmetric potential problems using the indirect boundary element method.pdf Метод граничных элементов (МГЭ) является эффективным средством решения задач механики сплошных сред, в частности, сводящихся к уравнению Лапласа с различными краевыми условиями. Это особенно актуально для исследования потенциальных течений жидкости со свободной поверхностью. Во многих случаях эта проблема имеет осесимметричный характер. При этом фундаментальные решения уравнения Лапласа выражаются через эллиптические интегралы первого и второго рода, что значительно усложняет их интегрирование, особенно в сингулярном случае. В [1] для вычисления сингулярных интегралов предлагается использовать специальную квадратурную формулу для интегралов с логарифмическими особенностями. В [2, 3] в целях облегчения вычисления интегралов описывается переход от эллиптических интегралов к функции Лежандра второго рода. В [4] приводится модификация МГЭ с использованием мультипольного разложения основных соотношений. Альтернативой непрямому методу граничных элементов выступает метод фундаментальных решений. При его использовании сингулярные интегралы не возникают, однако сложность заключается в выборе точек приложения фиктивных нагрузок. Этот метод для решения осесимметричных задач рассматривается в [5, 6]. В настоящей работе вычисление сингулярных интегралов производится с помощью аппроксимации эллиптических интегралов многочленами [7], разложения оставшейся части подынтегральной функции в ряд Тейлора с сохранением слагаемых первого порядка и последующего аналитического интегрирования. Постановка задачи Требуется найти решение уравнения Лапласа Au (х ) = 0, х eQ , (1) в области ограниченной осесимметричной поверхностью S с образующей Г (рис. 1) при следующих краевых условиях: u (х ) = g (х), х e S1; (2) q(х) = -ди(х1 = h (х) , х e S2 , (3) dn где S = S1 u S2, q - поток, n - внешняя нормаль к S, g (х) и h (х) - известные распределения потенциала и потока на границах S1 и S2 соответственно. Причем функции g (х) и h (х) не зависят от угловой координаты 9. Рис. 1. Область решения и система координат (r,9,z) Гранично-интегральная формулировка задачи Потенциал, создаваемый в точке х источником единичной интенсивности, сосредоточенным в точке 4, может быть найден из уравнения АО (х, 4 ) = -S (х - 4), где 5 (х - 4) - дельта-функция Дирака. Решением этого уравнения является ньютоновский потенциал О (х, 4 )= -П-, (4) 4пр где р - расстояние между точками х и 4. Тогда, если считать, что потенциальное поле создается фиктивными источниками, непрерывно расположенными на границе S, то суммарный потенциал, создаваемый в точке х всеми фиктивными источниками, можно получить интегрированием выражения (4), предварительно умножив его на плотность распределения фиктивных источников ф (4): u (х) = {ф(4)О(х,4)d% , (5) где обозначение dS^ означает, что при интегрировании именно точка £ «пробегает» поверхность S. В дальнейшем индекс V в подобных случаях будет опущен, однако подразумеваться будет то же самое. С учетом того, что dS^ = rvd9dГ , выражение (5) преобразуется к виду 2п u (х) = 11 ф(£)G(х,£)rvd0dГ , Г 0 где Г - образующая поверхности S, расположенная для простоты в плоскости r0z, при этом 9 x = 0 . Так как задача осесимметрична, то rv и ф (£) не зависят от угла 9 и могут быть вынесены за знак соответствующего интеграла: u (х) = |1 | G (х, £)d9 г4ф (£)dГ . Г V 0 J Введя обозначение 2п G (х, £)= | G (х, £)d9 , 0 выражение (5) окончательно можно записать в форме u (х) = |G(х, £)^ф(£)dГ . Г Как видно из рис. 1: V(zx - zv )2 + К + r2 - 2rxrv cos 9 где rx , zx, r, zv - радиальные и аксиальные координаты точек х и £ соответственно. Тогда G (^ £)= ^ К d+ . . (6 Пусть ' = (zx -zv) + rx2 + b = 2^, (7) 2b m = a + b Тогда формула (6) принимает вид 2п (8) G (х, £ ) = -1 f-^L = 4п 0 va - b cos 9 Wa + b 0 K ( ) ? d9 й й где K (m) = I . - полный эллиптический интеграл первого рода. 0 V1 - m sin2 9 Для рассматриваемого осесимметричного случая поток в точке х еГ в направлении внешней нормали n (х), обусловленный наличием единичного сосредоточенного источника в точке 4, определяется выражением W( n dG (dG +5G F (x, 4) =--= -I-nr +--nz dn ^ drx dzx z где nr и nz - компоненты внешней нормали n (x). Совершенно аналогично предыдущему: q(x) = JF(x,4)г?ф(4)dГ , где F (x, 4 )= | f (x, 4 )de. Нетрудно получить, что 2 2 r2 - rx2 + 1 1 Z - Zj nr + X , E (m)nz a - b F (x,4)= E (m)- K (m) (9) W a+b a - b 2r где E (m )= J Vl - m sin2 ede - полный эллиптический интеграл второго рода. о Таким образом, гранично-интегральная формулировка задачи с учетом условий (2) и (3) представляется в виде системы уравнений (10) (11) g(x) = JG(x,4)гф(4)dГ, x e ^i; г h(x) = JF(x,4)гф(4)dГ, x e S2, где G (x, 4) и F (x, 4) определяются выражениями (8) и (9), ф (4) - искомая плотность распределения фиктивных источников. После нахождения функции ф (4) становится возможным вычисление потенциала и потока в произвольном направлении в любой точке области Дискретизация границы и метод решения Для решения интегральных уравнений (10), (11) образующая Г разбивается на N граничных элементов. Используются «постоянные» граничные элементы, т. е. элементы считаются отрезками прямых, а величина ф внутри каждого элемента считается постоянной (рис. 2). Тогда интегралы (10) и (11) преобразуются в конечные суммы: N g (x )=Е ф, JG (-^4 )rjd г; (12) i=1 г,. N (13) h(x) = £ф, J F(x,4)rjdг , i=1 г,. где г, - i-й граничный элемент, ф, - значение функции ф (4) на ,-м граничном элементе. Теперь, последовательно располагая точки наблюдения х в узлах (т. е. серединах) граничных элементов и записывая для этих точек одно из уравнений (12) или (13) (в зависимости от принадлежности точки к S1 или S2), можно получить систему из N линейных алгебраических уравнений вида N gj = X Gj ф* ; (14) i =1 N h} = X F Ф* (15) i =1 относительно N неизвестных ф*, где gj и hj - известные значения потенциала или потока в узле j-го граничного элемента, Gj = fG (j, \)dГ , F = fF (j, \)dГ , гг Xj - узел j-го граничного элемента (рис. 2), а G (х, \) и F (х, \) определяются выражениями (8) и (9). Значения ф* находятся из системы (14), (15). После этого потенциал и поток в любой точке х eQ могут быть найдены по формулам NN u (х) = X ф* IG (х, ^ kd г, q (х) = X Ф* IF (х, ^ h dr. i =1 Г i =1 Г Вычисление интегралов Gj , Fj Из выражений (8) и (9) видно, что в интегралах Gj и Fj подынтегральные функции содержат эллиптические интегралы первого и второго рода, поэтому их аналитическое интегрирование не представляется возможным. В связи с этим их значения определяются численным интегрированием с помощью квадратурной формулы Гаусса. В отдельном рассмотрении нуждается случай, когда точка наблюдения х расположена на элементе, по которому происходит интегрирование. Эти интегралы являются сингулярными, поскольку фундаментальные решения имеют особенности в точке, где расположен источник. Для вычисления интеграла Gii используется координата n (рис. 3). Начало оси П располагается в узле элемента. Тогда от интегрирования по Г можно перейти к интегрированию по n в пределах от -Ц до Li, где Ц - половина длины i-го элемента. Рис. 3. Локальная координатная ось n Из рисунка видно, что Г - r = n cos a = -ner; (16) z^ - zi = n sin a = nez, (17) где ri, zi - координаты точки xi (узла i-го элемента), n - координата точки 4, er, ez - компоненты орта e. Выражения a, b и m, определяемые формулами (7), при x = xi с учетом (16) и (17) приобретают вид a = 2r2 + 2nrer + n2, (18) 2b b = 2r + 2nrer 4r + 4nrer m = a + b 4r2 + 4nrer + n2 и являются функциями переменной n. При подстановке (18) в (8) с учетом (16) интеграл Gii записывается как f 4rf + 4nrier ^ 4r2 + 4nrer + n2 •(ner + ri )n K 1 L G = - J П -Li er + n2 V4r2 + 4lir, Выражение, стоящее под знаком интеграла, является произведением двух функций: 1 L - Gn = - J K (n) A (n)dq , TT •> ( 4r/2 + 4щег ^ 4rf + 4щег + n2 ner + r/ ^r2 + 4щег + n2 Эллиптический интеграл (20) можно аппроксимировать многочленом [7]: K £ n - . K (n) = K A (n) = (20) (21) (22) где а функция A (n) после разложения в ряд Тейлора с учетом слагаемых порядка не выше n имеет вид 1 е A(n)~- +-n • V ' 2 4r (23) Тогда из (19) находится приближенное значение G// : 1 L ^ „ Л ^ „ Л т ( п e2 L2 L 1 + bh. - in L 1 er Ь L - + -n I dn =- 2 4r. J п er , n -n - ln- 2r/ 8r G/ (24) 8r. 12r Интеграл F. имеет значение 0.5 в соответствии с теоремой о разрыве производной потенциала простого слоя [8]. Решение тестовых примеров Для подтверждения работоспособности метода приводятся результаты численного решения трех тестовых примеров, имеющих аналитические решения. Погрешности оцениваются с использованием метрики пространства L2, т. е. по формулам Е = /-]Г()-и(А))2 и Е = ±]Г(Ч)-qj7)2 (Ч) .(Ч) где u и q. (ось цилиндра), r = 0.5 (потенциал вычисляется в 101 точке: zt = 0.01-[i -1], 1 < i < 101) и r = 1 (боковая поверхность, вычисление потенциала в узлах элементов). На рис. 5 приведены графики зависимости погрешности для потенциала от числа граничных элементов. Причем N означает число элементов на каждом из трех участков границы, т.е. общее число элементов - 3N. Видно, что на границе решение сходится медленнее, чем внутри области решения. Это может быть объяснено тем, что при вычислении потенциала в узлах необходимо вычислять сингулярные интегралы, которые, по-видимому, определяются с меньшей точностью, чем несингулярные. Вообще говоря, эти интегралы имеют слабые особенности порядка логарифмической и поэтому могут быть посчитаны непосредственно по квадратурной формуле Гаусса. Рис. 6 показывает незначительные отличия в результатах при вычислении этими двумя способами, что подтверждает правильность формулы (24). Также проведен расчет потока на прямой z = 1 (верхнее 0 10 20 30 40 N Рис. 5. Зависимости погрешности E для потенциала от числа элементов N для примера 1: точки - r = 0 , сплошные линии - r = 0.5 , пунктирные линии - r = 1. Черный цвет - формулы Гаусса с четырьмя точками, серый цвет - формулы Гаусса с восемью точками E 0.012 0.009 0.006 0.003 основание цилиндра, вычисление потока в узлах элементов). Значение потока для внутренних узлов (т. е. всех узлов, кроме ближайшего к оси и ближайшего к боковой поверхности) сходится к аналитическому решению, однако значение в при-осевом узле отдаляется от аналитического (рис. 7). В крайнем правом узле погрешность еще в 10-15 раз больше, чем в приосевом (погрешность в узле рассчитана по формуле E = (Ч^ - ). Такое поведение решения для потока обусловлено наличием особенности в угловой точке. 0 10 20 30 40 N Рис. 6. Зависимости погрешности E для потенциала от числа элементов N для примера 1: точки - при использовании формулы (24), сплошные линии - при использовании формулы Гаусса. Черный цвет -формулы Гаусса с четырьмя точками, серый цвет -формулы Гаусса с восемью точками 0 10 20 30 40 N Рис. 7. Зависимости погрешности E для потока от числа элементов N для примера 1: сплошные линии - во внутренних узлах, пунктирные линии - в приосевом узле. Черный цвет - формулы Гаусса с четырьмя точками, серый цвет - формулы Гаусса с восемью точками Пример 2. Рассматривается область между двумя концентрическими сферами. На внутренней, единичного радиуса, задан потенциал 0; на внешней, с радиусом 2, потенциал равен 1 (рис. 8). Точное решение для потенциала u = 211- 1 Vr 2 + z2 Рис. 8. Постановка задачи для примера 2 Численный расчет проведен на двух прямых: r = 0 (вычисление потенциала в 101 точке: zi = 1 + 0.01-[i -1], 1

Скачать электронную версию публикации