Введение. Ранее мы достаточно подробно рассмотрели малоизученный в численной реализации диаграммного метода расчёта железобетонных стержневых элементов вопрос о его точности (погрешности) [1]. Тесно с этим вопросом связано такое понятие, как сходимость. К сожалению, не нашлось ни одной работы как в отечественной литературе, так и зарубежной, посвящённой теории сходимости рассматриваемого метода. Цель данной публикация – в определённой мере восполнить этот пробел.

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

В статье [2] отмечается, что установить сходимость и оценить быстроту сходимости итерационного процесса не представляется возможным из-за немонотонности итерационного расчета диаграммным методом. Скорость сходимости итераций зависит от уровня нагружения. На начальных этапах итерационный процесс в пределах шага нагружения сходится в среднем за 10-15 итераций, по мере приближения к предельному состоянию скорость сходимости замедляется. На точность результата большое влияние оказывает формулировка условия прекращения итерационного процесса указанием величины минимального изменения итерационно уточняемых переменных и предельного числа итераций.

В монографии [3] расчёт статически неопределимых железобетонных конструкций с учётом физической нелинейности предложено вести методом последовательного уточнения жесткостей. Отмечается, что начальное приближение решения не влияет на сходимость итерационного метода, поэтому для формулирования критерия сходимости рекуррентная зависимость для нахождения искомого параметра (усилия в сечении стержня с учётом перераспределения) принята в виде простейшей цепной дроби вида . Согласно теореме А.Я. Хинчина для сходимости цепной дроби необходимо и достаточно, чтобы ряд  был расходящимся. Сформулировать условие сходимости процесса последовательного уточнения жесткостей удаётся лишь при рассмотрении простейших систем. Попытки сформулировать это условие для более сложных систем привели к усложнению и критериев сходимости, и аппарата, необходимого для их анализа [4].

В МКЭ, например, в [5] для итерационного уточнения на каждом шаге нагружения напряжениями используется метод переменных параметров упругости (модифицированный алгоритм Ньютона-Рафсона, в котором касательные модули заменяются на секущие). В качестве критерия сходимости применяется эвклидова норма вектора деформаций.

В Методическом пособии [6] нелинейные задачи железобетона решаются в самом общем виде в объёмной постановке методами упругих решений, переменных параметров, начальных напряжений и др. Интерес касаемо темы статьи представляет небольшой параграф этого пособия 11.7, в котором приведены общие соображения о критерии сходимости метода упругих решений. Критерий записывается через перемещения в трёх возможных вариантах:

где  – норма невязки; ε – константа (10^-2≤ε≤10^-6); 1≤i≤N – номер компоненты решения; Δr_i – разность компонент вычисления перемещений между текущей и предыдущей итерациями;  – эталонная величина перемещений.

Если процесс сходится по одной норме, то он сходится и по любой другой норме. Выбор нормы определяет скорость сходимости. При прочих равных условиях наибольшую скорость сходимости дает Евклидова норма, наименьшую – Чебышева норма.

Использовать вышеприведённые критерии для численного диаграммного метода также возможно, если вместо перемещений r_i использовать напряжения σ_i, деформации ε_i, предельный момент M_ult, момент трещинообразования M_crc, жёсткость D и др. Но более удобно, на наш взгляд, для железобетонных элементов, испытывающих изгиб (продольный либо поперечный), критерий сходимости строить через кривизну χ. Поскольку кривизна – это интегральный параметр, описывающий деформированное состояние сечения в целом, то на её основе может быть получен критерий сходимости по Чебышеву:

где  – кривизна оси стержня в расчётном сечении соответственно на текущей и предшествующей итерациях расчёта.

Ниже будет приведён строгий математический вывод этой формулы.

Из своего опыта построения алгоритмов расчёта НДС всевозможных железобетонных элементов численным диаграммным методом [7] сходимость для него может быть рассмотрена в двух формах:

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

– как сходимость дискретизации нелинейной деформационной модели нормального сечения при выполнении процедуры численного интегрирования (внутренний цикл алгоритма) нормальных напряжений, σ, секущих модулей, , жесткости по всем компонентам расчётного сечения: под сходимостью в таком случае следует подразумевать стремление значений решения дискретной модели к соответствующим значениям решения исходной задачи при стремлении к нулю параметра дискретизации (размера малых областей, Δ_h×Δ_b, на которые разбивается сечение).

В статье [8] сходимость численного интегрирования методом Рунге-Кутта определяется условием , где R – функция линейной устойчивости (например, полином 4-го порядка), Δ_h – шаг разбивки, λ – параметр, определяющий связь между функцией и её производной. В монографии [9] рекомендуется для увеличения скорости сходимости использовать формулу Симпсона (правило одной трети). Вопросам сходимости численных методов посвящены также и другие актуальные зарубежные работы [10–15].

Методика. Исследуем последовательно вначале сходимость итерационного процесса вычислений внешнего цикла алгоритма численного диаграммного метода, а затем перейдём ко внутреннему циклу и вопросу дискретизации расчётной модели.

В работе [7] для стержневых элементов, в которых помимо прочих деформаций возникает изгиб, в самом общем виде постановка задачи для решения её численным диаграммным методом состоит в том, чтобы решить уравнение вида:

с корнем t в интервале [a;b]. Здесь χ – кривизна оси стержня. При этом в качестве границ интервала, как показывает опыт расчётов, можно принять значения a=-0,1, b=0,1 практически для любых железобетонных элементов. Функция g предполагается непрерывной на этом интервале.

Отметим, что при центральном растяжении либо сжатии χ=0. В таком случае в качестве управляющего параметра итерационного алгоритма можно принять осевую жёсткость сечения D₁₁=ν_bE_bA_b+ν_sE_sA_s.

Вообще говоря, при решении задач, связанных с изгибом железобетонных элементов, в качестве управляющего параметра итерационного алгоритма численного диаграммного метода, вместо кривизны χ, могут быть использованы и другие, связанные с ней параметры: предельный момент – M_ult (либо момент трещинообразования – M_crc), относительные деформации в характерных точках сечения – , , , , , изгибная жёсткость сечения D₃₃, координата нейтральной линии y₀ и др. При этом, как показали собственные результаты расчётов, закономерности по сходимости алгоритма (приближения управляющего параметра к точному (аналитическому) своему значению) будут качественно идентичными. Причём  и т.д.

Результаты. Решим уравнение (3) численным методом простой итерации. Так, если известен какой-либо член последовательности χ_k, например, χ₀ [a;b], то χ_k+1 можно взять . Здесь k=0,1,2…m – соответственно номер текущей итерации и общее количество итераций. Тогда рекуррентная формула метода имеет вид:

Если существует конечный предел  и функция g непрерывна в точке z, переходом к пределу в равенстве (4) получим , то есть число z является корнем уравнения (1). Если z [a;b], то в силу единственности корня на отрезке [a;b] z совпадает с t.

Вычисления по формуле (4) проиллюстрированы на рис. 1.

Рис. 1. К численному методу простой итерации для отыскания кривизны

Построим графики функций из левой и правой частей уравнения (3), то есть линии y=χ и y=g(χ). Они должны пересекаться в точке с абсциссой t. Взяв некоторое число χ₀, вычислим g(χ₀) и получим на кривой y=g(χ) точку А₀. Линия проекции этой точки на ось Оу пересечёт прямую y=χ в точке В₁. Проекция В₁ на ось Оχ даёт χ₁. Из равенства треугольников ΔОВ₁χ₁ и ΔОВ₁g(χ₀) геометрически χ₁=g(χ₀). Проекция χ₁ на кривую y=g(χ) даёт точку А₁. Линия проекции этой точки на ось Оу пересечёт прямую y=χ в точке В₂. Проекция В₂ на ось Оχ даёт χ₂. Из равенства треугольников ΔОВ₂χ₂ и ΔОВ₂g(χ₁) геометрически χ₂=g(χ₁). Через какое-то количество итераций m величина χ_i = χ_m настолько близко подойдёт к t, что её можно буде считать ответом. Это количество шагов будет определять точность приближения >0 [1].

Далее запишем условия сходимости итерационного процесса решения. Пусть корень t уравнения (3) отделён на отрезке [a;b] длины h. Если на отрезке [c;d] = [a-h;b+h] функция g дифференцируема и найдётся число 0<q<1 такое что

при всех χ [c;d], то итерационная последовательность, предложенная формулой (4), сходится к корню t при любом выборе начального приближения χ₀ [a;b] (рис. 2).

Рис. 2. К вопросу сходимости решения

Геометрически условие (5) означает, что угол наклона касательной к кривой g в любой точке из интервала [c;d] должен быть меньше 45⁰ (tg45⁰=1), то есть меньше угла наклона прямой y=χ (а он равен ровно 45⁰). В этом случае линии y=χ и y=g(χ) будут иметь пересечение в искомой точке с абсциссой t.

При этом при всех k=1,2,…m числа χ_k [c;d] и верно неравенство (в теории численных методов оно принимается на основе доказательства соответствующей теоремы):

Это неравенство играет две полезных для нас роли. Во-первых, из него видно, что чем меньше число q, тем быстрее сходится последовательность приближений. Во-вторых, оно является основой для оценки погрешности итерационного метода, о чём чуть ниже.

Если известно, что значения g(χ) находятся в интервале [a;b], то выполнение условия (5) достаточно потребовать лишь на [a;b], тогда необходимость отрезка [c;d] отпадает.

Запишем ещё одно утверждение: если на отрезке [a;b] длины h функция g дифференцируема и

то определяемая формулой (4) итерационная последовательность не сходится к корню t [a;b] ни при каком χ₀ ≠t из этого отрезка.

Геометрически условие (7) означает, что если угол наклона касательной к кривой g в каждой точке из интервала [a;b] получается больше 45⁰, то есть больше угла наклона прямой y=χ (а он равен ровно 45⁰), то тогда линии y=χ и y=g(χ) не будут иметь пересечения в рассматриваемом интервале.

Как уже было сказано выше, оценить погрешность итерационного метода можно из неравенства (6), согласно которому справедлива также следующая формула:

Обозначим f(χ)=χ-g(χ). Функция f(χ) дифференцируема, причём

на [c;d]. Учитывая определение функции f и то, что f(t)=0, имеем . Применив сначала теорему Лагранжа с некоторым числом d_n между t и x_i, а затем неравенство (1), получим:

Модуль разности  можно оценить и сверху. Поскольку , на основе теоремы Лагранжа с числом p_n между  и неравенства (10), получим

Чтобы получить требуемую оценку (8), надо выделить  из неравенства (10) и учесть (11).

Таким образом, если задана точность приближённого корня Δ>0, то итерационный процесс необходимо закончить при выполнении условия

и взять .

Условие (12) и есть оценка погрешности численного метода простой итерации.

Число 0<q<1 может быть принято произвольно, удобно взять q=0,5, тогда из (12) следует

Мы получили строгое математическое обоснование для оценки точности (погрешности) численного диаграммного метода в случае изгиба или внецентренного сжатия. В относительных величинах условие (13) запишется:

При этом по смыслу это и есть критерий сходимости по Чебышеву – см. формулу (2).

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

В СП 63.13330 применяется формула правых либо левых прямоугольников (не уточняется), как наиболее простая.

Получить аналитически критерий сходимости можно на основе остаточного члена той или иной квадратурной формулы (15) в виде

где  заданная предельная погрешность [1].

Изменяя шаг дискретизации Δ_h расчётной схемы, можно добиться выполнения условия (16). При этом , когда .

Исследуем вопрос сходимости численного диаграммного метода на примере железобетонного изгибаемого элемента прямоугольного профиля, b×h=200×500 мм, изготовленного из тяжёлого бетона класса по прочности В10-В60 – var. Для бетона используется усовершенствованная криволинейная диаграмма Карпенко Н.И. [2], для арматуры – двухлинейная Прандтля по СП 63.13330. Площадь сечения нижней растянутой арматуры A_s – переменная (var), верхнее армирование – конструктивное (2Ø12 А500С). При этом общий процент армирования сечения μ=(A_s+A'_s)/(bh) меняется от 0,5 до 3 %. Привязка арматуры к граням бетона: a_s=a'_s=30 мм. Для удобства анализа решение уравнения (1) представим в виде известной зависимости , из неё получим  – предельный изгибающий момент, который способно воспринять рассматриваемое железобетонное сечение. Для определения M_ult,k воспользуемся ранее предложенным алгоритмом [2], но при этом будем вычислять его на каждой итерации, количество которых примем m=1…10 – var. Количество разбиений на элементарные площадки по высоте сечения примем n=20 – const. Начальное приближение кривизны и координаты уровня нулевой линии соответственно:  либо , . Результаты вычисления M_ult,k сведём в таблицы 1 и 2.

Таблица 1

Результаты вычисления M_ult,k, кН×м при

Таблица 2

Результаты вычисления M_ult,k, кН×м при

Представим численные данные графически в системе координат «  – k» – для табл. 1, и «  – k» – для табл. 2, где , ,  – предельный момент соответственно на текущей k-й, 1-й и 2-й итерациях.

Для табличных численных данных на 1-й и 10-й (последней) итерациях построим аппроксимирующую функцию двух переменных в виде полинома 2-й степени (двумерную полиномиальную регрессию порядка n_P=2): , где ,  либо  ( ). Для полиномиальной поверхности порядка n_P количество точек аппроксимации должно быть больше или равно  – условие выполняется. Средствами ПК MahtCAD, применяя встроенную функцию regress/interp, построим следующие поверхности 2-го порядка, представленные на рис. 3.

Для дальнейшего изучения сходимости на рассматриваемом примере уравнение (5) запишем так:

где q, напомним – произвольное число, удовлетворяющее неравенству 0<q<1. Примем q=0,5.

В выражении (17) производную слева представим численно через соответствующие приращения функции  и аргумента k+1-k=1, тогда:

Результаты расчёта левой части формулы (16) представим в таблицах ниже.

Таблица 3

Результаты вычисления левой части (16) при

Таблица 4

Результаты вычисления левой части (16) при

а)                                                                                  б)

Рис. 2. Графики зависимости «  – k» при  (а) и «  – k» при  
(б) для различных вариантов конструирования железобетонного сечения 1…9 (см. табл. 1, 2)

а)                                                                                                б)

в)                                                                                                г)

Рис. 2. Графики двумерной полиномиальной регрессии  для случая  (а, б) и  (в, г):
1 – на 10-й (последней итерации), 2 – на 1-й итерации (а, б) либо 2-й итерации (в, г)

Таблица 4

Результаты вычисления левой части (16) при

В таблицах 3 и 4 цветом выделены значения производной, удовлетворяющие неравенству (16).

Установлены следующие закономерности:

1 – для всех рассмотренных вариантов конструирования итерационный процесс вычислений сходится после 6-й итерации при  и после 4-й итерации при , при этом относительная погрешность расчёта составляет менее δ<1 %;

2 – с увеличением класса бетона по прочности отношения  и  падают, исключение составляют сильно армированные элементы при , для них  – растёт; соответственно в первом случае увеличение прочности бетона, а во втором – её уменьшение, приводит к необходимости увеличения числа итераций для обеспечения сходимости с заданной точностью;

3 – с увеличением процента армирования сходимость расчёта улучшается; например, при количестве итераций равном 4, погрешность при варианте конструирования В60, μ=0,5 % составляет 10,3 %, а при В35, μ=3,0 % – 0,98 %;

4 – для слабо- и среднеармированных элементов приближение к решению задачи (к корню уравнения ) происходит слева, следовательно, тангенс угла карательной здесь отрицательный, а функция  является убывающей, для сильно армированных элементов наоборот: приближение к решению задачи происходит справа, следовательно, тангенс угла карательной здесь положительный, а функция  является возрастающей; и при каких-то комбинациях значений параметров В и μ может быть достигнут минимум функции , в этом случае уже на 2-й итерации будет обеспечена нулевая погрешность расчёта;

5 – для улучшения сходимости в инженерных расчётах рекомендовано принимать начальное значение кривизны , а количество итераций – не менее 6, что обеспечит не превышение относительной погрешности расчёта δ<1 %.

Теперь перейдём к вопросу о влиянии числа разбиений n расчётного сечения железобетонного элемента на сходимость дискретизации нелинейной деформационной модели. На том же примере, рассмотренном выше, для характерных вариантов конструирования № 3 и 7 вычислим на 10-й итерации значение предельного момента  при числе разбиений равном m=2,5, 10, 15 и 20. Полученные результаты сведём в нижеследующие таблицу и график.

Таблица 5

Результаты вычисления M_ult,10, кН×м в зависимости от числа разбиений сечения n