Основы метода конечных элементов

Интегрируя уравнение (VI 1.39) на отрезке [хп xc+\]f получаем

Заменяя интеграл в правой части (VII.41) различными формулами численного интегрирования, получают различные методы численного решения задач Коши, обладающие теми или иными качествами. За­ мена интеграла по формуле прямоугольников приводит нас к явной

У( + 1 = yi + hf(xc+1, yi+ О

по формуле (VI 1.42) можно вычислить у19 у2 и т. д. до удг.

Для численного решения задачи неявным методом Эйлера, зная начальное условие (VI 1.43), на каждом шаге интегрирования можно использовать итерационную формулу

дится к соответствующему решению нелинейного уравнения. Началь­ ное приближение выбирают, используя явную формулу метода Эйле­ ра. Численные методы решения обыкновенных дифференциальных уравнений, которые для вычисления значения ус+\ используют зна­ чение уiy называют одношаговыми. Явный и неявный метод Эйлера принадлежит к классу одношаговых численных методов.

Ввычислительной машине реализуется не схема (VII.42), (VII.43),

анекоторая возмущенная схема

Для вычисленного значения уы справедлива оценка (см., например,

где О (Л) — погрешность дискретизации.

Из оценки (VI 1.44) следует, что при h О ошибка дискретизации стремится к нулю, а погрешность реализации схемы на ЭВМ — к бес­ конечности. Таким образом, при решении задач Коши на ЭВМ полу­ ченное машинное решение всегда будет отличаться от математического решения. Из (VII.44) следует также, что выбор шага интегрирования является принципиальным моментом. Шаг интегрирования определя­ ется соображениями не только близости решения разностного уравне­ ния к решению дифференциального, но и связанными с влиянием воз­ мущений, возникающих в процессе вычислений.

При численном решении задач с начальными условиями применя­ ют вычислительные схемы, которые имеют вид конечно-разностных уравнений, определенных на сетке соЛ. Чаще всего это линейное раз­ ностное уравнение с постоянными коэффициентами вида

тт

Этому уравнению ставят в соответствие характеристическое уравнение

Решение разностного уравнения (VII.45) называют численно устойчивым, если все корни характеристического уравнения (VI 1.46) по модулю меньше или равны единице и корни, по модулю равные еди­ нице, простые. Численная устойчивость разностной схемы решения задачи Коши не означает уменьшения начального возмущения от шага к шагу интегрирования, а только ограничивает рост этого возмуще­ ния, если он происходит. Решение разностного уравнения (VII.45) называют асимптотически численно устойчивым, если все корни его характеристического уравнения (VI 1.46) по модулю меньше единицы. И наконец, решение численно неустойчиво, если по крайней мере мо­ дуль одного из корней больше единицы.

Рассмотрим явную схему метода Эйлера (VII.42), (VII.43), приме­ ненную к тестовому уравнению (VI 1.34)

y {+ i = у \ + M iyt, i = 0, 1, 2,

Требование асимптотической численной устойчивости в этом случае приводит к условию

где а = Re X, р = Im X. Областью асимптотической численной устой­

чивости в этом случае будет внутренность круга единичного радиуса с центром ha = — 1, Яр = 0 на плоскости {ha, Яр). Для чисто мнимого

значения X (а = 0) неравенство (VI1.47) не выполняется ни при каком значении Я > 0. Отметим, что явный метод Эйлера численно устойчив. Действительно, если для него составить характеристическое уравне­ ние

то оно имеет корень г = 1, что показывает численную устойчивость явной схемы.

Рассмотрим теперь численную устойчивость неявного метода Эй­ лера, который, будучи применен к тестовому уравнению (VII.34), записывается в виде

Уi+i = yi + hly{+и

I d - Я Л ) - 11 < 1,

или

(1 — Л<х)2 + /гф2> 1.

Областью асимптотической численной устойчивости этого метода бу­ дет вся плоскость (ha, /ф), за исключением единичного круга с центром

(1, 0).

Численные методы решения обыкновенных дифференциальных уравнений, у которых область асимптотической устойчивости при ре­ шении тестового уравнения (VI 1.34) содержит всю левую полуплос­ кость плоскости (ha, /ф), называют А-устойчивыми [132]. Иными сло­ вами, А-устойчивыми^ называют численные методы, если их область асимптотической устойчивости включает всю полуплоскость Re (Нк) < ; < 0. Таким образом, неявный метод Эйлера является Л-устойчивым.

Отметим, что в жестких дифференциальных уравнениях влияние на решение погрешности в задании исходных данных убывает со вре­ менем. Применение для решения таких задач численно устойчивых методов, в которых может ограниченно расти погрешность машинной реализации, нецелесообразно. Уменьшение погрешности от шага к ша­ гу обеспечивается использованием /4-устойчивых численных методов.

Рассматриваемые вопросы возникают и для других схем, обладаю­ щих более высоким порядком точности численного интегрирования обыкновенных дифференциальных уравнений, а также для систем таких уравнений (см., например, [80, 126]).

Проиллюстрируем важность рассматриваемых нами теоретических вопросов на примере решения следующей модельной задачи Коши для системы обыкновенных дифференциальных уравнений:

= — 0,01 Зы2 — IOOOUJUJ — 2500ы1ы3, « х (0) = 0, = — 0,013ы2 — 1000ы1Ы2, и2(0) = 1,

- 2 5 0 0 ^ 3 . «,(<>)= 1, 0 < * < 5 0 .

Эта задача решалась на ЭВМ МИР-2 с длиной машинного слова 12 десятичных знаков по программе метода Кутта — Мерсона с автомати­

= 50 было получено за 102 ч работы машины, его относительная по­ грешность 2000 %. Шаг интегрирования всегда выбирался автомати­

чески не больше чем 10-3 При исследовании выяснилось, что эта систе­ ма принадлежит к классу жестких. Применение метода решения жест­ ких систем на этой же машине с прежней длиной машинного слова позволило получить всего за 20 мин значение у = [— 1,842 945 749 X

X 10_6, 0,547 654 343, 1,402 343 765]г в точке х = 50. Относительная погрешность этого машинного решения порядка тысячной доли про­ цента.

Решать все системы обыкновенных дифференциальных уравнений по программам решения жестких систем слишком дорого. Ведь на каж­ дом шаге, как мы видим, приходится решать нелинейные уравнения. Поэтому, с одной стороны, целесообразно различать, с жесткой или нежесткой системой нам приходится иметь дело, а с другой — необ­ ходимо учитывать, что в настоящее время разработаны и разрабаты­ ваются методы решения жестких систем, снимающие требования А- устойчивости (см., например, [137]) и, следовательно, уменьшающие затраты на получени решения.

3. Характеристика некоторых методов и программ решения. Разра­ ботка методов решения задач Коши для систем обыкновенных диффе­ ренциальных уравнений всегда была почти такой же популярной, как и разработка методов решения систем линейных алгебраических уравнений. Именно поэтому к настоящему времени имеется большое количество численных методов интегрирования задач Коши 14, 6, 99]. Все численные методы можно подразделить на одношаговые и много­ шаговые, на явные и неявные. Собственно выбор метода интегрирова­ ния обусловливается в основном характером решаемой задачи, пове­ дением ее решения (основную роль в этом играют свойства правых частей системы обыкновенных дифференциальных уравнений), а также желаемой точностью получения решения и опытом специалиста, решаю­ щего задачу До начала 70-х годов предпочтение при выборе численных методов решения рассматриваемых задач отдавалось явным одношаго­ вым методам решения Однако проведенные впоследствии различными авторами исследования показали, что в случае когда система диффе­ ренциальных уравнений относится к классу жестких, ее целесообраз­ но решать неявными методами. В связи с этим началось интенсивное исследование «старых» неявных методов и создание новых методов, позволяющих уменьшить время получения решения на одном шаге интегрирования [1, 90, 126, 137, 1541

При решении задач Коши для систем обыкновенных дифференци­

К настоящему времени создан ряд пакетов программ интегрирова­ ния задач Коши для систем обыкновенных дифференциальных уравне­ ний, как обычных, так и жестких. В основе почти всех существующих пакетов лежит разработанная Гиром программа DIFSUB [137]. Обыч­ ные системы уравнений здесь интегрируются методом Адамса перемен­ ного порядка, жесткие —~ методом Гира, порядок которого может меняться от первого до шестого. Возможность объединения указанных выше двух методов в одной программе обусловлена тем, что в обоих методах используется один и тот же предиктор.

В качестве контроля для выбора шага и порядка метода использу­ ется требование, чтобы ошибка метода на каждом шаге интегрирова­ ния не превышала е. При этом если решение является возрастающей функцией, то контролируется относительная ошибка метода, а при вы­ ходе решения на стационарное значение — абсолютная ошибка. Ошиб­ ка оценивается через поправку корректора.

В исходной информации для программы DIFSUB, кроме е, под­ программы вычисления правых частей системы, границ изменения шага интегрирования и некоторых других параметров, задается указа­ тель метода MF с возможными значениями 0, 1,2. Если MF = 0, то при интегрировании используется метод Адамса; если MF = 1, то используется метод Гира с аналитически заданной матрицей Якоби (при этом необходимо написать подпрограмму ее вычисления); если MF = 2, то используется метод Гира, но при этом матрица Якоби вычисляется автоматически с помощью численного дифференцирова­ ния правых частей.

Как уже упоминалось, на основе идей, заложенных в программу DIFSUB, и самой программы созданы многочисленные пакеты, на­ пример GEAR, GEARBI, GEARS, GEAR1B и др. Отметим, что пакет GEARBI ориентирован на решение жестких систем дифференциальных уравнений с матрицей Якоби, имеющей ленточную структуру или с матрицей Якоби, которая хорошо приближается матрицей ленточной структуры. Пакет GEARS интегрирует системы с разреженной мат­ рицей Якоби.

Как дальнейшее развитие пакетной проблематики создан ODEPACK — систематизированный набор программ для решения обык­ новенных дифференциальных уравнений [144]. В этот набор включено пять программ решения задач Коши для систем обыкновенных диффе­ ренциальных уравнений. Набор разработан в Национальной лабора­ тории им. Лоуренса в Ливерморе (штат Калифорния, США). Приведем состав набора с указанием названия программ и их назначения:

LSODE — предназначена для интегрирования жестких и нежест­ ких систем обыкновенных дифференциальных уравнений, объединяет возможности пакетов GEAR и GEARB\

LSODI— предназначена для интегрирования задач Коши для си­ стем уравнений вида A (t, у) = g (/, у), создана на основе пакета

GEARIB-

LSODES — предназначена для интегрирования разрешенных от­ носительно производных систем дифференциальных уравнений G

разреженной матрицей Якоби, в основе лежит пакет GEARS; системы линейных алгебраических уравнений решаются программой YSMP\ LSODA — предназначена для автоматического выбора метода

интегрирования;

LSODAR — предназначена для решения нелинейных функцио­ нальных уравнений.

Опишем еще один тип пакетов (см. работу [127]). Предназначенные для решения задач Коши для систем обыкновенных дифференциальных уравнений пакеты EPISODE, EPISODEB, DISPL появились из ре­ шения уравнений в частных производных после дискретизации урав­ нений по пространственным переменным. В пакетах EPISODE и EPISODEB используются формулы обратных разностей переменно­ го порядка с переменным шагом по времени для решения жестких систем обыкновенных дифференциальных уравнений и формулы Адам са переменного порядка с переменным шагом для уравнений, не разре шенных относительно производных. Пакет EPISODE предназначен для решения систем дифференциальных уравнений с полной матри­ цей Якоби, a EPISODEB — с матрицей Якоби ленточной структуры*:

В пакете DISPL используется метод Бубнова— Галеркина и Б- сплайны для дискретизации уравнений по пространственным пере­ менным.

Применение метода конечных элементов для дискретизации и после­ дующее решение дискретных задач дает возможность получить при­ ближенные решения ряда задач теории термоупруго-пластичности, а также задач на колебания и устойчивость, связанных с расчетом на прочность узлов, элементов и конструкций в целом.

Для решения реальных научно-технических задач не существует готовых рецептов. Однако, как показывает опыт, методика, изложен­ ная в монографии, облегчает решение практических задач. Действи­ тельно, можно многие научно-технические задачи вполне определенно отнести к тому или иному классу, описанному в данной работе (или исследованному в других работах), а следовательно, и найти путь ре­ шения конкретной прикладной задачи.

Несмотря на большие возможности современных ЭВМ, они все же имеют ограничения как по быстродействию, так и по объему запоминаю­ щих устройств. Об этих ограничениях не следует забывать при реше­ нии прикладных задач. Построение приемлемой математической мо­ дели, которая, с одной стороны, достаточно полно отражает физику протекающих процессов, а с другой — может быть решена на имею­ щейся ЭВМ, осуществляется, как правило, в несколько этапов, путем постепенного уточнения.

Для построения дискретных задач можно использовать как схемы МКЭ, когда в качестве базисных или допустимых функций использу­ ются линейные (полилинейные) полиномы, так и схемы высокого по­ рядка точности, где в качестве базисных или допустимых функций при­ меняются полиномы высших степеней.

Разбиение исходной области на элементы и выбор базисных функ­ ций производятся с учетом требований к точности искомого приближен­ ного решения и теоретических оценок погрешности метода конечных элементов. Следует заметить, что эти оценки устанавливают только ско­ рость сходимости в некоторой норме приближенного решения к точ­

ному (например, О (hy)), без конструктивного и достаточно точного опре­ деления константы при множителе /iv, а потому для первого конкрет­ ного выбора конечного элемента приближенное решение МКЭ может оказаться недостаточно близким к точному решению исходной задачи.

В ряде случаев расчет приходится выполнять на нескольких сетках (для нескольких разбиений), чтобы обеспечить нужную точность решения.

В связи с рассмотрением погрешностей МКЭ необходимо упомянуть также ошибки, возникающие при аппроксимации криволинейной гра­ ницы области определения задачи некоторой кусочно-полиномиальной границей. Поскольку в данной монографии подробно представлены исследования только для одномерных задач, интересующихся этим типом ошибок отсылаем к работам [50, 101] и другим литературным источникам.

В результате выполненной МКЭ дискретизации вопрос о построении приближенного решения в зависимости от исходной задачи сводится к решению либо системы алгебраических уравнений (линейных или не­ линейных), либо алгебраической задачи на собственные значения, либо задачи Коши для обыкновенных дифференциальных уравнений и т. п. Здесь также важно помнить о следующих источниках погрешностей вычисленного на ЭВМ приближенного решения.

Во-первых, это ошибки, возникающие при вычислении коэффици­ ентов и параметров указанных выше дискретных задач, например за счет численного интегрирования, замены начальной функции ее интерполянтом и т. п. Эти ошибки, которые можно рассматривать как на­

Во-вторых, это погрешности, связанные с накоплением ошибок округления в процессе машинной реализации численного метода реше­ ния дискретной задачи. Чтобы уменьшить их влияние (снизить до уровня, не превышающего погрешность МКЭ), необходимо выбирать соответствующую длину машинного слова. В частности, при решении систем линейных алгебраических уравнений МКЭ, для которых апри­ орные оценки числа обусловленности матрицы системы известны, та­ кой выбор, как указано в гл. II, не представляет особого труда.

Выбор численного метода и алгоритма для решения дискретной задачи — далеко не простая проблема, имеющая свои сложности и осо­ бенности. Классификация, учитывающая специфические свойства дис­ кретной задачи, помогает обеспечить эффективное вычисление искомого решения. В качестве иллюстрации коснемся лишь задачи решения си­ стем линейных алгебраических уравнений МКЭ и существующих для этого наиболее эффективных численных алгоритмов. Как известно,

блоков. Вид их определяется исходной задачей, выбором конечного элемента и нумерацией узлов сетки. С переходом к решению дифферен­ циальных и вариационных задач в трехмерных пространствах порядок и разреженность матриц существенно возросли. Для решения таких систем уравнений используются прямые и итерационные методы, раз­ личные алгоритмы которых ориентированы на специальный вид (и свойства) матриц системы. Для повышения эффективности использо­ вания алгоритмов прямых методов, в частности для экономии машин­ ного времени и памяти, широко применяются алгоритмы предваритель­

ного упорядочения матрицы системы. (Краткая характеристика этих алгоритмов и соответствующая библиография приведены в гл. VII.)

Повышение эффективности использования итерационных методов для решения больших разреженных систем успешно достигается в ря­ де случаев применением специальной процедуры предварительной под­ готовки матрицы системы, а именно преобусловливанием. Эта процеду­ ра использовалась в параграфе VI.4 при описании решения приклад­ ной задачи, а краткая библиография по.данному вопросу дана в гл. VII.

Для минимизации объема вычислительной работы при решении различных дискретных задач, возникающих в рамках МКЭ, в послед­ нее время широко применяются многосеточные итерационные процес­ сы. Для ознакомления с этими методами можно воспользоваться их краткой характеристикой и библиографией, указанной в гл. VII.

Как уже упоминалось, априорные теоретические оценки погрешно­ сти МКЭ практически малопригодны для контроля достоверности вы­ численных на ЭВМ решений. Поэтому особый интерес представляют удобные для машинной реализации достаточно надежные апостериор­ ные оценки, не требующие к тому же слишком больших вычислитель­ ных затрат В последние годы разработке именно таких оценок уделя­ ется все больше внимания. Эти оценки особенно важны и потому, что являются основным средством конструирования адаптивных проце­ дур для эффективного вычисления МКЭ решений требуемой точности или для получения решений наибольшей возможной точности в рам­ ках допустимых вычислительных затрат. Сейчас уже общепризнано, что эффективное решение современных прикладных задач практически неосуществимо без использования того или иного адаптивного про­ цесса.

Главные «рычаги» управления адаптивным процессом, предназна­ ченным для уменьшения погрешности метода конечных элементов (т. е. повышения точности аппроксимации искомого решения), следую­ щие: равномерное сгущение первоначально введенной сетки; изменение вида сетки (например, переход к неравномерной); повышение степени кусочно-полиномиальных базисных функций МКЭ или добавление к старому базису функций, описывающих поведение решения вблизи особых точек (т е. переход к другому виду конечного элемента).

Каждый из этих способов управления имеет свои достоинства и не­ достатки, поэтому представляется целесообразным их комбинирован­ ное использование. Однако теоретически обоснованное построение по­ добной эффективной комбинации (даже для линейных задач) остается еще нерешенной проблемой На практике обычно при организации адаптивного процесса ограничиваются пока использованием одного из упомянутых способов, чаще всего это измельчение (адаптация) сетки в зависимости от получаемых в ходе вычислений апостериорных оценок погрешности найденного МКЭ решения.

Весьма важным в рамках данной тематики является и вопрос о соз­ дании программного обеспечения (некоторой программной системы) для машинной реализации адаптивного применения МКЭ при решении

адаптивных процедур решения некоторых классов задач можно, на­ пример, по работе [153] и имеющейся в ней библиографии.

Данная монография посвящена алгоритмическим аспектам МКЭ и может служить основой для построения дискретных задач и эффек­ тивных численных методов их решения на ЭВМ.

Для массовых практических расчетов в прикладных организациях на основе теоретических исследований должны создаваться пакеты прикладных программ. Под пакетом программ понимается комплекс программных средств, включающий как набор функциональных мо­ дулей, характерных для данной прикладной области, так и наборы управляющих и вспомогательных модулей. Этот комплекс должен обес­ печивать пользователю возможность в интерактивном режиме постро­ ить эффективную программу и вычислить по ней решение конкретной прикладной задачи в соответствии со спецификой задачи, а также с учетом особенностей ЭВМ, на которой выполняется вычисление иско­ мого решения.

Пакет программ должен ускорить программирование и получение решения достаточно широкого класса математических и прикладных задач, дать пользователю не только решение задачи, но и некоторую достоверную информацию о точности вычисленного решения.

Пакеты программ, основанные на МКЭ, и создаваемые с их помощью программы должны иметь единый интерфейс с программными средства­ ми систем автоматизированного проектирования и автоматизации науч­ ных исследований, а также единую с этими системами информацион­ ную базу.

В работе [73] приведена некоторая характеристика зарубежных пакетов программ, которые используются в машиностроении, судо­ строении, ядерных и аэрокосмических исследованиях, на транспорте, в строительстве и других прикладных областях. Эти пакеты организо­ ваны на различных принципах, но в алгоритмической основе их ле­ жит метод конечных элементов. Использование таких пакетов приклад­