Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 104
Текст из файла (страница 104)
И„„=у... Ч у, у. В,т (8) Здесь А" означает, что матрица А вычислена для сетки л х л. Несложный анализ (хотя бы размерностей) показывает, что можно вычислить «уннверсальную» матрицу на сетке 1 х 1 один раз и после этою пользоваться формулой А" = УР 'А'. В универсальной матрице А' наиболее сложно вычисляются элементы Ау с малыми значениЯми )У! + )У'~ цЗ, так как именно в этих случаях сказывается сингулярность подынтегрзльной функции.
Для вычисления таких элементов были разработаны специальные программы аккуратного интегрирования, с помощью которых вычислены элементы А, для малых у, у (их можно найти в соответствующих работах). Остальные элементы А, . легко вычисляются по су простым асимптотнчсским формулам. При этом используется очевидное свойство симметрии: А,.
=А, =А,, позволяющее хранить в памяти толъко «четверть» матрицы (у, у=0, 1, ..., К, где уь — целое число, такое, что уи превосходит линейный размер СР). Итак, составление системы уравнений (8) проблем не содержит, если известны значения А, у нри малых ! У~ + !у~. Следующий вопрос — решение этой системы.
Число уравнений не так уж мало. В расчетах, которые обсуждаются ниже, число узлов было от 500 до Очевидно, число уравнений равно числу неизвестных. В такой форме можно считать рь м любым базисом, а не только тем, который быя описан выше. Мы получили систему уравнений (7) с четырехиндексной матрицей. Хранение ее в памяти ЭВМ, если число базисных функций не очень мало, — проблема, впрочем, для современных ЭВМ (серии ЕС, например) уже преодолимая. Но вот вычисление матрицы остается проблемой: ведь каждый элемент А — интеграл по 0 (хуже того, интеграл от функции с достаточно сильной особенностью), Именно это заставляет сузить общую галеркннскую конструкцию, используя специальные конечные элементы.
Заметим, что ядро интегрального уравнения (1) зависит не от четырех аргументов (х, у, х', у'), а только от двух (х — х', у — у'), что дает основание ожидать соответствующею характера матрицы: АВ '" = АВ В, у Это весьма полезное свойство действительно удается получить,но только за счет использования одинаковых (с точностью до сдвига аргументов) конечных элементов, Очевидно, рь,„(х, «) = р (х — х, у — у,„). 493 ПСВВЦОДНЬЬЕГЕНЦНАЛЬНЫЕ ГГ4ВНЕННН 1000, причем матрица системы полностью заполнена.
Применение стандартных программ решения линейных уравнений потребовало бы порядка 109 операций, 104 памяти, что уже представляет серьезную проблему для современных машин ЕС (типа 1040, 1045, 1050), не говоря уже о БЭСМ-б, для которой зто просто непосилышя задача (а именно эта ЭВМ применялась в расчетах). Достаточно эффективным средством решения системы (8) является простейший итерационный метод, известный под названием «релаксационный», Он состоит в поочередном пересчете и по формуле ...-ч' ~у.. -~',~,...„). (9) ь/ Здесь в сумме по 1, у' пропускается слагаемое 1= к, у' = ш. Используется также ускорение сходнмости методом сверхрелаксации.
По формуле (9) вычисляется «предварительное» значение й, окончательное значение находится по формуле иь = ив + «о(й„— и» ). Параметр гл определялся экспериментально. В табл. 22 приведена эволюция нормы невязки уравнения (8) в зависимости от номера итерации 9 при разных значениях ш. Видно, что оптимальное значение ю 1.4, при этом скорость сходимости достаточно высока.
Начиная с тривиального приближения и,. ш О, за 15 итераций можно получить невязку, существенно меньшую,г'. Такое приближенное решение системы (8) можно трактовать как точное, оютветствующее незначительно измененной правой части. Таблица 22 Возникает вопрос о точности расчета прн приемлемом числе узлов сетки. В самом деле, 1000 узлов на двумерную область — это не так уж много (примерно по 30 узлов на линейный размер области).
Уменьшить шаг л для повышения точности здесь не так просто: вы- 494 пгязлюханныв матовы вычислительной ьизики числение суммы в правой части (9) для каждого узла (Ь, т) «стоит» О(Ь з) операций. Число узлов О(Ь ~), да и число итераций можно оценивать как 0(Ь ') или, в лучшем случае, О(Ь пз). Вышеприведенные оценки сделаны по аналогии с хорошо изученным уравнением Лапласа. Там мы имеем дело с оператором второго порядка, для которого обусловленность матрицы конечномерной аппроксимации есть 0(Ь з), сходимосгь простых итераций имеет скорость 1 — О(Ьз) (сверхрелаксацня при оптимальном ю доводит коэффициент сходимости до ! — О(Ь)).
Мы же имеем дело с оператором первого порядка, с числом обусловленности матрицы' А порядка 0(Ь '). Эксперимент подтверждает это предположение. Итак, «стоимость» приближенного решения есть в лучшеы случае 0(Ь «л) операций и она резко возрастает прн незначительном уменьшении шага (например, раз в 25 при переходе от Ь к Ь/2). На БЭСМ-6 расчет трещины с числом узлов порядка 800 стоит около 5 мин машинного времени. Достаточно ли доступного нам числа узлов для получения решения с нужной для приложений точностью? Следует иметь в виду, что требования к точности не очень высоки: погрешность 1 +5;/ допустима для многих технических приложений.
Конечно, ответ зависит от дифференциальных свойств решения, последние, в свою очередь, — от гладкости правой части и свойств оператора в (1). В прикладных расчетах ~, как правило, — не очень сложная функция. Это понятно. Ведь трещина есть обьект малых размеров, находящийся «в глубине» некоторой конструкции. Функция У создается системой внешних нагрузок. «Создать» функцию У, сильно меняющуюся на малом расстоянии, не так-то.
просто, тем более что никто к этому специально не стремится, Псевдодифференциальный оператор в (1) по своей природе близок к «эллиптическому дифференциальному оператору первого порядка», т.е. при его обращении гладкость решения повышается на один порядок по сравнению с /. Поэтому есть основание ожидать, что искомое решение и(х, у)— достаточно просго устроенная функция. Действительно, расчеты простых задач (Π— эллипс, У = сопш) с известным точным решением показывают, что и(х, у) есть колоколообразная функция простого вида, и при числе узлов 20 - 30 на линейный размер области точность приближенного решения очень высока почти всюду в области О. Исключение составляет узкая (шириной (2 + 3)Ь) полоса узлов, примыкающих к контуру дб.
В этой полосе погрешность составляет 1б-ь 30 , вне ее — погрешность около 1 /; и чем дальше от границы, тем она меньше. Это обстоятельство не случайное, и ею можно было бы предвидеть, используя достаточно развитую математическую теорию уравнения (1). Основной результат, который необходимо учесть при построении численном метода, состоит в следующем. При гладкой правой 495 % зо1 псккдоди«»ггкнци»лькык гглккккия части у решение и есть непрерывная функция, обращающаяся в нуль на дб, гладкая всюду, кроме окресгности контура дб.
Известен и характер и(х, у) в этой окрестности. Вводя в точках контура локальную систему координат (»„Ч), где г — длина дуги на контуре, Ч вЂ” расстояние от контура по внутренней нормали, имеем для и асимптотику и(К Ч) = Л1® Ч + О(Ч ).
(10) Иными словами, и имеет на контуре корневую особенность. Теперь понятны источники погрешностей численного решения вблизи контура: негладкость самого решения и грубая аппроксимация контура <ступенчатой» линией. Именно такой является граница счетной области бк, на границе которой решение в форме (6) обращается в нуль. Казалось бы, можно пренебречь полученной погрешносп ю: ведь в подавляющем числе узлов точность приближенного решения высока. Однако она высока для величин, особенного интереса для приложений не представляющих. Дело в том, что наиболее важным для приложений является так называемый коэффициент концентрации напряжений на контуре трещины дб — это функция МД), входящая в асимптотику (10).
Именно ради определения Ф(г) затеваются расчеты, так как этой функцией определяется дальнейшая судьба трещины. Будет ли она расти и с какой скоростью (т.е. представляет ли трещина опасность для конструкции), через какое время ее рост приведет к разрушению тела с трещиной, и т.п.? Итак, получив численное решение, нужно проанализировать ход изменения ик вблизи контура и извлечь из него оценку М(г). Ясно, что при этом анализе нужно отступить от контура внутрь области на (2+ 3)л, однако там (если шаг Ь не очень мал) уже начинают сказываться величины порядка О(Чкп).
В трещинах простой формы такой анализ давал неплохие результаты, Но при переходе к трещинам «произвольной» формы ситуация осложняется и требуется уточнение описанной выше расчетной схемы. Кстати, поясним смысл краевого условия «и = 0 на дб». Если подействовать оператором (1) на функцию и с асимптотикой (1О), получим функцию /(х, у), непрерывную в С вплоть до границы. Если же и(х, у) непрерывна в б н не равна нулю на дб, /(х, у) обращается в бесконечность на дб. Метод граничных сеток.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
