Брэббия К., Теллес Ж., Вроубел Л. - Методы граничных элементов (1050609), страница 27
Текст из файла (страница 27)
Н.бг! л=! Ряд в выражении (4.63) сходится очень быстро при малых значениях а и медленно — когда а приближается к величине порядка единицы. Чтобы преодолеть это затруднение, интегрирование по отрезку, лежащему вблизи сингулярности, можно выполни~ь аналитически, поскольку при этом коэффициент а всегда меньше или равен единице (что гарантирует сходимость ряда с требуемой точностью при удержании не более шести членов ряда), а численное интегрирование с использованием стандартной квадратурной формулы Гаусса — по остальной части длины элемента, как если бы эти части были отдельными элементами. Длина Ь части элемента, допускающей аналитическое интегрирование, определяется с учетом выражения (4.61) из неравенства Ь ~ 4 (йег(„-)!ег ( !. (4.64) Подобный подход можно распространить как на линейные функции, так и-на функции более высокого порядка.
Следует отметить, что используя некоторую интерполнрующую по пространственным координатам функцию и повышая порядок аппроксимации по времени, можно получить лишь дополнительные регулярные вклады в граничные интегралы. Диагональные коэффициенты Ннгм в уравнении (4.40) содержат интегралы с особенностями более высокого порядка, интегрируемые только в смысле главного значения. Эти коэффициенты можно найти аналогично тому, как это делалось при рассмотрении уравнения Лапласа, воспользовавшись равномерно распределенным по всему телу потенциалом, однако наличие интегралов по области делает этот процесс неэффективным.
Поэтому оба слагаемых с! и Нп (Нм + Нп) необходимо вычислять непосредственно. Значения свободных членов с получаются, как и в случае уравнения Лапласа, с помощью предельного перехода, что дает е, о! = ! — — 1пп ) — ехр ь! — г! в е(8 = 1, г 1 Г вг ! и+ а! — аа е-~о а 4)е (!г — !о) 2п (4.65) В случае постоянных и линейных элементов члены Нп тождественно равны нулю, поскольку г и и ортогональны, что приводит к равенству !! = 0 в выражении (4.37).
В случае элементов 175 174 Глава 4 Задачи апории ззкплопроводпаспзи более высокого порядка это условие не выполняется и интегрирование следует выполнить аналитически (хотя бы по короткому прямолинейному отрезку, включающему сингулярность), с тем чтобы соответствующим образом выделить главные значения интегралов. 4.6. Шаговые по времени схемы Основное различие между представленными здесь двумя щаговыми по времени схемами лежит в способе, которым значения переменных, относящихся к данному моменту времени, учиты- Г ваются при получении решения для нового момента времени. l В 1-й схеме они учитываются в интеграле по области как псевдоначальные значения, тогда ° го,н„„„„, как во 2-й схеме их изменение узлы учитывается путем суммирова° внуузсннне ния граничных интегралов.
1-я схема использовалась в работах Рис. 4.4. Разбиение областийз и еегра. [2 — 6, 20, 21[; 2-я схема была иицы Гна граничные алементы и ячей-, предложена Талером и Муллером [221 в сочетании с непрямой формулировкой еще в 1970 г. и впоследствии использовалась в работах[13, 19, 23[. Начнем с более детального рассмотрения 1-й схемы. В начале процесса (при 2 = (о) устанавливаются начальные значения иа функции и в области й1. Область разбивается на /. (треугольных в данном случае) ячеек (рис. 4.4).
Начальные условия берутся в соответствии с результатами численного интегрирования по области и темя их значениями, которые задаются в некоторых рассматриваемых внутренних точках, Поскольку половина граничных значений функций и и д задана, для определения остальных граничных значений на первом шаге по времени (Р = 1) можно использовать уравнение (4.38), Отметим, что при использовании линейных (или более высокого порядка) интерполяций функций по времени должны быть также заданы значения функции у7 на границе Г (см.
уравкение (4.50)). В конце каждого шага по времени в предварительно выбранных внутренних точках снова определяются значения функции и, с тем чтобы использовать их в качестве начальных условий на следующем шаге. Это можно сделать с помощью уравнения (4.29), которое в матричной форме принимает (например, в случае постоянных элементов) вид й — б'4~о — Н'с/и+ В'й„,.
(4.66) 0,1 Коэффициенты матриц 6 и Н, входящих в уравнение (4.42), и матриц 6' и Н' нз уравнения (4.66) зависят от геометрии, свойств среды и величины шага по времени. Так, если применяется постоянный шаг по времени, их можно вычислить только один раз н записать в память. То же самое применимо для коэФфициентов матриц В и В', обусловленных интегрированием по ячейке. Для численного интегрирования по ячейке можно применить квадратурную схему Хаммера [24 [ и допустить, что значения функции // в каждой точке интегрирования вычис- 00 ляются непосредственно; тогда для коэффициентов матрицы В получаем следующие выражения: [2[з 2 / суй 1 В, = 4,З 'х Р ~ 4йау / ' 0,2 (4.67) 0.5 где [./[ — якобиан, гау2 — ве- 0,0 совой коэффициент в квадра- 1,0 турной схеме Хаммера.
С другой сгороны, можно 0 зз цз 00 ол предполозкить, что функция и изменяется по области, вани- Рис. 4 5. Зависимость ФУнкции и* от г ллн различных значений шага ио времаемой ячейкой, всоответствии с выбранной формой интегрирующей функции, и вычислить коэффициенты Воь используя, например, полуаналитическую схему интегрирования 1251. Хотя в литературе имеются сообщения об очень хороших результатах, полученных при использовании этой схемы, необходимо соблюдать осторожность при выборе шага по времени.
Поскольку при Гзг-ы 0 подынтегральная функция в интеграле по области (фундаментальное решение) становится все менее гладкой, превращаясь в пределе в дельта-функцию Дирака (рис. 4.5), трудность численного интегрирования функции с указанным характером поведения может породить проблемы численного характера (см. работы [21, 26[). Во 2-й схеме не требуется находить значения функции и во внутренних точках в конце каждого шага по времени. Интеграл по области, обусловленный учетом начальных условий при 1 = йш требуется только в том случае, если ио ~ О. Кроме того, если урзио = О, то интеграл по области можно преобразовать в эквивалентный граничный интеграл. Так же как и в большинстве практически важных задач, снижение размерности задачи в данном случае является важным достоинством.
Поскольку с ростом числа шагов по времени увеличивается количество граничных ин- 177 Задачи теории тенлопрооодноегли Глана 4 176 ~ иоив Нй = — „~ ( —, ехр ~ — д (,, ) ] ио— я г 1 (Г. 2 11 4/г(!р — 1о) 1 ~о) (4.71) тегралов, для эффективной работы численной процедуры необходимо использовать специально подобранные схемы численного интегрирования. Для того чтобы пояснить сказанное, вернемся к уравнению (4.40).
Как видно из этого уравнения, при вычислении неизвестных граничных значений в момент времени 1 = 1р требуется определять матрицы бур и Игр для / =- 1, 2, ..., Р. Матрицы от 6,„ до 61р г! р и от Нар до Н(р пи будут последовательно умножать заданные или вычисленные на предыдущем шаге по времени значения функций и и д, представленные в форме вектора с неизвестными компонентами. В силу характера зависимости от времени подынтегральных функций в выражении (4.41) допустимо прн подсчете вклада матриц, соответствующих начальным шагам, брать меньшее число точек в квадратурной формуле Гаусса. Заметим, что если в процессе численного решения задается постоянный шаг по времени, то на каждом шаге требуется определять только две новые матрицы, тогда как все остальные сохраняются в памяти ЭВМ. Важно отметить, что в методе граничных элементов благодаря его простоте не возникает проблем устойчивости решения.
Действительно, математическое доказательство равномерной сходимости и устойчивости численных расчетов по методу граничных элементов в приложении, например, к нелинейным двумерным нестационарным задачам теплопроводности недавно было дано в работе !27 1.
С помощью второго тождества Грина интеграл по области можно преобразовать в эквивалентные граничные интегралы для случая, когда ио — гармоническая функция: ~ иоЧ'[/г(й = ~ [чио — „— (/ ~ ' ) г(Г. (4.68) я Р Поскольку интеграл по области в соответствии с уравнением (4.34) имеет вид ~ и ив!114, (4.69) функцию (/ следует взять такой, чтобы выполнялось равенство (/в(/ = ив. Одну из таких функций можно легко найти: (/ =- ~ — ', Ц.
* Ь) / = —,'„Е, ( — „) . (4.70) Тогда из соотношения (4,68) получаем Таблица 4.2. Температура (нб) в момент времени и = 1,2 ч при шаге по времени 61= 0,1 ч Реи~ение но 2.му методу греннчных элементов Решение но !-м>' методу греиичннх элементов Решение оо методу ионеч- ных элементов Анелнтнчеевое решение .ет — 16,607 — 16,664 — !Б,830 — !7,089 — 17,4!б — !6,718 — 17,011 — 17,373 — 17,666 0,0 0,0 0,0 0,0 0,0 0,3 0,6 0,9 1,2 — 16,6Б2 — !6,716 — 1Б,87! — !7,114 — 17,423 — 16,766 — 17,042 — 17,383 — 17,666 — 16,673 — 16,726 — 16,881 — 17,122 — 17,428 — 16,777 — 17,050 — ! 7,389 — 17,667 0,0 0,3 0,6 0,9 1,2 0,3 0,6 0,9 1,2 — ! 6,771 — 16,820 — 16,963 — 17,186 — 17,467 — ! 6,867 — 17,119 — 17,430 — 17,682 Здесь г/е = див/дп, е( определяется из выражения (4.37).
Приведенные выше интегралы можно найти численно, за исключением сингулярных слагаемых, для которых интегралы вычисляются аналитически, как об этом говорилось в предыдущем разделе. Пример 4.2. Предметом исследования здесь является распределение температуры в квадратной области размером ЗхЗ м с начальной температурой и, =- — 1,111'С и коэффициентом теплопроводности, равным й =- 0,204 Вт/(м 'С), для которой задано граничное условие и — 17,77 'С на границе Г при 1 > 1,. Эти исходные данные были выбраны для того, чтобы иметь возможность сравнить полученные результаты с решением методом конечных элементов [28 1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
