CC BY
3Труды МАИ. Выпуск № 101_ http: //trudymai .ru/
УДК 517.9
Об обеспечении точности решения задач модального анализа
Коровайцева Е.А.
НИИ механики МГУ им. М.В. Ломоносова, Мичуринский проспект, 1, Москва, 119192, Россия e-mail: [email protected]
Аннотация
На примере формулировки задач модального анализа предлагается введение трех простейших канонических форм, позволяющее структурировать изложение алгоритмов решения задач и обоснованно минимизировать их круг. Для базового алгоритма, подразумевающего применение метода начальных параметров, предлагается использование интегрального подхода к контролю точности решения задачи. Эффективность применения данного подхода проиллюстрирована на примере решения задач анализа собственных изгибных колебаний шарнирно опертой балки и колебаний цилиндрической оболочки. Проанализированы отличия результатов решения данных задач различными методами, как реализованными в авторской программе, так и встроенными в математические пакеты.
Ключевые слова: краевые задачи, модальный анализ, метод начальных параметров, сопряженные уравнения.
Введение
Основным требованием к любым вычислительным задачам является существование их решения. При его удовлетворении на втором плане возникают вопросы, связанные с получением этого решения и, в частности, вопросы точности получаемых, как правило, с применением ЭВМ решений, обусловленности задач и сходимости алгоритмов их решения к точному решению. Задачи механики деформируемого твердого тела (МДТТ) после их должного представления формально являются задачами одного из разделов теории дифференциальных уравнений. При использовании процедур численного интегрирования дифференциальных уравнений важное значение имеет устойчивость интегрирования. С обеспечением устойчивости тесно связаны вопросы рационального выбора шага интегрирования. Простейшими приемами выбора являются либо повторение расчетов с уменьшенным шагом для сравнения результатов, либо применение двухсторонних формул различных алгоритмов, не требующих повторных расчетов [1-3]. Но эти алгоритмы часто приводят к неоправданно сильному затягиванию времени расчетов на ЭВМ.
Характерным для работ по анализу устойчивости методов интегрирования (например, [4-6]) и выбору шага является сосредоточение внимания не на исследовании условий устойчивости решения данным алгоритмом различных уравнений и их систем при конечном шаге интегрирования, а на анализе устойчивости собственно алгоритмов (при шаге, стремящемся к нулю). Сложность вопроса обеспечения устойчивости численного решения систем дифференциальных уравнений вынудила разработать частные приемы решения таких задач методами,
устойчивыми в целом ряде случаев лишь в техническом смысле, либо в виде сегментных методов [7], либо в виде прогонок [8-10].
Так как численное интегрирование представляет собой приближённый метод решения уравнений, то кроме вопросов устойчивости важную роль играют средства определения степени приближения к точному решению оценкой погрешностей численного интегрирования. Однако имеющиеся в литературе оценки, например, [11], [12], не являются достаточно эффективными для практического применения, давая в ряде случаев не только количественно, но и качественно неправильное представление о поведении погрешности.
При реализации метода конечных элементов построения решения задач МДТТ зачастую выбор сетки оказывается направленным только на уменьшение ширины ленты матрицы разрешающей системы линейных алгебраических уравнений, то есть на решение чисто вычислительных проблем линейной алгебры, но не на обеспечение точности решения исходной задачи. Таким образом, вопрос отыскания эффективных способов контроля точности и сокращения времени численного решения сложных систем уравнений, описывающих прикладные задачи МДТТ, остается и до сих пор весьма актуальным.
В данной работе впервые продемонстрировано применение авторского алгоритма контроля точности решения линейных краевых задач [13] на примере задач модального анализа простейших элементов конструкций летательных аппаратов. Предложена систематизация одномерных линейных краевых задач МДТТ, на основании которой могут быть разработаны алгоритмы решения задач модального анализа произвольных составных конструкций летательных аппаратов,
отличающиеся большей универсальностью, чем известные в доступной литературе [14-16]. Представленная структуризация задач позволяет использовать в разрабатываемых алгоритмах предлагаемый подход к контролю точности получаемых результатов для обоснованного построения конечно-элементной сетки.
Систематизация одномерных линейных краевых задач МДТТ Наряду с дифференциальными уравнениями, описывающими модели задач МДТТ и записанными в векторно-матричной форме, выделяются три типа дискретных соотношений, замыкающих задачу описания состояния элементов расчетной модели любой прикладной задачи: граничные условия, условия неразрывности обобщенных векторов, составленных из искомых векторов, и условия совместности состояния внутренних элементов расчетной модели. В соответствии с указанными типами соотношений можно выделить три простейшие канонические формы краевых задач МДТТ, математическая формулировка которых для задач модального анализа имеет следующий вид:
1) двухточечные хорошо обусловленные краевые задачи без дополнительных соотношений (каноническая форма I)
^- = А{х,]л,а)-у (1)
ах
Вх{хнГ^аУун= б 1П 2 *„□ (2)
Здесь и далее у(х,]л,{1н), ун{х,]л,].7Н) - векторы разрешающих переменных задачи, А(х,/7,а), 0с„, Д,,гу) - матрицы коэффициентов, /7Ц'), рп - функция и вектор исходных значений параметров задачи, 1 и 2, Н и К - индексы значений переменных
и аргумента в начальной и конечной точках интервала решения краевой задачи, т -собственная частота колебаний конструкции.
2) многоточечные неразветвленные хорошо обусловленные краевые задачи без
дополнительных соотношений (каноническая форма II)
^ = (3)
ах{
А(х1н,Дн^)-Я„=б 1П 2 х1нП ^ (4)
ЦЫ■ УЛхм) = со) ■ у]+,{х]+Хн) ./ е [1,А'-1] (5)
Здесь аналогично (1)-(2) для простоты записи в искомых векторных функциях опущены фактические их зависимости от полной совокупности функций и векторов исходных значений параметров задачи, а матрицы Б. - только числовые, условия
(4) являются граничными, а (5) - условия сопряжения суперэлементов краевой задачи.
3) многоточечные разветвленные хорошо обусловленные краевые задачи без
дополнительных соотношений (каноническая форма III)
^ = Ч1^] (6)
ах{
Б (х )-у + £> (х )-у =0 (8)
3^3 ,н/ т\ т, к / у т,к V /
(9)
0^к)-у^+От(хт^-ут^ = 0 (10)
1 р
Уи(х )-у +УО(х )-у +5 -С (х П т)-у +8 -С (х П ,а)-у =0 (11)
/ т\ т,н ; .Ути/! п\ п,к; у п,к дн д^ д,н> г*д,н> ' ' д,н дк д^ д,к> г*д,к? ' * д,к \ '
т=1 п=\
Здесь значения индексов ], т, I, р, д определяются топологией конкретной краевой задачи. Условия (7) являются граничными, (8)-(10) - условия неразрывности обобщенных перемещений, (11) - условия равновесия узлов. В соотношениях (11) символы 8 , 8 имеют следующие значения: 8 = 1, если в узле
ы2
д используется инерционное слагаемое ^ д2н , 8дк = 1 если в узле д используется
Ч
2,к иначе 8 = 0, 8= 0.
^2 у
, иначе 52н ~Чк
Предлагаемая систематизация соотношений задач модального анализа позволяет использовать для задач любой топологии единый алгоритм контроля точности решения.
Базовый алгоритм решения задачи модального анализа
Простейший алгоритм решения задач (1)-(2), (3)-(5) или (6)-(11) основан на использовании метода начальных параметров [17], согласно которому на /-том суперэлементе искомый вектор может быть представлен с использованием его значения на левой границе суперэлемента в следующем виде:
г е[1, Щ . (12)
В таком представлении матричная функция М^ (х1) является решением следующей задачи Коши:
^ = м1(0) = к. (13)
ахI
Для обеспечения необходимой точности расчета матрицы М^ (х1) применяется описанный в [13] интегральный подход, основанный на контроле ортогональности
нормированных интегральных матриц исходной и сопряженной систем дифференциальных уравнений М1 (х) и Ц (х):
МТ (х) • Ц (X) - Е
<8, (14)
где Е - единичная матрица, 8 - принятое допустимое значение нормы ортогональности нормированных интегральных матриц.
Данный подход применяется на этапе препроцессирования для получения рекомендаций по сегментации задачи. Таким образом, число условий сопряжения суперэлементов для канонических форм II и III увеличивается, а каноническая форма I переходит в каноническую форму II с автоматически определяемым числом сегментов.
Подстановка соотношений (12) в условия совместности деформирования суперэлементов (5) или (7)-(11) приводит к однородной системе линейных алгебраических уравнений относительно векторов уг н, блочная структура которой
определяется самими этими условиями. Условие нетривиальности решения данной системы приводит к частотному уравнению, из которого определяются собственные частоты колебаний конструкции.
Примеры
Тестирование предложенного алгоритма проведено для задач анализа собственных колебаний шарнирно опертой балки и цилиндрической оболочки.
Как указано в [18], «при расчете высших частот и форм собственных колебаний изгиба методом начальных параметров, независимо от того, в какой форме он реализуется, возникают вычислительные трудности». Расчеты частот и
форм изгибных колебаний шарнирно опертой балки показали, что отсутствие сегментации задачи в соответствии с условием (14) приводит к ошибочным результатам для десятой и более высоких частот.
В таблице 1 представлены результаты расчета безразмерной десятой собственной частоты колебаний шарнирно опертой балки квадратного сечения и удовлетворение граничным условиям на торцах балки. Частота ю масштабируется
величиной , перемещение w и координата х - величиной £, изгибающий
ЕЗ
момент М7 - величиной —2, где ЕЗ7 - жесткость сечения на изгиб, р - плотность
г Б 2
материала, Г - эталонная площадь сечения, £ - длина балки. Значения ^(0), М2(0) и w(l), М2 (1) представляют собой граничные условия на левом и правом торцах балки соответственно. Видно, что расчет частот и форм, основанный на непосредственном использовании соотношений (12)-(13) без предварительной сегментации задачи в соответствии с условием (14) приводит к погрешности расчета частоты А«8% и катастрофическому росту значений компонент вектора форм колебаний, что иллюстрируется невыполнением граничных условий.
Таблица 1
Влияние сегментации задачи на результаты расчета десятой собственной частоты изгибных колебаний шарнирно опертой балки
Теория _Расчет
Без С
сегментации сегментацией
ю 31,4159 34,3508 31,4163
w{0) 0 0 0
41) 0 - 2,64 ■ 1012 1,76 ■Ю-13
М2 (0) 0 -60 - 9,34 ■ 10-4
М2 (1) 0 -3,11 ■ 1015 - 2,47 ■ 10-12
На рис. 1 показано изменение нормы ортогональности нормированных интегральных матриц системы дифференциальных уравнений изгиба балки
. Видно, что уже при значении аргумента х > 0.45 решение системы будет
Мт ■ Ь
некорректно без использования специальных математических приемов.
Рис. 1. Изменение нормы ортогональности нормированных интегральных матриц
системы дифференциальных уравнений изгиба балки
При использовании сегментации задачи в соответствии с условием (14) везде в работе принято значение нормы ортогональности нормированных интегральных
матриц 8 = 1.05. Определяемое при расчете десятой собственной частоты изгибных колебаний балки на этапе препроцессирования число сегментов разбиения интервала интегрирования равно трем. Соответствующий график изменения значения 8 по длине балки представлен на рис. 2.
Рис. 2. Изменение нормы ортогональности нормированных интегральных матриц системы дифференциальных уравнений изгиба балки при сегментации задачи
Представленные в табл. 1 и на рис. 2 результаты получены при интегрировании системы (13) по методу Рунге-Кутта четвертого порядка точности с использованием авторской программы. Также были проведены расчеты при использовании функций системы МАТЬАВ, реализующих метод Рунге-Кутта (ode45) и методы, применяемые при решении жестких задач или при необходимости получить высокоточное решение (метод Адамса-Башфорта - функция ode23, метод
Гира - функция ode15s, метод Розенброка - функция ode23s, неявный метод трапеций с использованием свободной интерполяции - функция ode23t). Отмечено, при интегрировании системы (13) даже специализированными методами расчет без использования сегментации задачи приводит к ошибочным результатам (определение частоты с точностью порядка 10%, невыполнение граничных условий). Расчет же с использованием сегментации при применении для численного интегрирования любых встроенных функций MATLAB неэкономичен, так как затраты машинного времени на сегментацию задачи на порядок превышают затраты при использовании авторского программирования метода Рунге-Кутта. Последнее связано с невозможностью самостоятельного задания шага интегрирования в MATLAB. При использовании, к примеру, пакета Mathcad указанная проблема отсутствует. Кроме того, при использовании встроенных функций требуемое число сегментов выше, а точность расчета частоты и удовлетворения граничных условий ниже, чем в случае авторского программирования даже всего лишь метода Рунге-Кутта четвертого порядка точности. При этом уменьшение абсолютной и относительной погрешности вычислений по функциям MATLAB приводит к менее быстрому увеличению точности расчета, чем сегментация задачи. Соответствующие результаты расчетов для десятой собственной частоты представлены в таблице 2. Здесь N - требуемое число сегментов, tcesM - затраты времени на сегментацию задачи, Sc, % - погрешность расчета частоты. Для встроенных функций MATLAB результаты расчетов приведены для значений абсолютной погрешности вычислений Sabs = IQ"6 и относительной погрешности Sre[ = 10"3.
Таблица 2
Результаты расчета десятой собственной частоты изгибных колебаний шарнирно
опертой балки различными методами
Метод N 1сегм, мин. 8(о,°% "(0); М2 (0) "(1); м2 (1)
Рунге-Кутта 3 0:07 6,12 ■Ю-3 0; - 6,13 ■ 10-2 1,97-10"13; -5,07■Ю-11
ode45 7 1:10 6,12 ■Ю-3 0; - 5,38 ■ 10-2 -7,6 ■Ю"8; -6,95 ■Ю"5
ode23 12 0:56 4,01 ■Ю-3 0; - 3,67 ■ 10-2 -2,47 ■Ю"5; -3,10 ■Ю"2
ode15s 17 1:24 6,48 ■Ю-2 0; - 7,09 ■ 10-2 -1,00 ■Ю"4; -5,63 ■Ю"2
ode23s 4 5:09 3,7 ■Ю-1 0; - 3,69 ■ 10-2 -4,73 ■Ю-3; - 4,91
ode23t 3 2:32 3,7 ■Ю-1 0; -1,21 ■ 10-1 -1,61; -1,6 ■Ю3
Полученные результаты позволяют подтвердить точку зрения, что использование готовых программ расчета, будь то конечноэлементные комплексы или встроенные функции сред программирования, требует высокой аккуратности, глубоких знаний особенностей работы этих программ, и даже при этих услових может быть менее эффективно, чем самостоятельное программирование.
Отметим, что использование для интегрирования задач Коши метода Рунге-Кутта более высокого порядка точности не только не позволяет повысить точность решения задачи, а напротив, снижает ее. Так, при применении метода Рунге-Кутта десятого порядка точности без использования метода сегментации погрешность расчета девятой собственной частоты изгибных колебаний балки составляет 8( = 0.17%, что на три порядка выше, чем при применении метода Рунге-Кутта
четвертого порядка точности. Соответствующие формы колебаний для прогиба и угла поворота сечений балки показаны на рис. 3. Видно, что амплитудные значения указанных функций убывают по длине балки, что свидетельствует об ошибочности расчета.
и', г?
1.5
1
0.5 0
-0.5 -1
-1.5 _2
"О 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 X Рис. 3. Результаты расчета форм изгибных колебаний шарнирно опертой балки по методу Рунге-Кутта 10 порядка точности: 1 - прогиб, 2 - угол поворота сечения
балки.
Вторая иллюстрационная задача связана с расчетом низшей частоты собственных колебаний цилиндрической оболочки. В таблице 3 приведены результаты расчетов по уравнениям общей моментной теории оболочек [17] с учетом инерции по всем трем направлениям. Отметим, что использование условия (14) для сегментации задачи позволяет определить минимально необходимое число точек интегрирования на интервале по выполнению данного условия на первом
шаге интегрирования, что особенно актуально при расчете длинных и тонких оболочек.
В таблице используются следующие обозначения: тип граничных условий (Г) - оба края шарнирно оперты, с-с - оба края жестко защемлены, Ь - длина оболочки, Я - радиус, И - толщина, п - номер окружной гармоники, -минимально допустимое число точек интегрирования по длине оболочки, т -принятое число сегментов разбиения, ю - рассчитанная собственная частота, ю* -собственная частота, полученная аналитическим путем по теории пологих оболочек с учетом инерции по всем трем направлениям [19]. Частота масштабируется
величиной Е ■ —. Отметим, что расчет по теории пологих оболочек дает \Р Ь
завышенные значения частоты по сравнению с общей моментной теорией.
Таблица 3
Результаты расчета первой собственной частоты колебаний цилиндрической
оболочки
№ Г Ь/Я Я/И п Nmin т ю ю*
1 s-s 1 100 7 - 2 0,2392 0,2419
2 ^ 1 100 7 - 2 0,3052 0,3076
3 s-s 2 300 7 200 6 0,138 0,1398
4 ^ 2 300 8 200 6 0,1945 0,1964
5 s-s 8 100 3 400 15 0,226 0,245
6 ^ 8 100 3 400 15 0,326 0,341
7 s-s 8 300 3 700 51 0,142 0,1456
8 8 300 4 700 44 0,1984 0,2044
Представленные в таблице 3 результаты получены при интегрировании уравнения (13) по методу Рунге-Кутта в авторской реализации программы. Интегрирование, основанное на использовании встроенных функций MATLAB
численного интегрирования дифференциальных уравнений, подтвердило выводы, полученные при модальном анализе шарнирно опертой балки: расчет без использования сегментации даже коротких оболочек при применении любых встроенных функций, в том числе реализующих методы решения жестких задач, приводит к неверным результатам. При этом минимально необходимое число сегментов и затраты машинного времени выше, чем в случае авторского программирования метода Рунге-Кутта, а точность удовлетворения граничным условиям на несколько порядков ниже.
Также были выполнены расчеты первой собственной частоты цилиндрической оболочки при использовании метода ортогонализации Годунова [10] и выборе числа точек ортогонализации и интегрирования по соотношению (14). Отметим, что в доступной литературе известна только одна работа [20], в которой даются рекомендации по расчету числа точек ортогонализации и интегрирования. При этом подчеркивается, что предложенные соотношения являются приближенными. В таблице 4 приведены результаты для оболочки с геометрией и условиями закрепления, указанными в п.6 таблицы 3. В этом случае число точек ортогонализации, рассчитанное по рекомендациям работы [20], составляет т> 5, число точек интегрирования на всей длине оболочки Л^>160...240. Вместе с тем в соответствии с соотношением (14) число точек интегрирования необходимо выбирать из диапазона N > 400. В таблице 4 приведено сравнение результатов расчетов частоты и удовлетворения граничных условий на правом торце оболочки (заделка) при выборе минимально допустимого числа точек ортогонализации и
интегрирования, рассчитанного по соотношению (14) и в соответствии с рекомендациями работы [20]. В дополнение к обозначениям, принятым в таблице 3, введены следующие: и - меридиональное перемещение, V - окружное перемещение, w - перемещение по нормали к меридиану оболочки, 3 - угол поворота нормали.
Таблица 4
Сравнение расчетов первой собственной частоты цилиндрической оболочки для
разных алгоритмов сегментации задачи
Тип алгоритма сегментации Ытт т ю и(1); V(1); 41); 3(1)
Условие (14) 400 58 0,326 - 3,43; 6,38 • 10-13; -1,24 -10-14; -1,14 -10-10
Условие работы [20] 240 5 0,236 -4,07 -109; - 2,66 • 1011; 1,08 • 1010; 8,01 • 1013
Как следует из таблицы 4, принятие минимально допустимого разбиения отрезка интегрирования на сегменты по рекомендациям [20] приводит к неверным результатам, в отличие от минимально же допустимого разбиения на основании соотношения (14). Таким образом, соотношения работы [20], определяющие сегментацию задачи, являются недостаточно жесткими, и можно рекомендовать использование условия (14).
Заключение
Предложена структуризация одномерных линейных краевых задач модального анализа, которая может быть использована для одних и тех же исходных скалярных
соотношений двояко: во-первых, для описания алгоритма решения задачи, а во-вторых - и для структурного программирования. Приведен базовый алгоритм решения задачи модального анализа, подразумевающий использование интегрального подхода к контролю процесса решения. Применение данного подхода проиллюстрировано примерами, подчеркивающими преимущества авторского программирования при проведении расчетов, а также использования предложенного метода контроля точности решения для получения достоверных результатов.
Работа выполнена при финансовой поддержке РФФИ (проект № 16-38-60074
моладк).
Библиографический список
1. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: БИНОМ. Лаборатория знаний, 2008. - 636 с.
2. Горбунов А.Д., Шахов Ю.А. О приближенном решении задачи Коши для обыкновенных дифференциальных уравнений с наперед заданным числом верных знаков // Журнал вычислительной математики и математической физики. 1963. Т.3. № 2. С. 239 - 253.
3. Лозинский С.М. Недостаточные и избыточные методы численного интегрирования обыкновенных дифференциальных уравнений // Вестник Ленинградского университета. Сер.: Математика. Механика. Астрономия. 1967. № 7. В. 2. С. 74 - 86.
4. Hildebrand F.B. Introduction to numerical analysis. New York, Me Graw - Hill, 1956, 511 p.
5. Carr J.W. Error bounds for the Runge-Kutta singlestep integration process // Journal of the Association for computing machinery, 1958, IT 1, vol. 5, pp. 44 - 59.
6. Hamming R.W. Stable predictor - corrector methods for ordinary differential equations // Journal of the Association for computing machinery, 1959, vol. 6, no. 1, pp. 37 - 47.
7. Калнинс А. Исследование оболочек вращения при действии симметричной и несимметричной нагрузок // Прикладная механика. Серия Е. 1964. Т. 31. № 3. С. 112 - 122.
8. Абрамов A.A. 0 переносе граничных условий для систем линейных обыкновенных дифференциальных уравнений (вариант метода прогонки) // Журнал вычислительной математики и математической физики. 1961. Т. 1. № 3. С. 542 - 544.
9. Бидерман В.Л. Применение метода прогонки для численного решения задач строительной механики // Инженерный журнал. Механика твердого тела. 1967. № 5. С. 62 - 65.
10. Годунов С.К. О численном решении краевых задач для систем обыкновенных линейных дифференциальных уравнений // Успехи математических наук. 1961. Т. 16. В. 3. С. 171 - 174.
11. Демидович Б.П., Марон И.А., Шувалова Э.З. Численные методы анализа. Приближение функций, дифференциальные и интегральные уравнения. - М.: Наука, 1967. - 367 с.
12. Коллатц Л. Численные методы решения дифференциальных уравнений. - М.: Инлитиздат, 1953. - 460 с.
13. Коровайцев А.В., Коровайцева Е.А., Ломовской В.А. Решение прикладных одномерных краевых задач с автоматической точностью // Вестник МИТХТ. 2012. Т. 7. № 6. С. 41 - 45.
14. Алексушин С.В. О расчете собственных частот стабилизатора летательного аппарата на ранних этапах проектирования // Труды МАИ. 2014. № 73. URL: http://www.trudymai.ru/published.php?ID=48459
15. Сысоев О.Е., Добрышкин А.Ю., Нейн Сит Наинг. Аналитическое и экспериментальное исследование свободных колебаний разомкнутых оболочек из сплава Д19, несущих систему присоединенных масс // Труды МАИ. 2018. № 98. URL: http: //www.trudymai .ru/published.php?ID=90079
16. Гришанина Т.В., Тютюнников Н.П., Шклярчук Ф.Н. Метод отсеков в расчетах колебаний конструкций летательных аппаратов. - М.: Изд-во МАИ, 2010. - 180 с.
17. Бидерман В.Л. Механика тонкостенных конструкций. - М.: Машиностроение, 1977. - 488 с.
18. Бидерман В.Л. Прикладная теория механических колебаний. - М.: Высшая школа, 1972. - 416 с.
19. Прочность, устойчивость, колебания. Справочник в трех томах. // Под ред. И.А. Биргера и Я.Г. Пановко. - М.: Машиностроение, Т. 3. 1968.-567 с.
20. Кармишин А.В., Лясковец В.А., Мяченков В.И., Фролов А.Н. Статика и динамика тонкостенных оболочечных конструкций. - М.: Машиностроение, 1975. -376 с.
