Численно-аналитический алгоритм метода начальных параметров для расчета балок на упругом основании Текст научной статьи по специальности «Математика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2009. Вып. 2

УДК 624.04:519.6 А. В. Матросов

ЧИСЛЕННО-АНАЛИТИЧЕСКИЙ АЛГОРИТМ МЕТОДА НАЧАЛЬНЫХ ПАРАМЕТРОВ ДЛЯ РАСЧЕТА БАЛОК НА УПРУГОМ ОСНОВАНИИ

Введение. Метод начальных параметров используется в расчетной практике для решения одномерных краевых задач строительной механики в виде различных модификаций [1, 2]. Его идея заключается в замене исходной краевой задачи для системы обыкновенных дифференциальных уравнений, описывающей поведение конструкции, задачей Коши, в которой начальные условия определяются из заданных краевых условий. Основной процедурой метода, называемой процедурой переноса граничных условий (ГУ), является выражение компонентов вектора решения на одном из концов конструкции через их значения, определенные на другом конце конструкции. Наличие указанной зависимости позволяет построить систему алгебраических уравнений для нахождения таких начальных условий, при которых решение системы дифференциальных уравнений будет удовлетворять заданным краевым условиям задачи.

Алгоритм реализации процедуры переноса ГУ строится на основе разбиения промежутка интегрирования системы дифференциальных уравнений на интервалы, в которых переменные параметры системы могут быть аппроксимированы постоянными величинами, а само решение может быть с достаточной степенью точности аппроксимировано несколькими членами разложения в ряд Тейлора в окрестности левого конца интервала. Значение вектора решения на правом конце интервала вычисляется через начальные условия, заданные на его левом конце, в качестве которых берутся компоненты вектора решения, вычисленные на правом конце предыдущего интервала разбиения. Начав с первого промежутка, можно получить значения неизвестных величин на правом конце последнего промежутка через начальные условия, определенные на левом конце конструкции. Следует заметить, что метод начальных параметров обычно относят к классу численных методов решения краевой задачи, хотя в его основе лежит аналитическое представление решения в виде степенного ряда.

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

Матросов Александр Васильевич — доцент кафедры информационных систем факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 47. Научное направление: численно-аналитические алгоритмы механики деформируемого твердого тела. E-mail: avmatrosov@mail.ru.

© А. В. Матросов, 2009

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

Аналитические возможности этих систем позволяют вернуться к истокам некоторых «численных» методов и реализовать для них аналитические или численноаналитические алгоритмы, а реализованные в них алгоритмы рациональной и вещественной арифметик - преодолевать вычислительную неустойчивость алгоритмов.

Рис. 1. Расчетная схема балки на упругом основании Объяснение в тексте.

В настоящей работе представлен численно-аналитический алгоритм для метода начальных параметров и приведены результаты вычислительных экспериментов по исследованию сходимости используемых в нем рядов и анализу его вычислительной неустойчивости с помощью системы Maple.

Математическая модель. Рассмотрим тонкую балку длиной l с относительно небольшой шириной b и с переменной вдоль ее оси жесткостью на изгиб, лежащую на упругом основании Винклера с переменным коэффициентом постели c(x) и загруженную распределенной по всей длине нормальной к оси балки непрерывной нагрузкой q(x). Предполагается свободное без трения проскальзывание балки по плоскости контакта с основанием. Введем прямоугольную систему координат, совместив ее начало с левым концом балки (рис. 1). Дифференциальное уравнение четвертого порядка изгиба балки

[EI(x) w"(x)]'' + k(x) w(x) = q(x),

в котором w(x) - прогиб балки по оси z, E - модуль упругости материала балки, I(x) -момент инерции поперечного сечения балки относительно оси, проходящей через центр тяжести сечения, k(x) = c(x)b, запишем в виде системы обыкновенных дифференциальных уравнений в матричной форме

W'(x) = A(x) W(x) + B(x), (1)

где

W(x

w(x)

0(x)

M(x)

N(x)

; A(x

EI (x) ООО

-k(x

О

О

О

О

; B(x) = О

q(x)

9(x) - угол поворота сечения x; M(x) и N(x) - соответственно момент и перерезывающая сила в сечении x.

Завершает постановку задачи задание ГУ на концах балки, по два на каждом ее конце. Например, в случае свободных концов следует положить равными нулю момент и перерезывающую силу в сечениях x = О и x = l, а если какой-либо конец жестко защемлен, то нулю равны прогиб и угол поворота.

Модельная задача. Во всех экспериментах рассматривалась бетонная (E = 2О ГПа) балка длиной l = 2Б м, шириной b = О.1 м, высотой h = О.1Б м, лежащая на упругом основании Винклера с постоянным коэффициентом постели c(x) = 1ОО МН/м3. На левом конце балки действует нормальная к ее продольной оси сила N (О) = 1 кН, угол поворота сечения 0(О) = О. Правый конец балки свободен от нагрузок (M(l) = О, N(l) = О). Для исключения отрыва точек подошвы балки от основания к ней дополнительно была приложена равномерно распределенная по всей длине нагрузка интенсивности q(x) = 1 кН/м. По существу рассматривается правая половина балки длиной 50 м, нагруженной в центре сосредоточенной и по всей длине равномерно распределенной силой.

Все представленные в работе алгоритмы реализованы в системе аналитических вычислений Maple [4], позволяющей моделировать выполнение численного процесса с практически произвольным количеством цифр в мантиссе представления вещественных чисел.

Численный подход [5]. Традиционный алгоритм метода начальных параметров основан на разбиении балки по длине на n участков точками О = xq < xi < ... < xn-l < xn = l так, чтобы все Дk = xk — x,—i, k = 1,...,n, были достаточно малы по сравнению с длиной балки l и на этих интервалах можно было бы пренебречь изменениями матриц A( x) и B (x) из уравнения (1) и положить

A(x) = A(xk) = Ak, B(x) = B(xk) = Bk,

когда xk-l ^ x < xk.

Тогда на интервале x,—i ^ x < x, поведение механической системы может быть описано системой дифференциальных уравнений с постоянными коэффициентами

W'(x)= Ak W(x)+Bk.

(2)

Решение этой системы ищется в виде ряда Тейлора в точке х = хк-1. Тогда значение Wk вектора W(x) в точке х = хк может быть найдено по следующей формуле:

Wk

1

l=0

(l) д1 k-lAk,

(3)

где wkг)l - вектор 1-х производных компонентов вектора решения W(x) в точке х = хк-1, которые вычисляются последовательным дифференцированием системы (2):

W

(l)

kl

A—1 (Ak Wk-і + Bk), l =1,...ж, W

(0)

kl

W

k l.

Подставив рассчитанные значения производных в точке х = хк-1 в (3), получим выражение вектора решений W(x) в точке х = хк через его значение в точке х = хк-1

Wk = Л*к + В£.

(4)

Здесь введены следующие обозначения:

(5)

Используя представление (4) вектора решения W(х) на правом конце к-го интервала через его значение на левом конце, можно выразить вектор решения в точке х = хк через его значение Wo в точке х = хо

При к = п уравнение (6) называется уравнением переноса граничных условий с начала на конец интервала изменения независимой переменной х. Оно также позволяет определить вектор решения W (х) в любой точке хк как функцию, зависящую от параметров внешней нагрузки и значения вектора решения Wo = W(x)|ж=o на левом конце балки, который называется вектором начальных параметров.

Для определения четырех компонентов вектора Wo следует использовать заданные граничные условия задачи. Например, если оба конца балки свободны от нагрузок, то Woз = 0 и 04 = 0, а Wol и Wo2 получим из алгебраической системы уравнений

После нахождения компонентов вектора Wo определение компонентов вектора W(х) в точках разбиения балки по ее длине не представляет сложности - достаточно последовательно для к = 1,...,п применять соотношение (6).

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

Данный алгоритм обладает вычислительной неустойчивостью, связанной с накоплением ошибок округления, что проявляется в невозможности получения достоверных результатов на правом конце длинных балок. Для одного класса задач простое увеличение количества интервалов разбиения позволяет практически обойти указанную проблему, для другого требуется дополнительно применять процедуру ортогонализа-ции вектора начальных параметров на каждом интервале [3], а для третьего ни одно из предложенных решений не дает возможности найти достоверное решение в окрестности правого конца балки.

(6)

А+3,^0! + Л+3^02 + В + 3 = 0, А+4,^0! + Л+4^02 + В + 4 = 0.

Рис. 2. Результаты вычислительных экспериментов расчета балки по численному алгоритму метода начальных параметров

Объяснение в тексте.

На рис. 2,а показаны вычисленные прогибы балки при разбиении на 50 (кривая 2) и 500 (кривая 3) интервалов. Удерживалось 5 членов ряда в решении и длина мантиссы была 16 цифр. Точное решение представлено кривой 1. Этот вычислительный эксперимент демонстрирует именно тот случай, когда увеличение количества интервалов разбиения не устраняет вычислительной неустойчивости алгоритма.

Проблема решается использованием большего количества цифр в мантиссе представления вещественных чисел. На рис. 2,б приведены результаты расчета прогиба балки при 40 интервалах разбиения с мантиссой длиной 17 (кривая 2) и 18 цифр (кривая 3). Кривой 1 здесь также представлено точное решение. Видно, что с ростом количества цифр в мантиссе решение на правом конце балки все ближе подходит к точному решению. Полное совпадение получено при 53 интервалах разбиения с длиной мантиссы 26 цифр при удержании 5 членов ряда в решении. При 40 интервалах разбиения для получения точного решения с использованием мантиссы длиной 26 цифр пришлось повысить количество удерживаемых членов ряда до 6.

Таким образом, проведенные эксперименты показывают, что вычислительная неустойчивость алгоритма может быть преодолена как выбором длины мантиссы в представлении вещественных чисел, так и разбиением балки на определенное количество интервалов и удержанием соответствующего числа членов ряда (3).

Численно-аналитический подход [6]. На первый взгляд может показаться, что алгоритм метода начальных параметров позволяет получить только дискретное решение в точках, разбивающих балку по ее длине на заданное число интервалов. Для нахождения решения в промежуточных точках следует уменьшать интервалы разбиения и снова строить уравнение переноса граничных условий для большего числа точек разбиения. Однако, если использовать возможности системы аналитических вычислений выполнять операции со степенными рядами, то метод начальных параметров позволяет построить аналитическое (в форме степенного ряда) решение, определенное по всей длине балки, а не на дискретном множестве точек.

Предположим, что все функции в уравнении (1) могут быть представлены сходящимися рядами Маклорена на всей длине балки х € [0, /]:

О ОО

А(х) = ^ Агхг, В(х) = ^ Вгхг. (7)

г=0 г=0

Здесь Аг являются числовыми матрицами размерности (4 х 4), а В* - числовыми векторами размерности (4 х 1), компоненты которых определяются из вида соответствующих компонентов матрицы А(х) и вектора В(х).

Решение дифференциального уравнения (1) также ищем в виде ряда Маклорена

О

^х)=Е ^гх1 (8)

г=0

с неизвестными коэффициентами Wг - числовыми векторами размерности (4 х 1), которые следует определить из условия удовлетворения дифференциальному уравнению задачи и заданным граничным условиям.

Подставим в дифференциальное уравнение (1) разложения в степенные ряды вектора решения (8) и коэффициентов системы (7). Из условия равенства двух степенных рядов

ОО

]Г (з + 1) Wj+1xj = £ | £ Ай W; + В,- I х^

3=0 ,'=0 \к+1=, )

получим систему линейных уравнений относительно неизвестных коэффициентов в разложении решения W(х) приравниванием коэффициентов при одинаковых степенях переменной х

Wl = AoWo + В0,

= А1^^0 + А0^^1 + В1,

3Wз = A2Wo + AlWl + AoW2 + В2,

п—1

^2 Ап-1-1^1 + Вп

і=0

Неизвестные векторы Wi, і = 1,..., то, этой бесконечной системы могут быть представлены через вектор начальных параметров Wo путем последовательных подстановок в правые части уравнений выражений для входящих в них векторов Wi через вектор начальных параметров:

Wl = AoWo + Во,

1^2 = ^(Аі + А^о + ^АоВо + Ві),

Wз = - ^А2 + -А0Аі + 2^-0^ + з ^АіВо + -АдВо + -А0Ві

Теперь решение (8) системы дифференциальных уравнений изгиба балки (1) может быть записано в виде

(ОО \ ОО

YJW*ixi\ Wo + ^В*Xі, (9)

і=0 ) і=0

где матрицы W* размерности (4 x 4) и векторы B* размерности (4 x 1) представляются выражениями, в которые входят известные матрицы Ak и векторы Bk из (Т).

Завершает решение задачи определение четырех компонентов вектора W0 из системы четырех линейных уравнений, получающейся из условия удовлетворения искомого решения граничным условиям (по два на каждом конце балки). Например, если оба конца балки защемлены (перемещение и угол поворота равны нулю), то два неизвестных компонента вектора начальных параметров W0 в этом случае уже определены из начальных условий на левом конце балки Woi = w(G) = G и W02 = 0(G) = G. Оставшиеся неизвестными два компонента Woз и W04 вектора начальных параметров находятся из следующей системы линейных уравнений:

/ ОО \ / О \ О

EW* 1 зіЧ Woз + £ W* 1 4іЧ W04 + Е B* ,1і = g,

i=0 i=0 i=0

/ О \ / О \ О

Е W*2 зіМ Woз + Е W*2 4l0 W04 ^ B*21і = G.

i=0 i=0 i=0

Реализовать предложенный алгоритм средствами традиционного языка программирования высокого уровня (C, Pascal, Fortran, Java), конечно, можно, но для этого потребуется прежде всего разработать инструментальные средства (библиотеку объектов или процедур) работы со степенными рядами. С помощью системы аналитических вычислений, например Maple, это делается быстро - в ней уже реализован аппарат работы не только со степенными рядами, но и с любыми аналитическими выражениями.

Предлагаемый аналитический подход к реализации метода начальных параметров также при определенных параметрах задачи может быть вычислительно неустойчив. Однако используемая аналитическая среда Maple позволяет «настроить» ее так, чтобы указанный недостаток алгоритма не проявлялся.

На достоверность результатов расчета по методу начальных параметров влияют два фактора: сходимость рядов решения (9), используемая длина мантиссы в представлении вещественных чисел. Проведенный ряд вычислительных экспериментов для модельной задачи балки на упругом основании имел целью показать влияние указанных факторов на результаты, получаемые по численно-аналитическому алгоритму. В них исследовалось взаимное влияние длины мантиссы и суммируемости рядов решения на правом конце балки (x = l) на примере ряда

О

Y,W*i,ili (10)

i=0

в представлении перемещения

4 О О

Wi(x) = Е Е W*1 ,kx4 Wok + E B* ixi. k=1 \i=0 / i=0

На рис. 3,I,a представлена кривая прогиба w(x) балки при расчетах с 22 цифрами в мантиссе. Ряд (10) для этого случая вычислен с максимально возможной точностью. На рис. 3,1,6 показан график его частичных сумм в полулогарифмической шкале и отображен номер члена ряда, добавление которого к уже полученной сумме не изменяет ее значения. Хотя сумма ряда и определена с предельно возможной для данного случая

Рис. 3. Прогиб ,ш{х) балки (а) и частичные суммы (б) ряда (10) при проявлении вычислительной неустойчивости (I) и без нее (II)

Объяснение в тексте.

точностью, однако она далека от истинного значения, а это и приводит к недостоверным вычислениям на правом конце балки, что и наблюдается на рис. 3,1,а. Найденная сумма ряда 5 = 6806030019080498611941 • 1012, тогда как точное значение равняется 6806030019080479887086655469523254,9806. Погрешность составляет величину порядка 1020.

Рис. 4. Прогиб *ш(х) балки (а) и частичные суммы (б) ряда (10) при коэффициенте постели с(х) = 140 МН/м3 при проявлении вычислительной неустойчивости (I) и без нее (II)

Объяснение в тексте.

Увеличение длины мантиссы до 26 цифр приводит к достоверным результатам, что можно видеть из графиков рис. 3,11. В этом случае для вычисления суммы ряда необходимо учесть 164 члена ряда и ее значение равняется 5 = 68060300190804798871363998 х

X 108.

Однако с ростом коэффициента постели упругого основания до 140 МН/м3 расчеты с 26 цифрами в мантиссе снова приводят к проявлению вычислительной неустойчивости алгоритма (рис. 4,1), повторяя ситуацию вычислительного эксперимента при с(х) = 100 МН/м3 и расчетах с мантиссой длиной 22 цифры (см. рис. 3,1).

Увеличение длины мантиссы до 28 цифр и количества удерживаемых членов рядов решения до 176, приводит к достоверным результатам (рис. 4,II).

Заключение. В работе показаны преимущества использования аналитических систем для реализации вычислительно-неустойчивых алгоритмов. Возможность выполнения расчетов с мантиссой практически неограниченной длины (в Maple максимальную длину мантиссы можно устанавливать равной 268435448) по существу снимает проблему вычислительной неустойчивости численных и численно-аналитических алгоритмов. Однако за такое удобство приходится платить скоростью выполнения программ (системы аналитических вычислений интерпретируемы), а для ее увеличения наличием на компьютере оперативной памяти большого объема (для некоторых задач до 4 Гб).

Литература

1. Крылов А. Н. О расчете балок, лежащих на упругом основании. Л.: Изд-во АН СССР, 1931. 154 с.

2. Власов В. З., Леонтьев Н. Н. Балки, плиты и оболочки на упругом основании. М.: Изд-во физ.-мат. лит-ры, 1960. 491 с.

3. Годунов С. Г. О численном решении краевых задач для системы обыкновенных дифференциальных уравнений // Успехи мат. наук. Т. XVI, вып. 3. М.: Изд-во физ.-мат. лит-ры, 1961. С. 171—174.

4. Матросов А. В. Maple 6. Решение задач высшей математики и механики. СПб.: БХВ-Петербург, 2001. 528 с.

5. Постнов В. А. Численные методы расчета судовых конструкций. Л.: Судостроение, 1977. 279 с.

6. Матросов А. В. Метод начальных параметров: аналитический подход // Информационные технологии в образовании и науке: Материалы Междунар. науч.-практ. конференции ИТО «Поволжье-2007» / под ред. проф. Ю. Г. Игнатьева. Казань: Изд-во «Фолиант», 2007. C. 390—395.

Статья рекомендована к печати проф. Ю. М. Далем.

Статья принята к печати 25 декабря 2008 г.