Перейти к содержанию

Контурное интегрирование

Ориентация внешнего контура и отверстия

OpenCS умеет вычислять интегралы по полигональной области без генерации сетки фибр. Реализация находится в CScore/GreenIntegrator.cs, а вызывающие методы — в MaterialArea.ContourIntegral и MaterialArea.ContourSecantStiffness.

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

Идея через теорему Грина

Для функции f(x,y) вводится антипроизводная по y:

\[ Q(x,y) = \int_0^y f(x,t)\,dt. \]

Тогда двойной интеграл по области D переводится в интеграл по границе:

\[ \iint_D f(x,y)\,dA = -\oint_{\partial D} Q(x,y)\,dx. \]

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

Выбор dx удобен для полигональных контуров: вертикальные рёбра имеют dx = 0 и не дают вклада в эту форму записи — интегратор их просто пропускает.

Усилия сечения

Для напряжения σ(x,y) OpenCS за один обход вычисляет:

\[ N = \iint_D \sigma\,dA, \]
\[ M_x = \iint_D \sigma y\,dA, \qquad M_y = \iint_D \sigma x\,dA. \]

Внутренние антипроизводные имеют вид:

\[ Q_N(x,y)=\int_0^y \sigma(x,t)\,dt, \]
\[ Q_{M_x}(x,y)=\int_0^y t\sigma(x,t)\,dt, \]
\[ Q_{M_y}(x,y)=xQ_N(x,y). \]

Последнее равенство выполняется точно, а не приближённо: множитель x в подынтегральном выражении не зависит от переменной внутреннего интегрирования t и выносится за знак интеграла. Поэтому My не требует своей квадратуры и достаётся бесплатно — GreenIntegrator.IntegrateN_Mx_My возвращает тройку N, Mx, My за один проход по контурам, вычисляя всего две антипроизводные.

Контуры и отверстия

Внешний контур должен обходиться против часовой стрелки, отверстие — по часовой. Конструктор GreenIntegrator нормализует ориентацию через EnsureCCW и EnsureCW: ориентация определяется по знаку площади полигона, и при необходимости порядок вершин разворачивается.

flowchart TB
    Outer[Внешний контур CCW] --> Integrator[GreenIntegrator]
    Hole[Отверстие CW] --> Integrator
    Integrator --> Result[Внешний интеграл минус внутренний]

Contour в OpenCS обычно хранит замыкающую вершину, совпадающую с первой. Перед передачей в GreenIntegrator последняя точка исключается, а замыкание выполняется самим интегратором.

Ориентация важна

Если отверстие задано в той же ориентации, что и внешний контур, оно прибавится вместо вычитания. Если результат площади или N имеет неправильный знак, первой проверяйте ориентацию и замыкающую точку.

Квадратура Гаусса–Лежандра

На каждом ребре параметр t пробегает [0,1]. Интеграл по ребру переводится на стандартный интервал [-1,1], после чего вычисляется таблицей узлов и весов Гаусса–Лежандра.

Есть два уровня квадратуры:

  • outerGaussN — интегрирование вдоль ребра;
  • innerGaussN — вычисление Q по координате y.

Квадратура порядка n точно интегрирует многочлен степени 2n − 1. При n = 5 это многочлен девятой степени — с запасом достаточно для гладкого участка любой диаграммы из числа применяемых в СП 63.

Порядок квадратуры не настраивается пользователем

Конструктор GreenIntegrator принимает оба порядка в диапазоне 1…10 со значением по умолчанию 5, но расчётный путь OpenCS создаёт интегратор без явных аргументов. Поэтому в рабочем расчёте всегда действуют outerGaussN = innerGaussN = 5, и «подобрать порядок» из интерфейса нельзя. Изменить их можно только в коде или в тестах.

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

Разбиение по критическим деформациям

Квадратура Гаусса построена на предположении, что подынтегральная функция гладкая. Диаграмма материала это предположение нарушает: у неё есть изломы. Узел квадратуры, поставленный по обе стороны излома, даёт систематическую ошибку, которую не убрать повышением порядка.

Разбиение ребра контура по критическим деформациям

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

\[ \varepsilon(x,y)=e_0+\kappa_y y+\kappa_z x. \]

Линейность здесь принципиальна: она делает поиск точек разбиения точным, а не итерационным. Вдоль ребра деформация меняется линейно по параметру t, поэтому положение излома находится прямым делением:

\[ t^{*} = \frac{\varepsilon^{*} - \varepsilon(t=0)}{\varepsilon(t=1) - \varepsilon(t=0)}. \]

Интегратор получает массив critEps из Diagramm.GetCriticalStrains() и находит:

  1. параметры t пересечения ε(x(t),y(t)) = ε* с каждым ребром;
  2. координаты t пересечения ε(x,t) = ε* на внутреннем отрезке 0…y;
  3. интервалы между найденными точками.

Квадратура применяется отдельно на каждом интервале. Если производная деформации по направлению меньше 1e-14, деформация вдоль отрезка считается постоянной и разбиение не выполняется — искать пересечение не с чем.

Неполный список критических деформаций

Разбиение работает ровно настолько, насколько полон critEps. Если диаграмма имеет излом, не объявленный в списке критических деформаций, интегратор о нём не узнает и проинтегрирует через него. Это единственный способ получить систематическую ошибку контурного пути на корректной геометрии.

Геометрические моменты

IntegrateMonomials вычисляет шесть интегралов произвольной функции f:

\[ A_0=\int f, \quad A_x=\int xf, \quad A_y=\int yf, \]
\[ A_{xx}=\int x^2f, \quad A_{xy}=\int xyf, \quad A_{yy}=\int y^2f. \]

При f = E_sec(ε) эти величины преобразуются в SecantStiffnessMatrix. Поэтому контурный путь используется не только для усилий, но и для точной секущей жёсткости полигональной области.

Секущий модуль определяется как σ(ε)/ε. Вблизи нуля это отношение вырождается, поэтому при |ε| < 1e-20 берётся пробное значение деформации, а если и оно не даёт ненулевого напряжения — касательный модуль диаграммы.

Ограничения

  • Контурный путь работает с полигональной областью Hull и её отверстиями.
  • Точечные фибры сами по себе в GreenIntegrator не попадают.
  • Для сложного закона материала нужно передавать корректные критические деформации.
  • Криволинейные границы должны быть аппроксимированы полигоном с достаточным числом сегментов.
  • Контурный метод не заменяет визуализацию состояния отдельных фибр.

Быстрая проверка

Для прямоугольника с постоянным σ должны выполняться:

\[ N=\sigma A, \qquad M_x=\sigma A\bar y, \qquad M_y=\sigma A\bar x. \]

Для прямоугольника, центрированного относительно начала координат, Mx = My = 0. Добавление центрального отверстия должно уменьшить N, а при симметрии оставить моменты равными нулю.

Тот же контроль зафиксирован в CScore.Tests/GreenIntegratorTests.cs. Для прямоугольника с вершинами (±1, ±2) и постоянной функцией f = 3:

Величина Аналитика Ожидание теста
A0 3 · 2 · 4 24
Ax, Ay, Axy ноль по симметрии 0
Axx 3 · b³h/12 = 3 · 2³·4/12 8
Ayy 3 · bh³/12 = 3 · 2·4³/12 32

Обратите внимание: Axx и Ayy совпадают с аналитикой точно. Фибровая сетка на той же геометрии занижает их в отношении 1 − 1/n², где n — число делений (вывод и таблица — на странице фибрового метода). Это и есть основной практический аргумент в пользу контурного пути.