Проблемы с переходными нелинейными задачи. Часть 2: поиск источника проблем со сходимостью
В данной задаче в какой-то момент Comsol начинает сильно мельчить шаг по времени. Возникает два вопроса - в чем причина такого «поведения» и как ускорить решение?
- Может быть, тройная точка дает эти проблемы…
- Но когда вы решаете электростатику, у вас тоже встречаются в моделях тройные точки, тем не менее, те задачи сходятся хорошо.
- Можно начальный шаг попробовать поварьировать…
- Может быть, стоит измельчить сетку.
- Можно попробовать ограничить шаг по времени. Но с другой стороны, решение, наверное, тогда вообще перестанет сходиться…
- В матрице задачи очень разные по порядку величины значения - может быть, в этом дело? Вблизи электрода большая проводимость, а вдали от него очень маленькая. И матрицу из-за этого, наверное, тяжело обращать. Может быть, в скэйлинг-факторах поменять значения?
- Смотрите, когда возникает подобная проблема, обычно, как у нас здесь сейчас, появляется много идей - какие настройки можно покрутить. Потому что даже когда решается всего одно уравнение, настроек очень много. Интересно понять, как в такой ситуации действовать разумно и выявить проблемную точку, не перебирая все комбинации настроек. Поэтому попробуем найти причину этого замедления.
Есть простой способ такого анализа - стартовать с простой постановки, которая решается без проблем, и постепенно ее усложнять, обнаружив, добавление какого элемента создает проблему. Например, в нашем случае можно было бы начать с электростатической задачи без проводимости вообще (вы хорошо знаете из опыта, что такие задачи решаются «хорошо»), затем добавить фиксированную небольшую проводимость материала, затем ее увеличить и, наконец, задать зависимость проводимости от напряженности поля. Вы обнаружите, что именно на последнем шаге возникают проблемы со сходимостью.
Рис. 5. Зависимость обратного шага по времени от времени для задачи с фиксированной проводимостью и для задачи с проводимостью, зависящей от напряженности поля.
Можно обнаружить проблему и другим путем - обратить внимание на то, что первое время задача решается быстро, а момент, когда шаги по времени начинают уменьшаться соответствует моменту, когда максимум напряженности в изоляции превышает E0. Т.е. пока во всех точках изоляции напряженность меньше E0 и проводимость материала везде одинакова, задача решается быстро. Однако как только в какой-то области «начинает работать» правая часть зависимости σ(E), шаг по времени уменьшается.
Таким образом, можно «локализовать» проблему - найти, какой именно компонент модели «создает проблемы». Теперь хочется понять, что с этим делать.
Давайте сначала немного вспомним о том, как Comsol решает переходные задачи.
При помощи метода конечных элементов система уравнений в частных производных превращается в систему обыкновенных алгебраических уравнений. А что, если мы решаем переходную задачу?
Предположим, у нас всего один узел, и мы имеем обыкновенное дифференциальное уравнение:
(7)
Для решения таких уравнений используются схемы численного интегрирования. Простейшую из схем можно получить, заменив точную производную ее приближенным значением:
Тогда уравнение (7) превращается из дифференциального в обыкновенное:
В правой части стоит g(u,t). Здесь u - это u(t) или u(t+Δt)? Два этих варианта дают две разные схемы. В первом случае мы получаем:
Если мы знаем значение u(t) - на предыдущем шаге, мы легко находим значение на следующем шаге u(t+Δt):
Такие схемы - в которых значение искомой функции на следующем шаге явно выражается через значения функции на предыдущих, уже известных шагах - называются явными.
Во втором же случае мы получаем следующую схему:
(8)
Это уравнение уже не разрешается в общем виде относительно u(t+Δt). Явно разрешить его можно только если функция g(u,t) имеет совсем простой вид. Такие схемы называют неявными - они приводят к обыкновенным алгебраическим уравнениям.
Казалось бы, можно упростить себе жизнь и решать уравнения при помощи явных схем. Однако они обладают плохой точностью, и использовать их нежелательно.
Как решать уравнения, к которым приводят неявные методы? Для этого можно использовать тот же метод Ньютона, уже встречавшийся нам при разговоре о стационарных задачах в Comsol.
Напомню, что метод Ньютона - итеративный метод. Таким образом, решение идет как бы на двух уровнях: с одной стороны, есть дискретные шаги по времени, и решатель двигается шаг за шагом. С другой стороны, на каждом шаге по времени необходимо запускать метод Ньютона и итеративно искать решение - такие итерации будем называть «внутренними», имея в виду, что эти итерации происходят внутри одного шага по времени.
Рис. 6. Решение нелинейных переходных задач при использовании неявных схем численного интегрирования по времени движется как бы на двух уровнях: внутри каждого шага по времени выполняется необходимое число внутренних итераций.
Осторожно! Скучные выкладки:
Напомню в основных чертах схему решения при использовании метода Ньютона, на конкретном примере одного шага переходной задачи. Итак, на одном шаге переходной задачи нам нужно решить уравнение (8):
(8)
Значение u(t) нам известно, это значение искомой функции на предыдущем шаге, обозначим его u0=u(t) и будем также использовать как начальное значение (догадку) для метода Ньютона. Таким образом, получаем уравнение:
(9)
Искомая величина здесь: u(t+Δt). Для ее поиска мы производим итерации, стартуя с u0: u0, u1, u2, …. , uj → u(t+Δt). Сколько итераций понадобится, заранее предсказать нельзя, мы производим их, пока не будет выполнен критерий сходимости.
Уравнение (9) все еще неразрешимо в общем случае, поскольку в аргументе функции g стоит неизвестная величина - искомая функция u(t+Δt). В методе Ньютона для решения этой проблемы в правой части уравнения делаем замену (10):
(10)
Если бы приближенное равенство (10) выполнялось точно, мы бы решили уравнение (9) за одну итерацию. В общем случае, из-за того, что замена (10) неверна, нам необходимо провести более одной итерации. А именно, раз за разом решать уравнение:
Отсюда можно выразить значение оценки искомой функции на следующей итерации uj+1 через предыдущую итерацию uj:
(11)
Проблемы у этого итеративного метода могут возникнуть, когда знаменатель (11) обращается в ноль, а именно, если:
В этом случае метод Ньютона не сходится, итерации, вместо того чтобы давать все более близкие друг к другу значения uj, напротив, дают все более отстоящие друг от друга значения. Поэтому для больших Δt метод Ньютона не срабатывает. Чтобы избежать подобной проблемы, необходимо ограничивать шаг по времени, следить за тем, чтобы:
(12)
Если шаг по времени выбирается пользователем, сам пользователь и должен следить за выполнением условия (12). В Comsol длительность шага по времени Δt выбирается программой: если итерации метода Ньютона на данном шаге не сходятся, либо сходятся очень медленно, Comsol уменьшает шаг вдвое. Если сходимость, напротив, достигается за очень малое число итераций (например, две или три), Comsol «понимает», что можно увеличить шаг по времени и увеличивает его вдвое.
Отсюда видно, что вызывает замедление решения: если ∂g/∂u(uj,t+Δt) становится большим, Comsol вынужден уменьшать шаг по времени Δt, чтобы выполнялось условие (12).
Вернемся теперь к нашей задаче. Попытаемся провести аналогию и понять, какой порядок шага по времени можно в ней ожидать.
Осторожно! Скучные выкладки:
Для понимания мы рассматривали простое уравнение:
(7)
И пришли к выводу, что для такого уравнения Comsol будет держать шаг по времени порядка:
Мы же решаем уравнение в частных производных:
(6)
Попробуем путем грубых оценок по аналогии получить оценку для шага по времени для системы (6). По времени интегрируется первое уравнение в системе (6). Поэтому сопоставляя его с (7), u↔ρ, а g↔div(σ(E)E). Тогда можно ожидать, что:
Из второго уравнения системы (6) мы видим, что:
Также нам поможет следующее соотношение:
Получаем:
(13)
Оценим теперь по порядку величины каждое из слагаемых:
Здесь λ - пространственный масштаб изменения электрического поля.
Поскольку в скобках в (13) стоит сумма двух слагаемых, оценивая их по порядку величины, мы должны взять максимальное.
(14)
Таким образом, в зависимости от свойств материала оценкой для шага по времени, который будет выбран пакетом конечноэлементного моделирования, будет либо [σ/εε0]-1 (т.е. время Максвелловской релаксации, о котором речь шла выше), либо [dσ/dE•E/εε0]-1. Видно, что в случае сильной зависимости проводимости от напряженности поля происходит дополнительное замедление решения по сравнению с временем максвелловской релаксации.
Приведенные здесь выкладки могут показаться достаточно сложными. На практике можно получить полезную оценку гораздо проще, для этого заменим в системе дифференицальных уравнений в частных производных дифференциалы величин на сами величины, при этом dx→λ (пространственный масштаб изменения величины), а dt→τ (шаг по времени). Тогда из системы (6) получаем:
Отсюда получаем:
(15)
Конечно, более аккуратные выкладки привели нас к более содержательному результату (14) (в оценке (15) нет второй части оценки (14), связанной с производной dσ/dE). Однако и оценка (15) является полезной и дает близкую по порядку величины оценку шага по времени: действительно, рис. 7 показывает, что уменьшение шага по времени происходит синхронно с увеличением максимальной проводимости. А это означает, что простая оценка тоже позволяет выявить источник проблемы со сходимостью, и понять, что конкретно нужно поменять в модели, чтобы с этой проблемой справиться.
Рис. 7. Корреляция между шагом по времени и максимальной проводимостью в модели.
Итак, проведенные рассуждения позволяют по виду уравнения предсказать, какой можно ожидать шаг решения при расчете переходной задачи.