CC BY
44
УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м V 1974
№ 1
УДК 629.735.33.015
ДИНАМИЧЕСКИЙ РАСЧЕТ ТОНКОСТЕННЫХ ПОДКРЕПЛЕННЫХ ОБОЛОЧЕК С ПРИМЕНЕНИЕМ МЕТОДА КОНЕЧНОГО ЭЛЕМЕНТА
В. И. Зайцев
Рассмотрена задача получения динамических характеристик тонкостенных подкрепленных оболочек. Матрица жесткости для таких оболочек построена методом конечного элемента. В качестве основного рассмотрен плоский трапециевидный элемент, работающий на сдвиг. Масса всей конструкции дискретизирована и сосредоточена в узлах. Система дифференциальных уравнений приведена к системе алгебраических уравнений, которые решаются итерационным методом. В качестве примера приведены результаты расчета подкрепленной цилиндрической оболочки.
Расчетная схема конструкции. Матрица жесткости конструкции. Рассмотрим тонкостенную конструкцию, образованную двумя искривленными поверхностями, расстояние между которыми от точки к точке меняется. Оболочка подкреплена в двух направлениях, которые могут быть и не ортогональными.
Предположим, что напряженное и деформированное состояние всех частей оболочки подчиняется закону Гука. Напряженное состояние обшивки является безмоментным. Работу обшивки и стенок на нормальные напряжения учитываем путем присоединения некоторой части площади их поперечного сечения к площади поясов. В работе используем те же понятия .узла“, .перемещения в узле”, .отсек“ и .главные направления в узле“, что и в работе [2].
Предположим, что в пределах отсека все геометрические, массовые и жест-костные характеристики являются постоянными и равны средним значениям; •оси искривленных элементов заменяем прямыми линиями, а искривленные поверхности — плоскостями.
Масса элементов конструкции сосредоточена в узлах. Колебания оболочки предполагаются малыми около положения устойчивого равновесия. Гистерезис конструкции и материала отсутствует.
Матрицу жесткости всей конструкции получаем путем суммирования жесткостей отдельных элементов по главным направлениям в узлах [2]. Для этого матрицы жесткости отдельных элементов преобразуем с помощью конгруэнтного преобразования в матрицы жесткости по главным направлениям в узлах.
Если сгруппировать векторы усилий и перемещений по признаку их принадлежности к одной из линий сетки (условно назвав их группой), а в пределах группы расположить их по принадлежности к одному узлу и в порядке возрастания номера узла, то матрица жесткости конструкции будет иметь квазидевя-тиленточный вид. Каждая подматрица при этом будет квадратной и иметь порядок, равный числу неизвестных в узле.
Если оболочка замкнутая, то в верхнем правом и нижнем левом углу будут располагаться две дополнительные диагонали, которые учитывают взаимное влияние первой и последней групп.
Матрица масс системы. Масса элементов системы сосредотачивается в узлах. Основой для получения матрицы масс элементов служит кинетическая
энергия элемента, дифференцируя которую дважды по скорости перемещений, получим матрицу элемента
Сначала кинетическую энергию запишем в местных осях координат, а потом преобразуем к главным направлениям в узлах. Если считать, что движения масс, сосредоточенных в узлах, не зависят друг от друга, то, следуя вышеприведенной структуре векторов, получим диагональную матрицу масс. Однородные граничные условия могут быть наложены в любом из узлов по соответствующим направлениям. В векторах усилий и перемещений соответствующие компоненты исключаются, а из матриц масс и жесткости члены, соответствующие этим компонентам, удаляются. В данной работе, при использовании итерационного метода для решения системы уравнений, исходный вектор перемещений выбирался с учетом граничных условий. Матрица жесткости преобразовалась следующим образом: в строку, соответствующую закрепленным компонентам, заносились нули, а элемент на главной диагонали заменялся единицей. Диагональная матрица масс при этом может оставаться неизменной. Если использовать недиагональную матрицу масс, то над ней следует проделать ту же процедуру, что и над матрицей жесткости.
Решение системы уравнений. Систему алгебраических уравнений вида
Kr = R (1)
в данной работе решаем итерационным методом с неполной релаксацией [3].
Этот метод дает хорошую сходимость при небольшом числе итераций.
В данной работе используем метод неполной групповой релаксации, т. е. сразу изменяем не одну, а все компоненты перемещений, отнесенные к одному узлу.
Формула для вычисления компонентов результата п+1 цикла будет иметь вид
{/-i}("+1) = {n}(n)+ч (ад—1 (т + 2 [*у] (о- }(п+1) - 2 {^}(П)) >
V 1=1 j=i 1
где {г/} — вектор перемещений ¿-го узла; [А/у] — квадратная матрица жесткости,
порядок которой равен числу неизвестных в узле; {/?;}— вектор внешней на-
грузки, приложенной к /-му узлу; q — коэффициент релаксации.
Вопрос о выборе множителя q, приводящего к процессу с наибольшей быстротой сходимости, в общем случае не решен [3]. Неполная релаксация называется нижней, если 0<<7<1, и верхней, если В данной работе
были взяты четыре значения множителя q, равные 1; 1,2; 1,85; 1,95, причем последнее значение дало наиболее быструю сходимость. Заранее заданная точность вычислений обеспечивала в окончательно полученном векторе перемещений относительную погрешность между компонентами менее 1%.
Для реализации итерационного процесса необходимо систему уравнений динамического равновесия Кг + Mr = RB привести к виду (1).
Колебания оболочки под действием нагрузки, изменяющейся по гармоническому закону. Пусть вектор нагрузки изменяется во времени по гармоническому закону с частотой ш
Лв = Rbc = cos <at.
Эту систему уравнений можно привести к виду (1)
Kr = Rc, (2)
где
Rc — Rbc + “2 Mr.
Собственные колебания оболочки. Уравнения динамического равновесия в этом случае имеют вид
Кг+ Мг = 0. (3)
Решение этой системы будет иметь вид
г = г cos vi. (4)
Систему дифференциальных уравнений (3) с учетом (4) приведем к системе алгебраических уравнений вида
Kr = Rn,
где
Rn = v2 М г.
Результаты расчета. Для примера была рассмотрена круговая цилиндрическая оболочка, подкрепленная стрингерным набором. По торцам оболочка шарнирно закреплена.
Ниже приведены параметры оболочки.
Радиус 1?, м...........................................0,1775
Длина I, и.............................................0,988
Толщина обшивки 8, м................................... 0,464 • 10~3
Модуль Юнга Е, Па......................................20,0 • 10к>
Плотность материала р, кг/м3 .......................... 7840
Площадь стрингера Рс, м2............................0,105* 10“*
Коэффициент Пуассона V....................:............0,3
Для построения матриц жесткости и масс на оболочку были нанесены два семейства линий. Линии а — параллельны образующим и отстоят друг от друга на я/8. Семейство линий {3 параллельно направляющим линиям цилиндра. Число узлов 128, число перемещений в каждом узле равно пяти, число неизвестных равнялось 640. Узлы, расположенные на одной линии а, объединяются в одну группу. Число групп равнялось 16. По универсальной программе, составленной для ЭЦВМ БЭСМ-6, были сформулированы матрицы масс и жесткости системы, а потом получено решение системы уравнений типа (1) итерационным методом с неполной релаксацией.
Определение динамической реакции оболочки на гармоническую силу. Рассмотрим оболочку,подкрепленную в про -дольном направлении четырьмя стрин-
герами. В одной из точек по направлению нормали к оболочке приложена сила Р— 100 кг с частотой 10 Гц. Формы перемещений, полученные с помощью вышеприведенного алгоритма, приведены на фиг. 1 и 2. Кривая 1 соответствует статическому решению, полученному с помощью технической теории [1], если число членов ряда равно 33. Кривая 2 соответствует статическому расчету по приведенному алгоритму, ш = 0, а кривая 3 соответствует динамическому расчету, <о = 10 Гц. Число приближений для статической и динамической задачи 14 и 16 соответственно. Время получения вектора перемещений для статической задачи равно 10 мин, для динамической —11 мин.
Заметно, что максимальное значение прогиба, полученное по технической теории и с помощью метода конечных элементов, отличается незначительно. В остальных точках это отличие довольно значительное. Это можно объяснить тем, что решение для технической теории в двойных тригонометрических рядах довольно медленно сходится. Если проследить изменение этого решения (кривая 1) в зависимости от числа членов ряда, то можно заметить, что кривая, соответствующая этому решению, с увеличением числа членов ряда от чистой синусоиды постепенно приближается к кривой 2
В окружном направлении для окружности, где приложена сила (см. фиг. 2), отличие в прогибах, полученных по технической теории и методом конечных элементов, незначительное.
Расчет собственных форм и собственных частот. Для оболочки, имеющей подкрепления в виде 4, 8 и 16 стрингеров, были получены четыре первые собственные формы и соответствующие им частоты. Форма колебаний в продольном направлении соответствует низшей частоте (т. е. число волн в продольном направлении т = 1). Форма в окружном направлении не соответствует низшей частоте. Число волн в окружном направлении' п, соответствующее низшей частоте, равно пяти (фиг. 3). Форма с числом волн в окружном направлении, равным четырем, имеет частоту выше, чем форма, содержащая пять волн. Частоты, полученные для оболочки с 4 и 8 •стрингерами, сравнивались с частотами, полученными экспериментальным путем [4]. Различие между величинами частот, полученными из расчета и экспериментально, довольно незначительное(см. таблицу). •
Следует обратить внимание на то, что формы в окружном направлении не совсем гармоничны, т. е. в местах расположения стрингеров амплитуда меньше, чем в неподкрепленных местах. Это, вероятно, можно объяснить подкрепляющим влиянием стрингера на оболочку.
С возрастанием числа стрингеров эта негармоничность уменьшается.
В результате- сопоставления собственных частот, полученных для оболочки с 4, 8 и 16 стрингерами, замечаем, что частоты собственных колебаний практически не зависят от числа стрингеров
(см. таблицу). Фиг. 3
Число волн Частоты собственных колебаний (в герцах) при числе стрингеров, равном
в окружном направлении п 4 8 16
расчет экспери- мент расчет экспери- мент расчет
■ - 4 103,6 103 101,1 100 99,2
5 95.7 94 94,3 93 93,4
6 108,2 106 106,4 105 103,2
7 147.6 143 143.3 138 137,6
Частота собственных колебаний с увеличением числа стрингеров даже «есколько падает. Это можно, наверное, объяснить тем, что изменение кинетической энергии системы стрингеров при установке больше, чем изменение потенциальной. Для получения каждого тона требуется от 14 до 16 приближений. Время расчета всей задачи 30 мин.
ЛИТЕРАТУРА
1. Авдонин А. С. Прикладные методы расчета оболочек и тонкостенных конструкций. М., Машиностроение, 1969.
2. Иванов Ю. И. Расчет подкрепленных тонкостенных конструкций методом конечного элемента. „Ученые записки ЦАГИ“, т. 111,
№ 1, 1972.
3. Фадеев Д. К., Фадеева В. Н. Вычислительные методы
линейной алгебры. М., .Наука“, 1966. _
’ 4. Игл, Сьюолл. Исследование свободных колебаний цилиндрических оболочек с ортогонально расположенными подкрепляющими элементами, рассматриваемыми как дискретные. Ракетная техника и космонавтика, 1968, № 3. .
Рукопись поступила 22/V 1973 г.
