Независимое уточнение сетки, области и положения

Чему верить в численном решении

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

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

Погрешность дискретизации и порядок сходимости

Порядок $p$ определяется двумя обстоятельствами: степенью полиномов в функциях формы и гладкостью точного решения. Для линейных треугольников и гладкого решения теория даёт следующие априорные ожидания [Сегерлинд; Сильвестер].

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

Ёмкость по энергии и по заряду

Ёмкость можно вычислить и иначе — через заряд электрода, то есть поток смещения через замкнутый контур вокруг него: $C_Q=Q_h/U$, где $Q_h=\oint\uD_h\cdot\un\,\dd l$ (на единицу глубины). Здесь используется само разрывное поле $\uD_h$, а не его квадрат, и условие ортогональности не помогает. Результат зависит от того, где проведён контур: через элементы, прилегающие к электроду, или дальше — в точной задаче поток один и тот же, в дискретной — нет, ибо $\nabla\cdot \uD_h$ отличен от нуля на рёбрах. Поэтому $C_Q$, как правило, сходится медленнее, чем $C_W$, может приближаться к пределу с любой стороны, а разность $|C_Q-C_W|$ служит внутренней оценкой погрешности, не требующей знания точного ответа. Обе величины сходятся к одной ёмкости — их расхождение на конечной сетке есть мера того, насколько дискретное решение ещё не удовлетворяет закону Гаусса локально.

Оценка порядка по нескольким сеткам

Если точное решение известно, порядок находят по двум сеткам: $$p=\frac{\ln(e_1/e_2)}{\ln(h_1/h_2)} .$$ При $e_1=4\,\%$, $e_2=1\,\%$ и $h_1/h_2=2$ получаем $p=\ln4/\ln2=2$.

Точное решение обычно неизвестно — иначе и считать было бы незачем. Тогда применяют экстраполяцию Ричардсона [Самарский].

Экстраполяция Ричардсона применима только в асимптотической области, где старший член $Ch^p$ действительно доминирует. Признак того, что область достигнута, — стабильность наблюдаемого $p$ при переходе к следующей тройке сеток. На графике погрешности от $h$ в двойном логарифмическом масштабе асимптотика видна как прямая с наклоном $p$; загиб на грубых сетках означает, что асимптотика не достигнута.

Проверки баланса

Дискретное решение должно удовлетворять законам сохранения: полный поток $\uD$ через замкнутую поверхность вокруг всех электродов равен их суммарному заряду; сумма токов, втекающих в проводник, равна нулю; в стационарной тепловой задаче подведённое тепло равно отведённому плюс выделенному. Эти проверки обязательны, и их нарушение указывает на ошибку: неверное граничное условие, незамкнутую область, неточно решённую систему.

Однако выполнение баланса не гарантирует точности. Если потоки вычисляются согласованно с дискретной формулировкой — через те же уравнения, что и решение, — то метод Галёркина обеспечивает баланс с точностью решателя на любой, даже грубой, сетке: сумма уравнений Галёркина по всем узлам электрода есть не что иное, как равенство его полного потока заряду. Поток, найденный прямым интегрированием разрывного поля $\uD_h$ по контуру, удовлетворяет балансу лишь приближённо, и степень этого приближения также слабо связана с погрешностью ёмкости. Баланс выполняется на сетке из трёх элементов — ёмкость при этом может быть неверна вдвое. Баланс — необходимое условие, сгущение сетки — путь к достаточному.

Особые точки

Априорные порядки из таблицы предполагают гладкое решение. Вблизи острых кромок проводников, входящих углов и точек смены граничного условия решение уравнения Лапласа негладко, и сходимость ухудшается.

Следствия для численного решения таковы.

• Локальное поле у кромки не сходится вовсе. Какую бы сетку ни взять, вычисленное $|\uE|$ в элементе у вершины конечно и растёт при сгущении примерно как $h^{-1/3}$. Величина «максимальная напряжённость» для геометрически острой кромки не имеет смысла — реальные кромки скруглены, и предельное поле определяется радиусом скругления, который следует включить в модель, если он важен.
• Интегральные величины сходятся, но медленнее. На равномерной сетке погрешность энергии и ёмкости убывает как $O(h^{2\pi/\theta})$ вместо $O(h^2)$; для кромки с $\theta=3\pi/2$ это $O(h^{4/3})$. Наблюдаемый порядок $1{,}3$–$1{,}4$ на задаче с углом — не ошибка, а признак особенности.
• Локальное сгущение восстанавливает порядок. Если размер элементов убывает к вершине по геометрической прогрессии (градиентная сетка), погрешность интегральных величин снова убывает как $O(h^2)$ при том же общем числе узлов. Это главный аргумент в пользу неравномерных сеток: сгущать нужно там, где решение негладко, а не везде.

Ещё один источник погрешности, не связанный с $h$ напрямую, — численное интегрирование. Для линейных элементов с постоянными свойствами интегралы берутся точно; при переменных свойствах и в осесимметричных задачах с весом $2\pi r$ квадратуры Гаусса вносят погрешность, которая обычно много меньше погрешности дискретизации, но может стать заметной для элементов, примыкающих к оси.

Лаборатория

Лаборатория решает уравнение Лапласа методом конечных элементов с линейными треугольниками для двух задач. Первая — четверть коаксиального конденсатора ($a=10$ мм, $b=30$ мм, $U=100$ В): на дугах заданы потенциалы $U$ и $0$, на двух радиальных отрезках — условия симметрии. Все величины погонные и относятся к четверти сечения, поэтому аналитический эталон $C'=2\pi\varepsilon_0/\ln(b/a)=50{,}64$ пФ/м берётся с множителем $1/4$: $12{,}66$ пФ/м; поле в контрольной точке $r=a+(b-a)/4=15$ мм равно $|\uE|=U/(r\ln(b/a))=6068$ В/м. Вторая — электрод с острым углом (квадрат $[0;\tfrac12]^2$ в единичной области): со стороны диэлектрика угол входящий, решение имеет особенность вида $r^{\pi/\theta-1}$ (для прямоугольной кромки — $r^{-1/3}$), аналитического эталона в замкнутой форме нет, и за эталон ёмкости принято значение, экстраполированное по Ричардсону.

Четыре уровня сгущения сетки $h$, $h/2$, $h/4$, $h/8$ (на нулевом уровне — 4 или 8 ячеек на сторону). Для каждого уровня показаны ёмкость по энергии $C_W=2W_h/U^2$, ёмкость по заряду $C_Q=Q_h/U$ (поток $\uD_h$ через поверхность внутреннего электрода), $|\uE_h|$ в контрольной точке — для коаксиала в узле, где поле восстановлено усреднением по прилегающим элементам, для кромки в элементе у вершины, — погрешности относительно эталона, а также график погрешности от $h$ в двойном логарифмическом масштабе с оценкой наблюдаемого порядка. Карта потенциала и сетка показаны для каждого уровня.

На что обратить внимание:

• Для коаксиальной задачи наклон прямой погрешности $C_W$ равен $2$, и $C_W$ приближается к эталону сверху. Погрешность $C_Q$ убывает лишь с порядком около $1$ (в этой задаче — снизу): поток разрывного поля через поверхность электрода наследует локальную погрешность поля.
• Поле в узле, восстановленное усреднением, сходится быстрее, чем предсказывает оценка $O(h)$ для поля внутри элемента: наблюдаемый порядок близок к $2$. На регулярной сетке погрешности соседних элементов имеют разные знаки и при усреднении частично гасятся (суперсходимость в узлах); на нерегулярной сетке и вблизи границ на это рассчитывать нельзя.
• Сравните $C_Q$ и $C_W$ на каждом уровне. Как изменяется их разность при сгущении и какая из величин ближе к эталону?
• Для задачи с кромкой $|\uE|$ в элементе у вершины не стабилизируется, а растёт от уровня к уровню примерно в $2^{1/3}\approx1{,}26$ раза, как и предписывает особенность $r^{-1/3}$. Ёмкость $C_W$ при этом сходится, но наблюдаемый порядок около $1{,}4$ — ниже $2$ и близок к теоретическому $4/3$.
• По первым трём уровням выполните экстраполяцию Ричардсона вручную и сравните с четвёртым уровнем (для коаксиала — и с эталоном).

Сравнение с аналитическими решениями

Аналитическое решение — самый сильный эталон, потому что его погрешность равна нулю. Задач с такими решениями мало, но для каждого типа поля в курсе есть по крайней мере одна.

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

Верификация и валидация

Два слова, которые часто путают, обозначают разные и одинаково необходимые виды проверки [Oberkampf].

Структура полной погрешности результата: $$\text{результат}-\text{истина}=\underbrace{(\text{модель}-\text{истина})}_{\text{погрешность модели}}+\underbrace{(\text{точное решение модели}-\text{дискретное})}_{\text{дискретизация}}+\underbrace{(\text{округление, решатель})}_{\text{обычно пренебрежимо}} .$$ Первое слагаемое включает и неопределённость исходных данных: если теплопроводность изоляции известна с точностью $\pm10\,\%$, расчётная температура неопределённа примерно в той же мере, и сгущать сетку до $0{,}1\,\%$ бессмысленно. Разумная цель верификации — довести численную погрешность до уровня, заметно меньшего неопределённости данных, и не далее. Оценку чувствительности результата к данным получают, варьируя свойства в пределах их разброса, — это часть анализа неопределённости, а не сходимости.

Итоги курса

Этот урок завершает курс, и его итоги — итоги всего пройденного пути.

• Математический язык (М1). Поля описываются градиентом, дивергенцией и ротором; интегральные теоремы Остроградского–Гаусса и Стокса связывают локальные и интегральные формулировки; метрические множители криволинейных координат — не формальность, а размерность.
• Уравнения поля (М2). Уравнения Максвелла и уравнение теплопроводности в статическом и квазистатическом приближениях сводятся к краевым задачам для потенциалов; материальные уравнения и условия на границах раздела замыкают постановку; единственность требует опорного потенциала.
• Электростатика, ток, магнитостатика, гармонические поля (М3–М6). Одна и та же структура — уравнение Пуассона или Гельмгольца с коэффициентом $\varepsilon$, $\sigma$, $1/\mu$ или $\lambda$ — описывает поля разной природы; различие в источниках, граничных условиях и в том, какие интегральные величины из них извлекают.
• Интегральные величины (М7). Ёмкость, индуктивность, сопротивление, энергия и сила — функционалы поля с явно указанной областью интегрирования, ориентацией и допущениями; энергетические определения точнее потоковых.
• Тепловые поля (М8). Закон Фурье и баланс энергии дают уравнение теплопроводности; граничные условия теплообмена — модели внешнего мира; нестационарность вносит характерное время.
• Связанные задачи и численное моделирование (М9). Потери — источники тепла, свойства зависят от температуры; метод конечных элементов строится на слабой форме и функциях формы; симметрия понижает размерность ценой явных допущений; сходимость, порядок, особые точки, верификация и валидация — инструменты, отличающие расчёт от картинки.

Главный урок курса состоит в следующем: численная модель отвечает ровно на тот вопрос, который ей задан — с той геометрией, теми материалами, теми граничными условиями и той сеткой, которые выбраны. Понимать поле означает понимать, какой вопрос был задан, и уметь оценить, насколько ответ достоверен.