Слабая форма, линейный треугольник и независимые источники погрешности

От устройства к числу

Инженер редко решает уравнения Лапласа ради них самих. Ему нужно знать ёмкость кабеля, силу притяжения якоря, температуру обмотки. Между реальным устройством и этими числами лежит цепочка решений, каждое из которых вносит свою погрешность и свои допущения. Численное моделирование — это не «нажать кнопку и получить картинку», а осознанное прохождение этой цепочки с пониманием, что именно теряется на каждом шаге.

Цепочка выглядит так:

1. Физическая система — реальный объект со всеми деталями: фасками, зазорами, неоднородными материалами, крепёжными отверстиями.
2. Концептуальная модель — решение, что учитывать: какие явления существенны (статика или гармонический режим, нужна ли связь с тепловой задачей), какие детали геометрии влияют на искомую величину, а какие нет.
3. Математическая модель — уравнение в частных производных, расчётная область, материальные уравнения и граничные условия. Это краевая задача из урока о граничных условиях и единственности.
4. Дискретная модель — сетка и система алгебраических уравнений, приближающая краевую задачу.
5. Решение — узловые значения потенциала, найденные решателем системы.
6. Постобработка — вычисление поля (градиента), а затем интегральных величин: заряда, потока, энергии, силы, потерь.
7. Интерпретация — ответ на исходный инженерный вопрос с указанием границ достоверности.

Ошибки, сделанные на шагах 2–3, никакое сгущение сетки не исправит: если выбрана неверная модель, её точное решение так же бесполезно, как и грубое. Поэтому начнём с модельных решений, а уже затем перейдём к математике дискретизации.

Концептуальная и математическая модель

Расчётная область и её усечение

Поле, строго говоря, занимает всё пространство. Расчётная область конечна, и её внешнюю границу приходится придумывать. Для электростатической задачи о двух электродах окружающий «воздушный ящик» обычно берут в несколько характерных размеров системы и задают на нём либо $\varphi=0$ (условие Дирихле), либо $\partial\varphi/\partial n=0$ (условие Неймана). Оба условия искусственны. Практический приём: решить задачу с обоими типами условий на внешней границе и увеличивать область до тех пор, пока результаты не перестанут заметно различаться. Если расхождение велико, область мала, а не сетка груба.

Упрощение геометрии

Чертёж изделия содержит детали, несущественные для поля, но дорогие для сетки: мелкие скругления, тонкие покрытия, резьбу. Решение об их удалении зависит от искомой величины. Скругление радиусом 0,1 мм не изменит ёмкость конденсатора с зазором 10 мм, но определит максимальную напряжённость на кромке электрода — для расчёта пробивного напряжения его убирать нельзя. Тонкий слой лака толщиной в десятки микрометров сетка разрешит только ценой огромного числа элементов; его заменяют поверхностным условием (контактным сопротивлением в тепловой задаче, поверхностным импедансом — в задаче о вихревых токах). Симметрия сокращает область в два, четыре и более раз: на плоскости симметрии ставится условие Неймана или Дирихле в зависимости от чётности поля относительно этой плоскости.

Материальные данные

Проницаемость, проводимость, теплопроводность берут из справочников, протоколов материалов или измерений. Все они известны с ограниченной точностью: проводимость меди меняется на 0,4 % на каждый кельвин, кривая $B$–$H$ электротехнической стали различается от партии к партии на десятки процентов в области насыщения, а коэффициент теплоотдачи $\alpha$ при естественной конвекции известен в лучшем случае с точностью 20–30 %. Погрешность исходных данных часто превышает погрешность численного метода, и бессмысленно добиваться четвёртого знака от сетки, когда третий не гарантирован материалом.

Граничные условия как модельные решения

Граничное условие — не свойство уравнения, а утверждение о том, что происходит за пределами области. Условие $\varphi=U$ на электроде означает, что электрод идеально проводящий; условие $\partial T/\partial n=0$ на поверхности — что она идеально теплоизолирована; конвективное условие $-\lambda\,\partial T/\partial n=\alpha(T-T_\infty)$ заменяет всю гидродинамику окружающего воздуха одним коэффициентом. Каждое из них — упрощение внешнего мира, и его следует проверять так же критически, как геометрию и материалы.

Метод конечных элементов: идея

Пусть требуется решить уравнение $$-\nabla\cdot (\varepsilon\,\nabla \varphi)=\rho$$ в области $\Omega$ с границей $\Gamma$. Точное решение — функция бесконечного числа «степеней свободы». Метод конечных элементов (МКЭ) заменяет её функцией, определяемой конечным набором чисел — значениями в узлах сетки.

Для линейного треугольника с вершинами $(x_1,y_1)$, $(x_2,y_2)$, $(x_3,y_3)$ и площадью $S$ функции формы линейны по координатам: $$N_i=\frac{a_i+b_i x+c_i y}{2S},\qquad a_1=x_2y_3-x_3y_2,\quad b_1=y_2-y_3,\quad c_1=x_3-x_2,$$ остальные коэффициенты получаются циклической перестановкой индексов $1\to2\to3\to1$. Удвоенная площадь $2S=(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)$ положительна при обходе вершин против часовой стрелки.

Из третьего свойства следует важное: напряжённость $\uE_h=-\nabla \varphi_h$ в линейных элементах кусочно-постоянна и терпит разрыв при переходе через рёбра. Сам потенциал $\varphi_h$ непрерывен, поскольку на общем ребре двух треугольников он определяется одними и теми же двумя узловыми значениями.

Слабая формулировка

Приближение $\varphi_h$ — ломаная поверхность; её вторые производные не существуют, и подставить её в уравнение Пуассона напрямую невозможно. Выход — заменить уравнение равносильным ему интегральным соотношением, в котором производных на порядок меньше.

Поверхностный интеграл раскрывает природу граничных условий. Пусть граница состоит из части $\Gamma_D$, где потенциал задан ($\varphi=\varphi_0$), и части $\Gamma_N$, где задана нормальная составляющая смещения ($-\varepsilon\,\partial\varphi/\partial n=D_n$).

• На $\Gamma_D$ значение решения известно заранее, поэтому пробные функции выбирают равными нулю на $\Gamma_D$, и соответствующая часть интеграла исчезает. Такое условие называют главным (существенным): его накладывают явно на узловые значения.
• На $\Gamma_N$ в интеграл подставляется известная величина $D_n$. Такое условие называют естественным: оно входит в правую часть само, без каких-либо ограничений на неизвестные. Однородное условие Неймана $\partial\varphi/\partial n=0$ вообще не требует никаких действий — оно выполняется «по умолчанию» на любой границе, где ничего не задано.

Это объясняет известный практический факт: свободная граница сетки ведёт себя как линия симметрии поля (нормальная составляющая $\uD$ равна нулю).

Метод Галёркина для уравнения Пуассона равносилен минимизации энергетического функционала $\tfrac12\int\varepsilon|\nabla \varphi_h|^2\dV-\int\rho\varphi_h\dV$ на множестве функций вида $\sum u_jN_j$ [Сильвестер]. Иными словами, МКЭ находит среди всех ломаных поверхностей ту, которая ближе всех к точному решению по энергии. Этот факт будет ключом к пониманию скоростей сходимости в уроке о сходимости.

Элементная матрица и сборка

Интеграл $K_{ij}$ распадается на сумму интегралов по элементам, причём в элемент $e$ вклад дают только три его узла. Поскольку градиенты функций формы постоянны, интегрирование сводится к умножению на площадь: $$K^{e}_{ij}=\int_{e}\varepsilon\,\nabla N_i\cdot\nabla N_j\,\dS =\varepsilon\,S\,\frac{b_ib_j+c_ic_j}{(2S)^2} =\frac{\varepsilon\,(b_ib_j+c_ic_j)}{4S},\qquad i,j=1,2,3 .$$ Вектор нагрузки от постоянной в элементе плотности заряда распределяется между вершинами поровну: $f^e_i=\int_e\rho N_i\,\dS=\rho S/3$, так как интеграл линейной функции по треугольнику равен произведению площади на её среднее значение, а среднее значение $N_i$ по треугольнику равно $1/3$.

Сборка. Глобальная матрица получается суммированием элементных: вклад $K^e_{ij}$ прибавляется к элементу глобальной матрицы на пересечении строки и столбца, соответствующих глобальным номерам узлов $i$ и $j$ элемента. Поскольку каждый узел связан лишь с несколькими соседями, в строке глобальной матрицы отлично от нуля лишь 5–10 элементов из тысяч — матрица разреженная, и хранятся только ненулевые элементы.

Наложение условий Дирихле. В собранной системе узлы с заданным потенциалом исключают: соответствующие неизвестные полагают равными заданным значениям, их столбцы переносят в правую часть, а строки вычёркивают. Если бы условий Дирихле не было вовсе, матрица оказалась бы вырожденной (сумма строк равна нулю, как в примере): потенциал определён с точностью до постоянной, и без опорного потенциала задача не имеет единственного решения — численное выражение теоремы единственности.

Решение системы

Полученная матрица симметрична и положительно определена. Для её решения применяют два класса методов.

Прямые методы (разложение Холецкого, $LU$-разложение) дают решение за конечное число операций с точностью до ошибок округления. Их недостаток — заполнение: в ходе разложения нули разреженной матрицы превращаются в ненулевые элементы, и для больших трёхмерных задач памяти не хватает. Для плоских задач с сотнями тысяч узлов прямые методы остаются предпочтительными.

Итерационные методы строят последовательность приближений, сходящуюся к решению. Для симметричных положительно определённых матриц основным является метод сопряжённых градиентов: на каждом шаге требуется лишь умножение матрицы на вектор, память не растёт, а точность задаётся порогом невязки $\|\mathbf K\mathbf u-\mathbf f\|$. Скорость сходимости зависит от обусловленности матрицы, которая ухудшается при сгущении сетки и при сильном контрасте материалов (сталь — воздух), поэтому метод применяют с предобусловливанием.

Решатель находит решение дискретной системы с заданной точностью — скажем, $10^{-8}$ по невязке. Это не означает, что решение краевой задачи получено с такой точностью: погрешность дискретизации, то есть разница между $\varphi_h$ и $\varphi$, обычно на много порядков больше и не зависит от того, насколько точно решена система.

Нелинейные задачи

В магнитостатике с ферромагнетиком проницаемость $\mu$ зависит от индукции $B$, которая сама подлежит определению. Матрица становится функцией решения: $\mathbf K(\mathbf u)\mathbf u=\mathbf f$. Два основных подхода:

• Метод простых итераций (Пикара): задаются начальным $\mu$, решают линейную задачу, по найденному $B$ в каждом элементе определяют новое $\mu(B)$ из кривой намагничивания, собирают матрицу заново и повторяют до стабилизации. Сходится надёжно, но медленно, особенно при глубоком насыщении; часто требует релаксации — усреднения новых и предыдущих значений проницаемости в элементах.
• Метод Ньютона: уравнение линеаризуется относительно поправки с учётом производной $\dd\mu/\dd B$ (дифференциальной проницаемости). Сходимость квадратичная вблизи решения — число верных знаков удваивается на каждом шаге, — но при неудачном начальном приближении метод может расходиться; на практике первые шаги делают по Пикару, затем переключаются на Ньютона.

Численная модель магнитной цепи в лаборатории урока о магнитных цепях использует метод Ньютона с касательной проницаемостью (простые итерации для насыщенной стали сходятся медленно или расходятся), а лаборатория сходимости — линейные задачи без итераций.

Численное интегрирование

Для линейных треугольников интегралы в $K^e_{ij}$ вычисляются точно. При переменных в элементе свойствах ($\varepsilon(\ur)$, $\lambda(T)$), криволинейных элементах или элементах высшего порядка подынтегральные функции перестают быть полиномами низкой степени, и интегралы берут численно — квадратурами Гаусса: значение функции вычисляют в нескольких специально выбранных точках элемента и суммируют с весами. Квадратура из $m$ точек на отрезке точна для полиномов степени $2m-1$; недостаточное число точек — ещё один, обычно малый, источник погрешности.

Постобработка

Решение системы — набор узловых потенциалов. Всё остальное получается из него дифференцированием и интегрированием, и здесь кусочная природа $\varphi_h$ проявляется в полной мере.

Локальное поле. В каждом линейном элементе $\uE_h=-\sum u_i\nabla N_i$ — постоянный вектор. На ребре между элементами касательная составляющая $\uE_h$ непрерывна (потенциал на ребре одинаков с обеих сторон), а нормальная составляющая $\uD_h$ терпит скачок, хотя в точной задаче при отсутствии поверхностного заряда она непрерывна. Величина скачка — непосредственная мера локальной погрешности. Поле в узле обычно получают усреднением по примыкающим элементам; пространственная «ступенчатость» карты $|\uE|$ на грубой сетке — не физика, а артефакт аппроксимации.

Интегральные величины оказываются существенно точнее локального поля. Энергия $W=\tfrac12\int\varepsilon|\uE_h|^2\dV=\tfrac12\mathbf u^{\mathsf T}\mathbf K\mathbf u$ вычисляется из той же квадратичной формы, которую минимизирует метод Галёркина, и локальные ошибки поля в ней частично компенсируются. Ёмкость $C=2W/U^2$, индуктивность $L=2W/I^2$, потери $P=\int\sigma|\uE_h|^2\dV$ наследуют эту точность. Напротив, заряд, вычисленный как поток $\uD_h$ через контур вокруг электрода, зависит от выбора контура и использует разрывное поле напрямую, поэтому сходится медленнее. Количественно эти различия разобраны в уроке о сходимости.

Распространённые заблуждения

Итоги

• Численное моделирование — цепочка: физическая система → концептуальная модель → математическая модель → дискретная модель → решение → постобработка → интерпретация. Ошибки модели не исправляются сеткой.
• Усечение области, упрощение геометрии, материальные данные и граничные условия — модельные решения, каждое со своей погрешностью; граничное условие описывает то, что находится вне области.
• МКЭ ищет решение в виде $\sum u_iN_i$; для линейного треугольника $N_i=(a_i+b_ix+c_iy)/(2S)$, функции формы обладают свойством Кронекера и образуют разбиение единицы, их градиент постоянен в элементе.
• Слабая форма получается умножением на пробную функцию и интегрированием по частям; условия Дирихле — главные, Неймана — естественные. Метод Галёркина даёт систему $\mathbf K\mathbf u=\mathbf f$ с $K^e_{ij}=\varepsilon(b_ib_j+c_ic_j)/(4S)$, равносильную минимизации энергии.
• Решатель точно решает дискретную систему; погрешность дискретизации остаётся. Поле кусочно-постоянно и разрывно; интегральные величины точнее локальных.

Далее — размерность моделей: как трансляционная и вращательная симметрия сводят трёхмерную задачу к плоской или осесимметричной и как правильно пересчитывать результаты на реальное устройство.