Известия РАН. Механика твердого тела, 2023, № 1, стр. 115-128

1. Введение. Среди аналитических методов, применяемых для решения задач теории упругости, можно выделить такие как метод угловых суперпозиций [1, 2], метод расширения границ [3], метод возмущений [4–6], лучевой метод [7–9], разложения в ряды [10–14] и по функциям Фадля–Папковича [15, 16], двухшаговый метод последовательного возмущения параметров [17], метод быстрых разложений [18–20]. Последний, из перечисленных методов, имеет следующие положительные качества, которыми в совокупности не обладает ни один из известных методов:

1. Доказана сходимость и получена оценка погрешности метода [21].

2. Показана быстрая сходимость используемых рядов Фурье, что позволяет ограничиваться в расчетах небольшим количеством членов ряда и, как следствие, проводить расчет на ЭВМ с высокой точностью при большой экономии времени [22].

3. Метод быстрых разложений применим для решения как линейных [18–20], так и нелинейных задач [23–25], а также задач для криволинейных областей [24, 25], с подвижными границами [25] и фазовыми превращениями [25].

Классическую тригонометрическую интерполяцию чаще всего используют для улучшения качества обработки изображений [26] и восстановления периодических дискретных сигналов конечной длительности [27]. Для решения инженерных задач интегро-дифференциального типа применение классической тригонометрической интерполяции на конечном отрезке проблематично из-за невозможности ее дифференцирования в общем случае и большой ошибки между интерполяционными точками. Устранение подобных недостатков можно осуществить при использовании быстрой тригонометрической интерполяции. Быстрая тригонометрическая интерполяция использовалась при исследовании напряжений в клине [20], расчете траектории космических кораблей [23], решении нелинейного уравнения теплопроводности для криволинейной области с условиями Дирихле [24] и решении двухфазной задачи Стефана с внутренним источником [25]. В данной работе будет исследована погрешность быстрой тригонометрической интерполяции в зависимости от порядка граничной функции и количества членов в рядах Фурье на примере решения задачи о напряжениях в брусе.

2. Постановка задачи. В условиях плоской деформации проекции вектора перемещений материальных точек бруса зависят только от координат x, y:

Компоненты тензора напряжения будут иметь вид:

Запишем уравнения равновесия Ламе для перемещений с учетом массовых сил $X\left( {x,y} \right),$ $Y\left( {x,y} \right)$

К уравнениям (2.3), (2.4) необходимо добавить граничные условия. Будем считать, что упругий брус имеет прямоугольное сечение $\Omega = (0 \leqslant x \leqslant a,\;0 \leqslant y \leqslant b)$. На сторонах бруса зададим условия Дирихле в общем виде

Функции, входящие в граничные условия (2.5) и (2.6), следует подбирать с учетом условий их согласования

выполнение которых позволит найти непрерывное решение задачи (2.3)–(2.6).

В качестве примера функции из (2.5), (2.6) зададим следующим образом

Массовые силы в (2.3), (2.4) запишем выражениями

Зависимости (8) и (9) подобраны так, что задача (3)–(6) имеет точное решение

где K – константа, регулирующая величину перемещений.

Точное решение (2.10) позволит провести исследование погрешности решения краевой задачи (2.3)–(2.6) путем сравнения с приближенным аналитическим решением, полученным методом быстрых разложений. При сравнении будут вычислены: относительная погрешность компонент перемещений (2.1), компонент тензора напряжений (2.2), невязка уравнений равновесия Ламе (2.3), (2.4) и невязка граничных условий (2.5), (2.6).

3. Метод и построение решения. Для решения используем приближенный метод быстрых разложений [21], в соответствии с которым представим $U = U(x,y)$ и V = V(x, y) в виде суммы граничных функции $M_{{2p}}^{U}\left( {x;y} \right),$ $M_{{2p}}^{V}\left( {x;y} \right)$ и ряда Фурье по синусам

Здесь N₁ – число учитываемых членов в рядах Фурье. Граничные функции $M_{{2p}}^{U}(x;y)$ и $M_{{2p}}^{V}\left( {x;y} \right)$ порядка 2p определяются равенствами

где ${{A}_{i}}\left( x \right)$ и ${{B}_{i}}\left( x \right)$, $i = 1 - 2p + 2$ – коэффициенты граничных функций, ${{P}_{i}}\left( y \right)$, $i = 1 - 2p + 2$ – быстрые полиномы [21].

Для исследования влияния порядка граничной функции на точность решения краевой задачи (2.3)–(2.6) в быстрых разложениях (3.1) будем использовать граничные функции второго, четвертого и шестого порядков, которые получаются из (3.2) при p = 1, p = 2 и p = 3 соответственно, т.е.

Быстрые полиномы ${{P}_{i}}\left( y \right)$ и коэффициенты ${{A}_{i}}\left( x \right),$ ${{B}_{i}}\left( x \right)$ определяются равенствами

Для возможности выполнения выражений (3.4) необходимо потребовать, чтобы U = U(x, y) и V = V(x, y) удовлетворяли условию гладкости $\left( {U,V} \right) \in {{C}^{{\left( 6 \right)}}}\left( \Omega \right)$.

Неизвестными в (3.1) являются функции, зависящие только от одной переменной х

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

В (3.6) обозначено через ${{N}_{2}}$ – число учитываемых членов в рядах Фурье. Граничные функции $M_{{2p}}^{{A\left( i \right)}}\left( x \right),$ $M_{{2p}}^{{B\left( i \right)}}\left( x \right),$ $M_{{2p}}^{{u\left( m \right)}}\left( x \right),$ $M_{{2p}}^{{{v}\left( m \right)}}\left( x \right)$ определяем равенствами

где $a_{k}^{{\left( i \right)}}$, $b_{k}^{{\left( i \right)}}$, $u_{k}^{{\left( m \right)}}$ и ${v}_{k}^{{\left( m \right)}}$, $i = 1...2p + 2,$ $p = 1...3,$ $m = 1 - {{N}_{1}}$ – коэффициенты граничных функций вторичных разложений; ${{P}_{k}}\left( x \right)$, $k = 1 - 2p + 2$ – быстрые полиномы [21].

При p = 3 выражения для коэффициентов граничных функций $a_{k}^{{\left( i \right)}}$, $b_{k}^{{\left( i \right)}}$, $u_{k}^{{\left( m \right)}}$, ${v}_{k}^{{\left( m \right)}}$ и быстрых полиномов ${{P}_{k}}\left( x \right)$ имеют вид

Если p = 2, то для описания граничных функций (3.7) из (3.8) следует взять первые шесть полиномов, а в (3.9)–(3.12) учесть первые шесть коэффициентов. При p = 1 граничные функции (3.7) будут содержать первые четыре полинома из (3.8) и первые четыре коэффициента из (3.9)–(3.12).

Таким образом, краевая задача (2.3)–(2.6) сведена к определению 2(2p + 2 + + N₁)$(2p + 2 + {{N}_{2}})$ неизвестных коэффициентов

Значения восьми коэффициентов

входящих в (3.13), находятся при помощи значений компонент перемещений U = = $U(x,y)$ и $V = V(x,y)$ в угловых точках прямоугольной области (см. формулы (3.4), (3.9), (3.10)). С учетом условия согласований (2.7), коэффициенты (3.14) определяются равенствами

Для нахождения остальных $2\left( {2p + 2 + {{N}_{1}}} \right)\left( {2p + 2 + {{N}_{2}}} \right) - 8$ коэффициентов из (3.13) используем быструю тригонометрическую интерполяцию, апробированную в работах [20, 22–25]. Для этого подставим $U = U(x,y)$ и $V = V(x,y)$ из (3.1) в дифференциальные уравнения (2.3), (2.4) и граничные условия (2.5), (2.6). Полученные таким образом выражения в статье не приводим из-за их громоздкости.

Из граничных условий (2.5), (2.6) линейные алгебраические уравнения получим следующим образом. Промежуток [0, b] равномерно разобьем точками y = y_s = = ${{sb} \mathord{\left/ {\vphantom {{sb} {\left( {{{N}_{1}} + 2p + 1} \right)}}} \right. \kern-0em} {\left( {{{N}_{1}} + 2p + 1} \right)}},$ $s = 0,1,..,{{N}_{1}} + 2p + 1$ на ${{N}_{1}} + 2p + 1$ отрезков и запишем уравнения, полученные из граничных условий (2.5) при подстановке $U = U(x,y)$ и $V = V(x,y)$ из (3.1), в каждой внутренней расчетной точке $y = {{y}_{s}},$ $s = 1,..,{{N}_{1}} + 2p$. Будем иметь $4\left( {{{N}_{1}} + 2p} \right)$ линейных алгебраических уравнений. Аналогично, промежуток [0, a] равномерно разобьем точками $x = {{x}_{s}} = {{sa} \mathord{\left/ {\vphantom {{sa} {\left( {{{N}_{2}} + 2p + 1} \right)}}} \right. \kern-0em} {\left( {{{N}_{2}} + 2p + 1} \right)}},$ $s = 0,1,...,{{N}_{2}} + 2p + 1$ на ${{N}_{2}} + 2p + 1$ отрезков и запишем уравнения, полученные из граничных условий (2.6) при подстановке $U = U(x,y)$ и $V = V(x,y)$ из (3.1), в каждой внутренней расчетной точке $x = {{x}_{s}},$ $s = 1,...,{{N}_{2}} + 2p$. Тем самым, будем иметь еще $4\left( {{{N}_{2}} + 2p} \right)$ линейных алгебраических уравнений.

Из дифференциальных уравнений (2.3), (2.4) линейные алгебраические уравнения запишем следующим образом. На область прямоугольника $x \in \left[ {0;a} \right],$ $y \in \left[ {0;b} \right]$ равномерно нанесем сетку в ${{N}_{2}} + 2p + 2$ точках $x = {{x}_{s}} = {{sa} \mathord{\left/ {\vphantom {{sa} {\left( {{{N}_{2}} + 2p + 1} \right)}}} \right. \kern-0em} {\left( {{{N}_{2}} + 2p + 1} \right)}},$ $s = 0,1$, ..., N₂ + + 2p + 1 и в ${{N}_{1}} + 2p + 2$ точках $y = {{y}_{s}} = {{sb} \mathord{\left/ {\vphantom {{sb} {\left( {{{N}_{1}} + 2p + 1} \right)}}} \right. \kern-0em} {\left( {{{N}_{1}} + 2p + 1} \right)}},$ $s = 0,1,...,{{N}_{1}} + 2p + 1$. Для составления системы линейных алгебраических уравнений используются только внутренние точки, образующие сетку из $\left( {{{N}_{1}} + 2p} \right)\left( {{{N}_{2}} + 2p} \right)$ внутренних точек $\left( {{{x}_{s}},{{y}_{s}}} \right)$. Затем, уравнения (2.3), (2.4) при подстановке в них $U = U(x,y)$ и $V = V(x,y)$ из (3.1) запишем в каждой расчетной точке $\left( {{{x}_{s}},{{y}_{s}}} \right)$. В итоге, получаем $2\left( {{{N}_{1}} + 2p} \right)\left( {{{N}_{2}} + 2p} \right)$ линейных алгебраических уравнений. В результате приходим к замкнутой системе $2\left( {{{N}_{1}} + 2p} \right)\left( {{{N}_{2}} + 2p} \right) + 4\left( {{{N}_{1}} + 2p} \right) + 4\left( {{{N}_{2}} + 2p} \right)$ линейных алгебраических уравнений относительно оставшихся $\left( {2\left( {2p + 2 + {{N}_{1}}} \right)\left( {2p + 2 + {{N}_{2}}} \right) - 8} \right)$ неизвестных из (3.13). Данная система уравнений решена в среде Maple. После чего, найденные неизвестные (3.13) подставлены в быстрые разложения (3.1). Тем самым, построено приближенное аналитическое решение краевой задачи (2.3)–(2.6).

4. Полученные результаты и их анализ. В вычислительных экспериментах будем использовать граничные функции (3.3) второго, четвертого и шестого порядков. Количество членов в рядах Фурье первого (3.1) и второго (3.6) быстрых разложений примем одинаковыми, т.е. ${{N}_{1}} = {{N}_{2}} = N$, и выполним расчеты при N = 2…6. В качестве материала бруса выберем тяжелый бетон B30 с характеристиками [28] $E = 32.5$ × 10⁹ Па, ${v} = 0.2$. Тогда коэффициенты Ламе будут равны $\lambda = 9.03$ × 10⁹ Па, $\mu = 1.35$ × 10¹⁰ Па. Величину параметра K и размеры сечения примем равными $K = {{10}^{{ - 6}}},$ $a = 1$ м, $b = 1$ м.

Приближенное аналитическое решение (3.1) сравнивается с точным (2.10). Относительная погрешность перемещений (2.1), тензора напряжений (2.2), невязка уравнений равновесия Ламе (2.3), (2.4) и граничных условий (2.5), (2.6) вычислялась по формуле

где $\Delta $ – абсолютная погрешность,  f_max– максимальное значение исследуемого объекта.

Покажем на рис. 1–4 относительную погрешность расчетов, выполненных при использовании граничной функции четвертого порядка ${{M}_{4}}$ ($p = 2$) для $N = 4$. Из рисунков видно, что максимальная относительная погрешность ${{\delta }_{{\max }}}$ у компонент тензора напряжений ${{\sigma }_{{xx}}}$, ${{\sigma }_{{zz}}}$, ${{\sigma }_{{xy}}}$ и невязки дифференциального уравнения (2.3) достигается в точке (1; 1), а у компонент перемещений U и V, компоненты тензора напряжений ${{\sigma }_{{yy}}}$ и невязки дифференциального уравнения (2.4) – в ее окрестности. Поэтому в табл. 1–3 приведены значения ${{\delta }_{{\max }}}$ всех исследуемых объектов в этих точках для граничных функций второго ${{M}_{2}}$ (p = 1), четвертого ${{M}_{4}}$ (p = 2) и шестого ${{M}_{6}}$ (p = 3) порядков при $N = 2...6$.

Из табл. 1–3 можно увидеть, что точнее всего определяются компоненты перемещений U и V (искомая функция). По сравнению с компонентами перемещений точность нахождения компонент тензора напряжений ${{\sigma }_{{xx}}}$, ${{\sigma }_{{yy}}}$, ${{\sigma }_{{zz}}}$, ${{\sigma }_{{xy}}}$ (первые производные искомой функции) падает на порядок, а точность вычисления невязки уравнений равновесия Ламе (2.3), (2.4) (вторые производные искомой функции) падает на два порядка. Подобная тенденция наблюдается при выборе любого порядка граничной функции и для любого N, а также согласуется с многочисленными вычислительными экспериментами авторов, например [15].

Табл. 1 показывает, что даже при использовании граничной функции второго порядка ${{M}_{2}}$ и небольшом количестве членов в рядах Фурье $N = 3$ достигается точность вычислений компонент перемещений и компонент тензора напряжений (${{\delta }_{{\max }}} < 1\% $) выше, чем точность входных данных из справочников. С увеличением порядка граничной функции и/или количества членов в рядах Фурье (см. табл. 2 и 3) эта точность быстро возрастает. Так, если в быстрых разложениях (3.1) и (3.6) использовать граничную функцию второго порядка ${{M}_{2}}$, то увеличение N количества членов в рядах Фурье на четыре (с двух до шести) ведет к увеличению точности вычислений компонент перемещений, компонент тензора напряжений и невязки дифференциальных уравнений (2.3), (2.4) примерно на один порядок (см. табл. 1). Такое же увеличение точности вычислений можно достичь, используя граничную функцию четвертого порядка ${{M}_{4}}$ (см. табл. 2) и увеличивая N в рядах Фурье на три (с двух до пяти). Если же применить граничную функцию шестого порядка ${{M}_{6}}$ (см. табл. 3), то повышение точности вычислений на порядок можно достичь увеличением N в рядах Фурье на два (с двух до четырех).

Из табл. 1–3 видно, что граничные условия ${{\left. U \right|}_{{x = 0}}}$, ${{\left. U \right|}_{{y = 0}}}$, ${{\left. V \right|}_{{x = 0}}}$, ${{V}_{{y = 0}}}$ выполняются точно при выборе любого порядка граничной функции и при любом количестве учитываемых членов в рядах Фурье. Точность выполнения граничных условий ${{\left. U \right|}_{{x = a}}}$, ${{\left. U \right|}_{{y = b}}}$, ${{\left. V \right|}_{{x = a}}}$, $V{{{\text{|}}}_{{y = b}}}$, а также дифференциальных уравнений (2.3), (2.4), и точность вычисления компонент перемещений и компонент тензора напряжений зависят от выбора порядка граничной функции и количества членов в рядах Фурье. Если увеличить порядок граничной функции (со второго до четвертого или с четвертого до шестого), а N во всех расчетах взять одинаковым, то относительная погрешность ${{\delta }_{{\max }}}$ у компонент перемещений и компонент тензора напряжений уменьшится на два порядка. Также на два порядка уменьшится невязка граничных условий ${{\left. U \right|}_{{x = a}}}$, ${{\left. U \right|}_{{y = b}}}$. Для граничных условий ${{\left. V \right|}_{{x = a}}}$, $V{{{\text{|}}}_{{y = b}}}$ невязка уменьшается на три порядка.

Анализируя поведение максимальной невязки дифференциальных уравнений (2.3), (2.4) можно сделать следующий вывод: увеличение порядка граничной функции со второго до четвертого ведет к уменьшению невязки на один порядок, а при увеличении порядка граничной функции с четвертого до шестого на возрастание точности еще влияет число N. Так, при $N = 2$ уменьшение невязки будет составлять один порядок, при $N = 6$ уже два порядка.

С помощью табл. 4 и 5 можно ответить на вопрос, что более эффективно для уменьшения погрешности вычислений – увеличение порядка граничной функции или увеличение количества членов в ряде Фурье. В табл. 4 представлено количество линейных алгебраических уравнений в системах, которые необходимо решить для расчета данных, указанных в табл. 1–3. В табл. 5 записано время расчета программы в секундах на персональном компьютере с процессором Intel Core i3-4160 и ОЗУ 8 ГБ. Из табл. 4 видно, что одинаковая трудоемкость вычислений получается при выборе следующих комбинаций параметров p и N:

1) p = 1 и $N = 4$ или p = 2 и $N = 2$–120 линейных алгебраических уравнений;

2) p = 1 и $N = 5$ или p = 2 и $N = 3$–154 линейных алгебраических уравнения;

3) p = 1 и $N = 6$ или p = 2 и $N = 4$ или p = 3 и $N = 2$–192 линейных алгебраических уравнения;

4) p = 2 и $N = 5$ или p = 3 и $N = 3$–234 линейных алгебраических уравнения;

5) p = 2 и $N = 6$ или p = 3 и $N = 4$–280 линейных алгебраических уравнений.

Учитывая данные табл. 1–3 можно сделать вывод, что в каждом из пяти указанных вариантов комбинаций параметров p и N точность расчетов выше (значение ${{\delta }_{{\max }}}$ в среднем падает на порядок) в случаях, соответствующих более высокому значению параметра p, т.е. более высокому порядку граничной функции. Принимая во внимание данные табл. 5, отметим, что повышение порядка граничной функции приводит и к сокращению времени расчета. Для пяти комбинаций параметров p и N, соответствующих одинаковой трудоемкости расчетов, уменьшение времени расчета при большем значении параметра p составляет от 2.5% (третья комбинация) до 15% (первая комбинация).

При изучении свойств поля напряжений в брусе наибольший интерес представляет исследование влияния геометрических размеров его сечения на величину интенсивности напряжений $\tilde {\sigma }$ [29]

и расположение точки с наибольшим ее значением $\tilde {\sigma }$.

Размеры сечения в расчетах подобраны таким образом, чтобы площадь сечения бруса оставалась постоянной. Вычислительные эксперименты показали, что качественный вид профилей интенсивности напряжений $\tilde {\sigma }$ квадратного и прямоугольного сечений бруса будут различны. Виды подобных профилей для случаев $a = b = 1$ и $a = 2$, $b = {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}$ изображены на рис. 5, из которого видно, что точка с максимальной интенсивностью напряжений ${{\tilde {\sigma }}_{{\max }}}$ для квадратного сечения расположена на стороне y = 1, а для прямоугольного – находится в угловой точке (2; 0). Значения ${{\tilde {\sigma }}_{{\max }}}$ для различных значений a и b, и координаты точки, в которой достигается ${{\tilde {\sigma }}_{{\max }}}$ записаны в табл. 6. Анализируя данные этой таблицы, можно сделать вывод, что для прямоугольного сечения наибольшая интенсивность напряжений ${{\tilde {\sigma }}_{{\max }}}$ всегда расположена в крайней точке длинной жестко защемленной стороны при x = a. Также из табл. 6 видно, что напряжения в брусе тем выше, чем больше превосходство длины одной стороны над длиной другой стороны прямоугольного сечения. Наименьшая интенсивность напряжений среди всех сечений наблюдается в брусе с квадратным сечением.

5. Заключение. В статье показана эффективность быстрой тригонометрической интерполяции при решении задач теории упругости на примере задачи о напряжениях в брусе прямоугольного сечения. При использовании в быстрых разложениях граничной функции второго порядка и небольшом количестве членов в рядах Фурье (от двух до шести) максимальная относительная погрешность δ_max компонент перемещений и компонент тензора напряжений составляет менее 1% , что является приемлемой погрешностью для большинства технических расчетов. С увеличением порядка граничной функции и/или количества членов N в рядах Фурье погрешность δ_max быстро уменьшается. Увеличение порядка граничной функции является более эффективным способом уменьшения погрешности δ_max по сравнению с увеличением количества членов в рядах Фурье.

При исследовании интенсивности напряжений $\tilde {\sigma }$ в брусе с различными габаритными размерами прямоугольного сечения, но одинаковой площадью всех сечений выяснилось, что наименьшее значение ${{\tilde {\sigma }}_{{\max }}}$ среди всех сечений наблюдается в брусе с квадратным сечением. В этом случае точка с ${{\tilde {\sigma }}_{{\max }}}$ расположена на стороне $у = b$ не далеко от угловой точки $\left( {a;b} \right)$. Для прямоугольного сечения точка с максимальной интенсивностью напряжений ${{\tilde {\sigma }}_{{\max }}}$ будет расположена на длинной стороне в точке $\left( {a;0} \right)$. При этом, чем больше превосходство длины одной стороны прямоугольного сечения над длиной его другой стороны, тем выше ${{\tilde {\sigma }}_{{\max }}}$.