Решение Нелинейных Статических Задач Конечных Элементов

19/11/2013

Здесь мы начинаем обзор алгоритмов используемых для решения нелинейных статических конечно-элементных задач. Этот материал демонстрируется на примере очень простой 1D задачи и основывается на предыдущей статье Решение Линейных Статических Моделей Конечных Элементов.

Система пружины, прикрепленной к жесткой стене

Рассмотрим систему, показанную ниже: пружина одним концом закреплена на жесткой стенке, а к другому ее концу приложена некоторая сила. Жесткость пружины есть функция расстояния, на которое она растягивается k(u)=exp(u). То есть жесткость пружины увеличивается экспоненциально при растяжении.

Пружина закрепленная на жесткой стенке и приложенной силой

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

f(u)=p-k(u)u=p-exp(u)u

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

Нарисуем график этой функции, имея в виду, что мы пытаемся найти u такое, что f(u)=0.

График, представляющий функцию для нелинейной конечно-элементной задачи

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

 Одиночная итерация задачи по Ньютону-Рафсону

Как можно заметить, мы опять стартуем с начального приближения u_0=0, и оцениваем функцию, f(u_0), и ее производную, f'(u_0). Это дает нам точку u_1. Проверяя, мы убеждаемся, что эта точка не является решением, так как f(u_1) \ne 0. Но если мы продолжим делать итерации по Ньютону-Рафсону, как показано ниже, то становится очевидным, что мы приближаемся к решению задачи. (Для более подробного изучения этого алгоритма, вы можете использовать эту ссылку Метод Ньютона.)

Последовательные итерации Ньютона-Рафсона

Таким образом, нахождение решения нелинейной задачи по существу идентично решению линейной задачи, за исключением того, что теперь мы делаем многократные последовательные шаги по Ньютону-Рафсону, чтобы получить решение. На самом деле, мы можем продолжать делать итерации и сколь угодно близко подойти к решению, но в этом нет необходимости. Как обсуждалось ранее, рано или поздно мы достигнем машинной точности представления чисел на компьютере, так что всегда есть практический предел насколько можно приблизиться к решению. Давайте посмотрим на результаты нескольких итераций:

i u_i |f(u_i)| |u_{i-1}-u_i| |f(u_{i-1})-f(u_i)|
0 0.000 2.000
1 2.000 12.77 2.000 10.77
2 1.424 3.915 0.576 8.855
3 1.035 0.914 0.389 3.001
4 0.876 0.104 0.159 0.810
5 0.853 0.002 0.023 0.102
6 0.852 0.001 0.001 0.001

Мы видим, что после шести итераций разница между двумя последовательными значениями f(u), и u, а также абсолютное значение f(u), уменьшилась до 0.001 или менее. После шести итераций по Ньютону-Рафсону начиная с u_0=0, решение сошлось с погрешностью менее 0.001. При решении нелинейных задач, мы применяем этот алгоритм до тех пор, пока решение не сойдется с желаемой точностью. Имеется и второй критерий для завершения: решатель должен сделать не более определенного количества итераций. Вне зависимости от того, какой критерий выполнится первым — относительная погрешность или число итераций — решатель остановит процесс. Также, необходимо помнить о проблеме численного масштабирования, которую мы обсуждали в статье блога о решении линейных статических конечно-элементных задач. Критерий относительной погрешности применяется к масштабированному, а не к абсолютному вектору решения.

Хотя показать это значительно сложнее, но алгоритм используемый для решения нелинейных задач, в случае когда u является вектором остается тем же, что и для обычной задачи. Однако, при решении задачи с сотнями, тысячами или, даже, с миллионами степеней свободы, желательно, чтобы число итераций Ньютона-Рафсона было минимально возможным. Напомним, что нам необходимо решить задачу \mathbf{u}_{i+1}=\mathbf{u}_{i}-[\mathbf{f}'(\mathbf{u}_{i})]^{-1}\mathbf{f}(\mathbf{u}_{i}), в которой вычисление обратных производных является наиболее ресурсозатратной процедурой. Чтобы избежать вычислений в области значений, заведомо не содержащих решения, и минимизировать число итераций Ньютона-Рафсона, в среде COMSOL используется релаксационный множитель (коэффициент релаксации). Посмотрим снова на начальную итерацию Ньютона-Рафсона и ее график, и заметим, что для нее |\mathbf{f}(\mathbf{u}_{i+1})|>|\mathbf{f}(\mathbf{u}_{i})|. Так что для этой итерации, мы делаем слишком большой шаг. Когда это происходит, COMSOL осуществляет поиск вдоль интервала [\mathbf{u}_{i},\mathbf{u}_{i+1}] точки \mathbf{u}_{damped}=\mathbf{u}_i+\alpha(\mathbf{u}_{i+1}-\mathbf{u}_i) такой, что |\mathbf{f(u}_{damped})|<|\mathbf{f(u}_{i})|. После чего, схема итераций Ньютона-Рафсона возобновляется с этой точки.

Релаксационный множитель (коэффициент релаксации) в среде COMSOL Multiphysics

Параметр \alpha известен, как релаксационный множитель (коэффициент релаксации) и изменяется в интервале 0< \alpha \le 1. При \alpha \rightarrow 0 мы говорим, что релаксация возрастает, а \alpha = 1 значит, что задача не срелаксирована. Этот метод достаточно привлекателен, поскольку требует от среды COMSOL только оценки значения \mathbf{f(u}_{damped}), и вычислительные затраты при этом гораздо ниже, чем вычисление производных \mathbf{f'(u}_{i}) и обратных к ним величин [\mathbf{f}'(\mathbf{u}_i)]^\mathbf{-1}.

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

Почему Решение Нелинейной Задачи Может Не Сойтись

Нелинейные задачи по своей сути сложно решить, поскольку существует множество причин, из-за которых рассмотренная выше процедура решения может не сойтись. Хотя имеется много причин, по которым the Newton-Raphson method can fail (метод Ньютона-Рафсона терпит неудачу) , мы ограничимся обсуждением следующих практических случаев.

Случай 1: Начальное Приближение Слишком Далеко от Решения

Прежде всего, рассмотрим ту же самую нелинейную задачу, но с другой стартовой позицией, например, u_0=-2. Как можно заметить из графика ниже, если мы выберем начальную точку в области u_0\le-1, метод Ньютона-Рафсона не сможет найти решение, так как производные f(u) не дают точки пересечения, приближающейся к решению. Решение не может быть найдено слева от u_0=-1, так что эти начальные точки находятся вне радиуса сходимости метода Ньютона-Рафсона. Т.е. из-за выбора начального приближения метод Ньютона-Рафсона может не сойтись, даже если решение существует. Таким образом, в отличие от линейного случая, когда корректно поставленная задача всегда будет решена, сходимость нелинейных моделей в сильной степени зависит от выбора стартовой позиции. Позже мы вернемся к вопросу, как наилучшим способом выбрать начальное приближение.

Начальное приближение нелинейной модели располагается слишком далеко от решения, что приводит к сбою в сходимости

Случай 2: Задача Не Имеет Решения

Нелинейный решатель также потерпит неудачу, если задача сама по себе не имеет решения. Вернемся опять к задаче, рассмотренной выше, но с жесткостью пружины k(u)=\exp(-u). Другим словами, с растяжением жесткость пружины уменьшается. Если мы построим график f(u) для нагрузки p=2, мы увидим, что решение не существует. К несчастью, алгоритм Ньютона-Рафсона не может определить это самостоятельно; алгоритм просто будет пытаться найти решение и завершит свою работу после заданного пользователем числа итераций.

Жесткость уменьшается при растяжении пружины

Случай 3: Задача Не Гладкая и Недифференцируемая

Наконец, рассмотрим случай, когда свойства материала изменяются скачкообразным образом. Например, рассмотрим такую же систему, как и перед этим, но с жесткостью пружины, принимающей различные значения в различных диапазонах растяжений, значение k=0.5 для u\le1.8, значение k=1 для 1.8<u<2.2, и k=1.5 для u\ge2.2. Если мы построим график f(u) для этого случая, мы увидим, что он недифференцируем и разрывен, а это нарушает требования метода Ньютона-Рафсона. Легко убедиться непосредственной проверкой, что если стартовая точка выбрана не в интервале 1.8<u<2.2, то итерации Ньютона-Рафсона будут осциллировать за пределами этого диапазона.

Материал со скачкообразными изменениями свойств

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

  • Начальное приближение выбрано слишком далеко от решения
  • Поставлена задача, которая не имеет решения
  • Задача негладкая и недифференцируемая

Интерпретация Log-Файла среды COMSOL

Log-файл (системный журнал)

В ближайшее время мы обсудим пути решения всех этих вопросов, но прежде, давайте посмотрим на log-файл типичной нелинейной задачи конечных элементов. Посмотрите на приведенный ниже текст log-файла (дополненного номерами строк) геометрической нелинейной задачи из строительной механики:

1) Stationary Solver 1 in Solver 1 started at 10-Jul-2013 15:23:07.
2) Nonlinear solver
3) Number of degrees of freedom solved for: 2002.
4) Symmetric matrices found.
5) Scales for dependent variables:
6) Displacement field (Material) (mod1.u): 1
7) Iter ErrEst Damping Stepsize #Res #Jac #Sol
8) 1 6.1 0.1112155 7 3 1 3
9) 2 0.12 0.6051934 1.2 4 2 5
10) 3 0.045 1.0000000 0.18 5 3 7
11) 4 0.012 1.0000000 0.075 6 4 9
12) 5 0.0012 1.0000000 0.018 7 5 11
13) 6 1.6e-005 1.0000000 0.0015 8 6 13
14) Stationary Solver 1 in Solver 1: Solution time: 1 s
15) Physical memory: 849 MB
16) Virtual memory: 946 MB

Разъяснения

  • Строка 1 сообщает о типе вызываемого решателя (вычислителя) и времени начала выполнения задачи.
  • Строка 2 сообщает, что вызванное программное обеспечение — это решатель нелинейных систем.
  • Строка 3 сообщает о размере задачи в терминах числа степеней свободы.
  • Строка 4 сообщает о типе матрицы конечных элементов.
  • Строки 5-6 сообщают о параметрах масштабирования. В данном случае масштаб поля смещения есть 1 м, что является подходящим значением для ожидаемой величины решения.
  • Строки 7-13 сообщают, что шесть итераций по Ньютону-Рафсону потребовалось для получения сходящегося решения. Первый столбец сообщает о номере итерации, а во втором приведены сведения об ошибке, используемые для определения сходимости решения. По умолчанию, критерий сходимости равен 0.001. Третья колонка показывает, что для первых двух итераций использовалась релаксация, но шаги 3-6 остались не срелаксированными.
  • Строки 14-16 сообщают о времени получения решения и требованиях к объему памяти.

Теперь вы приобрели понимание того, как решаются нелинейные статические задачи в среде COMSOL, а также знания, как интерпретировать записи в log-файле.


Комментарии (0)

Оставить комментарий
Войти | Регистрация
Загрузка...
РУБРИКАТОР БЛОГА COMSOL
РУБРИКИ
ТЕГИ