CC BY
50
УДК 539.3, 519.63
Численное моделирование деформирования анизогридных конструкций с применением высокоточных схем без насыщения
12 1 © С.К. Голушко ' , Б.В. Семисалов
1 Конструкторско-технологический институт вычислительной техники СО РАН, Новосибирск, 630090, Россия
2 Институт вычислительных технологий СО РАН, Новосибирск, 630090, Россия
Рассмотрен класс перспективных анизогридных конструкций, представляющих сетчатые оболочки из углепластика. Приведен краткий анализ существующих подходов к моделированию сетчатых конструкций. Для достоверного описания сложного поведения анизогридных конструкций при воздействии различных нагрузок предложены математическая и вычислительная модели. Высокая степень точности и устойчивости вычислительной модели, основанной на разложениях неизвестных функций по базису Фурье и базису, состоящему из полиномов Чебы-шева, обусловлена отсутствием насыщения таких методов приближения. Эффективность предложенных моделей и методов показана на примере решения тестовых краевых задач и задачи осевого сжатия анизогридной цилиндрической оболочки.
Ключевые слова: анизогридная конструкция, цилиндрическая оболочка, углепластик, континуальная модель, схема без насыщения, базис Фурье, полином Чебышева.
Введение. Прогресс в авиационной, ракетно-космической, судо-и машиностроительной отраслях промышленности в значительной степени связан с разработкой и применением новых композиционных материалов, обладающих повышенными физико-механическими характеристиками. Одним из перспективных решений может стать применение гибридных и анизогридных углепластиковых сетчатых конструкций [1].
Для решения задач расчета напряженно-деформированного состояния (НДС) сетчатых оболочек из углепластика, анализа потери их устойчивости и оптимизации параметров конструкций целесообразно использовать современные методы и средства математического и численного моделирования. При этом важно учитывать специфические особенности рассматриваемых конструкций: разносопротивляе-мость и нелинейность деформирования полимерных композиционных материалов при растяжении, сжатии и сдвиге [2], а также наличие малых геометрических и механических параметров и сложные условия эксплуатации с различными видами закрепления и на-гружения конструкции, воздействие высоких и низких температур, влагонасыщения, радиации.
Для достоверного описания поведения анизогридных оболочек с учетом указанных факторов целесообразно исходить из уравнений
пространственной теории упругости. Численное решение краевых задач для таких уравнений, описывающих НДС сетчатых оболочек, представляет существенную сложность и требует конструирования новых вычислительных алгоритмов, обладающих повышенной точностью и устойчивостью.
В настоящей работе предложены математическая и вычислительная модели для расчета НДС анизогридной цилиндрической оболочки при осевом сжатии. В основу математической модели деформирования положены линейные уравнения механики деформируемого твердого тела и континуальный подход, позволяющий поставить задачу в перемещениях для сплошного аналога сетчатой конструкции [2, 3]. В основе вычислительной модели лежат идеи нелокальных аппроксимаций, схем без насыщения и спектральных методов, позволяющие достичь точности наилучших приближений для любых классов гладкости искомых решений [4-6].
О моделировании поведения сетчатых конструкций. Формально сетчатую упругую конструкцию можно рассматривать как единое упругое тело. С позиций механики сплошной среды поведение такого тела может быть описано системой алгебродифференци-альных уравнений, включающей уравнения движения (равновесия) тела, кинематические и физические соотношения. Рассмотрение сетчатой конструкции как единого многосвязного тела представляется малоперспективным из-за наличия большого количества параметров. В связи с этим при моделировании сетчатых конструкций принято использовать один из двух основополагающих подходов: дискретный либо континуальный [3].
Дискретный подход. В рамках данного подхода конструкция представляется в виде совокупности стержневых, балочных или плоских элементов. Для описания взаимодействия этих элементов применяют классические метод сил, метод перемещений и смешанный метод строительной механики.
Другие дискретные модели, ориентированные на конструкции регулярной структуры, используют метод «склейки» [7-9]. Суть данного метода состоит в разбиении конструкции на изолированные элементы, анализе поведения каждого элемента в отдельности и использовании геометрических условий сопряжения элементов в конструкции при действии внешних сил и заданных перемещений. В контексте метода «склейки» можно выделить идею микроподхода, когда конструкция разбивается на элементы минимальной протяженности — ребра и узлы их стыковки, и идею макроподхода, применимую при наличии элементов, простирающихся на всю конструкцию. Метод непрерывной автоматической намотки, используемый для изготовления анизогридных оболочек, предполагает, что конструкция будет состоять именно из таких макроэлементов, поэто-
му макроподход может оказаться весьма перспективным для моделирования рассматриваемых оболочек.
Идеи дискретного подхода могут применяться для расчета сетчатых структур на ЭВМ, в частности, с использованием метода конечных элементов (МКЭ) [10]. При этом в рамках каждого ребра конструкции делается независимое приближение искомых функций посредством многочленов малой степени и для поиска решения записывается система линейных алгебраических уравнений.
Континуальный подход. Суть данного подхода состоит в замене реальной сетчатой конструкции эквивалентным сплошным анизотропным телом. Отметим, что данная замена является неоднозначной как с точки зрения выбора континуального эквивалента, так и с точки зрения задания способов перехода к нему. Причина неоднозначности кроется в более сложном механическом поведении сетчатой конструкции, нежели континуального тела в рамках классической теории упругости, в утрате тех или иных деформативных свойств конструкции.
Методы математической континуализации основаны на анализе разрешающих уравнений, которые могут быть получены, например, в рамках дискретного подхода. Суть данных методов состоит в переходе в разрешающих уравнениях от дискретных операций к непрерывным (например, от конечных разностей к производным).
Отметим две группы моделей, основанных на идее математической континуализации. К первой группе относят модели микрополярных континуумов [11, 12], представляющие конструкцию в виде множества ячеек, для которых записываются уравнения механики. Эти уравнения связывают характеристики соседних ячеек и в совокупности являются аналогами конечно-разностных соотношений, заданных на всей конструкции. По ним восстанавливается дифференциальная задача с непрерывными характеристиками. Вторая группа включает модели осреднения периодических структур [13-17], основанные на разложении решений в ряды по малому параметру — размеру ячейки, многократно повторяющейся в структуре конструкции. Такие модели позволяют с помощью точных периодических решений выделять асимптотически усредненное поведение системы. Методы математической континуализации, несмотря на высокую строгость и обоснованность, редко удается применить для решения прикладных задач. Математический аппарат данных методов эффективен, как правило, только в случаях модельных задач, когда материал конструкции имеет простую структуру и свойства, а сама конструкция находится под действием обычного нагружения.
Методы физической континуализации базируются на предположении о характере связи дискретных деформаций и внутренних сил конструкции с соответствующими непрерывными распределениями
деформаций и сил в ее сплошном аналоге. По сути, задача физической континуализации состоит в выводе физических соотношений для сплошного аналога сетчатой конструкции исходя из ее структуры и свойств материала, причем полученные соотношения должны обеспечивать максимальное соответствие поведения конструкции и ее аналога. Вывод физических соотношений можно осуществить с помощью энергетических принципов [18] или посредством усреднения характеристик сетчатой структуры. В последнем случае речь идет о так называемой концепции «размазывания», когда каждое семейство однонаправленных ребер структуры заменяется упругим слоем со специфическими свойствами [2]. Физические соотношения задаются исходя из условия эквивалентности внутренних сил в семействе ребер и тождественном слое. Отметим, что идея усреднения характеристик материала по объему прошла апробацию при построении структурных моделей механики композитов, учитывающих различные виды армирования [19, 20].
В данной работе рассматривается континуальная модель анизо-гридной конструкции, основанная на концепции «размазывания». Такой подход позволяет достаточно эффективно строить модели сетчатых конструкций и находить приближенные решения соответствующих прикладных задач (некоторые приложения описаны в [21]).
Континуальная модель анизогридной конструкции. Используем идею концепции «размазывания» для описания поведения анизогридной оболочки в рамках пространственной теории упругости. В [2, 22] описан стандартный подход к записи физических соотношений на основе концепции «размазывания» для тонких слоистых и сетчатых пластин. Такой подход основан на предположении о малости главных компонент тензоров напряжений и деформаций по направлению, ортогональному к плоскости пластины. При этом делают дополнительные предположения о том, что ребра работают только в продольном направлении (нитяная модель), и деформациями сдвига можно пренебречь. Данные предположения при анализе деформирования анизогридных оболочек представляются недостаточно обоснованными по следующим причинам:
1) в условиях сложного нагружения ребра сетчатой пластины будут работать не только на растяжение-сжатие, но также на изгиб, причем не только в плоскости пластины, но и в направлении нормали к ее поверхности. При изгибе в ребрах неизбежно возникнут деформации сдвига;
2) распределение напряжений по толщине ребра из углепластика может существенно отразиться на картине деформирования. Углепластик — композиционный материал на полимерной основе, обладающий сложным нелинейным поведением, которое проявляется при растяжении, сжатии и сдвиге.
В статье предлагается использовать уточненные соотношения, полученные на основе концепции «размазывания» и позволяющие рассчитывать НДС сетчатых конструкций из углепластиков в пространственной постановке. Будем считать углепластик ортотропным квазиоднородным материалом с физическими соотношениями, записанными в виде обобщенного закона Гука, связывающего деформации е1з 82, 83, у12, у13, у23 с напряжениями оь о2, о3, т12, т13, т23 :
Здесь Е, О — модули упругости и сдвига в направлении волокон; ЕТ, ОТ — модули упругости и сдвига в направлениях 2, 3; V характеризует влияние напряжений, возникающих в направлениях 2, 3, на деформации в направлении 1; Vт — влияние напряжений в направлении 1 на деформации по осям 2, 3; V ь — влияние напряжения в направлении 2 на деформацию по оси 3 и обратно. Соотношения (1) записаны в предположении, что свойства ребра в направлениях 2, 3 тождественны и отличны от свойств ребра в направлении волокон 1 (рис. 1, а). Далее используем выражения напряжений через деформации в виде
= £ау8у, 1 = 1,-,э; т12 = °У12; т13 = °У13; т23 = °у23. (2)
(1)
= т13 = т23
У13 = "77"; У 23 = •
3
У=1
а
б
Рис. 1. Ортотропное ребро из углепластика:
а — в системе координат 1, 2, 3, связанной с ребром; б — в системе координат х, у, г, связанной с конструкцией
Здесь коэффициенты aij выражаются из (1) по формулам
a11 = (1 — vT )E / д; a12 = a13 = vE / д;
(1 — vvT ) ET
a21 = a 31 = vLET / д; a22 = a33 =-л 1 * ; (3)
Д(1 + vT)
a23 = a32 = (vt ++vvL)Et ; д = 1 — 2vvL —vT > 0.
д(1 + vT)
Уравнения (2) записаны в системе координат 1, 2, 3, связанной с ребром (см. рис. 1, а). Для вывода физических соотношений для сплошного аналога сетчатой конструкции введем систему координат x, y, z, связанную с конструкцией (рис. 1, б), и перейдем от системы координат 1, 2, 3 к системе x, y, z с помощью операции поворота на угол ф. Пусть s = sin ф, c = cos ф, тогда напряжения и деформации преобразуются по следующим формулам:
ox = a1c2 + o2 s2 — 2t12cs; sx = s1c2 +s2 s2 — y12cs;
Oy = СТ^2 +<2c 2 + 2t12cs; < =<3; Sy =s1s2 +s2c2 +Y12cs; Sz =S3; (4)
Txy = — (<1 — <2)cs + T12(c 2 — s 2); Yxy = 2(s1 — ^2)cs + Y12(c 2 — s 2); Txz = T13c — T23s; Tyz = T13s + T23c; Yxz = Y13c — Y23s; 1 yz = Y13s + Y23c
Используя (2)-(4), запишем физические уравнения для ребра в системе x, y, z:
<x = A11Sx + A12sy + A13sz + A14Yxy; Txy = A41sx + A42sy + A43sz + A44Yxy; < y = A21S x + A22S y + A23S z + A24 Y xy; T xz = A55 Y xz + A56 Y yz; (5)
<z = A31Sx + A32Sy + A33Sz + A34Yxy; Tyz = 4s5Yxz + АбYyz.
Коэффициенты Ak¡ имеют следующий вид:
A11 = a11c4 + (a21 + a12)^c2 + a22^ + 4G^s2; n ОЛ
(1,2) (6)
9 9 A A 99 \ j / \ /
A12 = (a11 + a22)s c + a21s + a12c — 4Gc s ,
3 3 3 3 9 9
A14 = a11c s + a21cs — a12c s — a22cs — 2G (c — s )cs;
3 3 3 3 2 2
A24 = a11cs + a21c s — a12cs — a22c s + 2G (c — s )cs,
(12,21), (14,41), (24,42) (7)
2 2
^23 = a13s + a23c ; A43 = cs(a13 — a23),
(13,31), (23,32), (43,34) (8)
2 2 2 2 A13 = a13c + a23s ; A23 = a13s + a23c ; A43 = cs(a13 — a23)
A33 = a33; A44 = (a11 - a21 - a12 + a22 )c2^ + G (c2 - s2 )2; (9)
A65 = A56 = (G-GT)cs; A55 = Gc2 + GTs2 (13,23),(55,66),(G,GT). (10)
Здесь числа и обозначения в скобках обозначают циклические перестановки индексов, например: (1, 2) означает, что если в соответствующих уравнениях заменить индексы 1 на 2, 2 на 1, то получатся верные равенства.
Сетчатая конструкция состоит из семейств ребер, каждое из которых имеет номер j (j = 1,..., N) и определяется углом наклона ф.
оси 1 ребра к оси х (рис. 2), толщиной ребер 5. и расстоянием aj
между ними. В соответствии с концепцией «размазывания» [2, 3] заменим каждое семейство ребер сплошной поверхностью с физическими соотношениями вида (5). Вместо коэффициентов Ak¡ в физических соотношениях слоя j-го семейства возьмем Aj. Для записи AJ обозначим символами aJkl, EJ, ET, vj, vJL, vJT, GJ, GT, s. = sin ф., Cj = cos ф. соответствующие характеристики ребер j-го семейства. Соотношения (6)—(10) после введения индекса j не изменятся. Коэффициенты aj
в них вычисляются по формулам, аналогичным (3), в которых характеристики материала ребер усреднены по поверхности конструкции, т. е. имеют вид
Ej = kJEJ; ET = kJET; vJ = kJvJ; vj = kJvj; vL = kJvJL;
. . . . . 5.
GJ = kJGJ; GT = kJGT; kJ = -..
Рис. 2. Замена ребер j-го семейства эквивалентным сплошным слоем
В таком случае равные массовые силы, действующие в семействе ребер и соответствующей сплошной поверхности, будут приводить к эквивалентным деформациям. В указанном смысле метод расчета коэффициентов Луц корректен для каждого семейства ребер. Связь напряжений а х, а у , а г, т^, т ж, ту2 и деформаций в х, ву , в г, у ^, уХ2, Уу2 в континуальном эквиваленте анизогридной конструкции имеет вид (5), где
N
Лк1 = Х Лк!. (11)
у=1
Поскольку построенный континуальный аналог анизогридной конструкции имеет эквивалентные деформативные свойства, дальнейшие рассуждения будем проводить, отождествляя сетчатую конструкцию и ее сплошной аналог.
Вывод разрешающих уравнений и постановка задачи деформирования сетчатой цилиндрической оболочки. Для моделирования процесса деформирования цилиндрической конструкции воспользуемся соотношениями пространственной теории упругости. При этом будем исходить из линейных кинематических и физических уравнений, вполне удовлетворительно описывающих поведение тела при малых нагрузках и перемещениях [21].
Осуществив переход в цилиндрическую систему координат (г, ф, г), где 0 < г < /; 0 <ф<2п; Я - И /2< г < Я + И /2; Я — радиус срединной поверхности конструкции; / — ее высота (рис. 3), приходим к задаче в канонической области О. Запишем для конструкции уравнения равновесия и кинематические соотношения:
д (гаг)-аф+^Ъ^ + (гт2г) + ¥гг = 0; дг дф дг
даф д д
~дф~ + дк (гТгф ) + дг (гТгф ) + Тгф + ^рг = 0; (12)
1Т (га г) + 1Г (г Тгг) + ^ + ^г = 0;
дг дг дф
ди 1 ду и с№ 1
В г — ; Вф - I ; В г — •
дг г дф г дг '
(13)
1 ди V | ду 1 дw дw ди
У гф — ^ +г I I; У фг — ^ + ^ ; У гг — ^ + ^ , г дф дг ^ г) дг г дф дг дг
где и(г, ф, г), у(г, ф, г), ^(г, ф, г) — перемещения точек тела по направлениям г, ф, г соответственно; в г, Вф, в г, у гф, Уф2, у ^, аг, аф,
аг, тгф, тфг, т^ — деформации и напряжения в континуальном аналоге сетчатой конструкции; Гг, , ^ — внутренние силы.
\в
ч = (Яг, «V Чт)
Рис. 3. Замена цилиндрической сетчатой конструкции сплошным анизотропным телом и переход к канонической области О
Физические соотношения (5) в системе координат г, ф, 2 примут следующий вид:
аг = А118г + А128ф + А138г + А14Угф; Хгф = А418г + А428ф + А438г + А44Угф; аф = А218г + А228ф + А238г + А24УТгг = А5зУгг + А5бУфг; (14)
аг = А318 г + А328ф + А338г + А34 У гф; тфг = Аб5У гг + Абб Уфг •
Сделаем преобразование растяжения по осям г, г: г ^ г = п(г - Я + 0,5И)/ И; г ^ г = пг //; 0 <г, г <п, и соответствующие замены переменных в уравнениях (12), (13) (далее знак тильды над г и г опускаем). Подставляя (13) в (14), (14) в (12) и полагая, что ¥х = ¥у = ^ = 0, получаем уравнения для перемещений точек
в области О:
Л, и = А
2 а2
п д и
Ми = А33 И2 дг2 + 4)6
д2 2
д и . п д и
—^ + А,
дф2
155 ,2
/2 дг2
= 11 - Л I Абб
д 2 и
дф2
/Л Л \П2 д2^ 1
-(А31 + А55)^^~ + -
И/ дгдг г
п
д^
(
Л А21 А31^ +(А22 А32)
/ дг г
ду дф
- + и
п
ди
+ , (А23 А33) п +(А24 А34)
И дг
' п ду + 1 д^ ^ / дг г дф
-1А
п д2у п ди 1 ду 1 ^
32
+--------и
И дфдг И дг г дф г
у
- Л
V д 2у
34
1 дч 1 д2 ч Л
- + —
И/ дгдг г дф г дфдг
Л65
г
Гп д2ч п д2и } 1 , (п д2
И дфдг / дфдг
-1 Л.
Л66
г
1 дv
Л
56
п д2и
- + г-
2 2 п 5 V п ду
И дфдг г дф
/ дфдг И/ дгдг / дг
(15)
л -л п2 5 2У д \ п2 д \ = ( 1 V д2у = 4)6 И2 дг2 + Л22 дф2 + Л44 /2 дг2 = I1 - г21Л22 дф2
2 я2
п д и
- Л43 , , ^ ~ - Л44
И/ дгдг
Л65^1 + Л661 Ж2
г
-Л -Л42 /г
^ д2у ди
дфдг дг
п д2 ч - 2 /г дфдг г
^ п ( п д2ч п д2и ^ Л55 2 +
И ^ И дг2 / дгдг
- Л,
2 я2
п д ч 41д^
1 ди п
)
д 2 и
)
+ + Лл^
г дф И дгдф
+ + Л
2
п д V
24
2
1 д ч
- + —
/ дгдф г дф2
п д2ч
Л21 / дгдф
-1 ЛМ;
г
(16)
2 я2 я2
. п д ч .о ч Л3ч = Л55 ТГ^Т +
-Л -
Л12 , 1г
(
И2 дг2
д2v ди +
дф2
+Л _(1 О Л
+ Л11 /2 дг2 Ч г2 ) Л44 дф2
дфдг дг
2 2
. п д ч . п
- Л13 и дг2 - Л14 у
^пд2у д 2 чЛ
п2 д2и 1 .
-----Л
И/ дгдг г
56
2 я2
п г д V 3 И2 дг2
41
д 2 ч дфдг
2
+ Л42
д V ди + ■
дф2 дф
/ дг2 дфдг
-1 [ Л55^1 + Л56^2 ]-
2 2 . п д и пд V
+ Л43 Т - - + Л44 '
И дгдф
/ дфдг
п д2и
(17) п дv
п дч п ди ди п дv
Здесь ^ =--+--; Ж2 =--V +—г —; Ж3 -
И дг / дг дф И дг И дгдф И дг
1 ди 1
---+ -V.
г дф г
Будем искать решения системы уравнений (15)—(17) в классе достаточно гладких 2п -периодических по координате ф функций. Рассматривая задачу об осевом сжатии цилиндрической конструкции под действием нагрузок р = (0, 0, рг), q = 0, выпишем граничные условия при г - 0, п:
=0,п = Рг =>
с№
/
дг пА
11
А12 Г ду
дф
+ и
^ А13п ди И дг
- А
14
Гп ду + 1 д^ ^ / дг г дфу
ду
1г=0
п= Рф= 0 => ^Г=--
дг
п
А41 п д^ А42 1 Г ду ^ 41 — + и
А44 1 дг А44 г
дф
У
А43п ди дw
--+ —
А44И дг дф
ди
1гф \г=0,п = рг = 0 => - =
дг п
А5б Г1 ди п ду 1 ^ п дw
5б--+----у +--
И дг
55
г дф И дг г у
(18)
и граничные условия при г = 0, п :
. . ди И
аг \г=0,п = Чг = 0 => -Т- = --
дг
пА
33
п
Г
А31 , - + А32
/ дг г
ду дф
- + и
+А
-34
г п ду + 1 д№ ^ / дг г дфу
ду
И
Сга \г=0,п = Чф = 0 => п =
дг п
И
1гг \г=0,п = Чг = 0 => - =
дг п
Аз5 гп5^ + п ди^ + 1 ди 1 у Абб ^ И дг / дг У г дф г
А5б Г1 ди п ду 1 ^ п ди
5б--+----у +--
/ дг
55
г дф И дг г у
(19)
Задачу (15)—(19) можно переписать в более компактном виде:
V = /1(и, у, w); Л2у = /2 (и, у, w); Л3w = /3(и, у, w); г, г е [0, п], фе [0,2п], и, у, w - 2п-периодические функции по ф;
^ = gl(u,уw); ^ = g2(u,уw); ^т = gз(u,уw) пРи г = 0п; (20)
дг дг дг
^ = ^4(иу,w); ^ = g5(u,у,w); ^ = g6(u,у,^ пРи г = 0п. дг дг дг
Здесь и, у, w — 2п-периодические функции по ф; / - /3, g1 - g6 — функции, включающие операторы первых и вторых производных, вид которых определяют равенства (15)-(19).
Разработка численного метода на основе приближений без насыщения. Система уравнений в частных производных (20) содер-
2
жит малые параметры: геометрические величины (И / п) , 5 / а■ и
механические характеристики жесткости углепластика Ет / Е, входящие в коэффициенты системы. Кроме того, отметим, что граничные условия записаны в неявной форме с использованием функций gl - g6. Поставленная задача имеет высокую вычислительную сложность и требует создания оригинальных численных методов, позволяющих строить устойчивые алгоритмы и находить решения и, V, ч с достаточной точностью за приемлемое время.
Для построения таких алгоритмов воспользуемся идеями К.И. Ба-бенко о методах приближения без насыщения [4]. За основу возьмем алгоритм решения 2В-краевых задач для уравнения Пуассона и его развитие и обобщение для случая произвольной размерности [6, 23]. Данный алгоритм использует итерационный метод установления и набор нестационарных регуляризаций, что делает его весьма эффективным при решении как линейных, так и нелинейных жестких краевых задач. В основе алгоритма лежат приближения неизвестных функций интерполяционными полиномами с узлами Чебышева, не имеющие насыщения:
Здесь /(х) — приближаемая неизвестная функция, х е [-1,1]; хт —
узлы интерполяции; М — их количество. При расчетах будем использовать модификацию (21) для приближения решений по направлениям г, 2, позволяющую задать значения решения и его производных в граничных точках /(±1), /'(±1):
Рм(х /) =
соб (М агееоБ х);
(21)
Рм (х, /) = 7тЕ Рт (х)/(Хт ) + 4 (х)/ '(1) + (х)/ '(-1) +
1 М
т=1
(22)
Тм (х); Тм (х) = соб (М агееоБ х);
х) = Тм(х)(1 + х)2 {1 + (1 -х)(1 + М2)}; л (х) = Тм (х)(1 - х2)(1 - х) , v (х) = Тм(х)(1 - х2)(1 + х)
v-l(x) = 4(-)м ; Vl(x) = 4 ,
при этом к узлам хт, т = 1, ..., М, добавляют точки -1, 1 и преобразуют отрезок [-1, 1] в отрезок [0, п], так как в (20) г, г е [0, п]. Для приближения неизвестных 2п-периодических функций по координате ф используют полином с ядром Дирихле, который может быть получен из ряда Фурье элементарными преобразованиями:
2 2 N-2
Рф (х, /) = т—Г Е / (х1)DN-1(х - хЛ);
2N-1 ■=о (23)
х = 2\Л7, ■ = 0,1,...,2N-2,
где DN-1 (х) = —°5)х — ядро Дирихле. Аппроксимации функ-2б1п( х / 2)
ций и, v, ч представляют собой произведения двух полиномов (22) и полинома (23) Рт, где М — количество узлов интерполяции в трехмерной области О. Такие приближения на гладких функциях обладают асимптотикой погрешности наилучших приближений [4]:
II/ - РМDCaМ,
где а — степень гладкости /; D — константа, ограничивающая значения производных /; Са зависит только от степени гладкости / Данные приближения существенно минимизируют погрешность расчетов и объем памяти, необходимый для работы программы, позволяя справиться с проблемой малых параметров и получить достаточно точное решение с минимальными временными затратами.
Будем использовать для решения задачи (20) вычислительный алгоритм, разработанный в [6]. Он состоит в последовательном выполнении следующих операций:
1) введение дополнительной переменной I и запись регуляризо-ванных уравнений (20). В случае использования параболической регуляризации данные уравнения включают дополнительные производные по I и имеют вид
Щ = А^ - /1 (и, V, ч); vt = А^ - /2 (и, V, ч); ^ = Аз^ - /3 (и, V, ч); (24)
2) поиск решений исходной системы как предела решений (24) при t ^ да с помощью метода установления. Обоснование сходимо-
сти метода установления при использовании двух регуляризаций для
некоторых линейных задач дано в [23];
3) аппроксимация производных по t разностными отношениями, производных по г, ф, г — производными от полиномов (22), (23);
4) запись задач линейной алгебры в виде аппроксимирующих уравнений, включающих трехмерные массивы, и решение этих задач с помощью спектрального разложения матриц, аппроксимирующих производные [6];
5) переход от предыдущей итерации метода установления к следующей. Такие переходы осуществляются до тех пор, пока решения не установятся, т. е. пока норма разности решений, полученных на следующей и предыдущей итерации, не станет достаточно близкой к нулю.
Тестовые численные эксперименты. Для апробации предложенной вычислительной схемы и демонстрации ее высокой точности и устойчивости в задачах с малыми параметрами и граничными условиями, записанными в неявном виде, рассмотрим две тестовые задачи. В качестве первой возьмем уравнение Пуассона для функции и (х, у, г)
Ли =
д2и д2и д2и
дх2 ду2 дг2
2(х2 - пх) + 2(у2 - пу) + (х2 - пх)(у2 - пу)
ег,
х, у, г е [0, п] (25)
с граничными условиями неявного вида при г = 0, п
и = 0 при х, у = 0, п; иг = и при г = 0, п (2б)
и известным точным решением иех (х, у, г) = (х2 - пх)(у2 - пу)ег. Найдем с помощью предложенного метода численное решение задачи (25Х (2б) иар (х, у], гк), I = 1,..., N ] = 1,..., М; к = 1,..., К, при различных количествах узлов интерполяции Ы, М, К и приведем график десятичного логарифма погрешности 8= иех -иар /||иех|(рис. 4).
Здесь полиномы вида (21) применялись для приближения функции по каждому направлению.
Видно, что всего при десяти узлах сетки в каждом из направлений получаем приближенное решение, отличающееся от точного только в десятом знаке.
В качестве второго эксперимента рассмотрим задачу о расчете прогиба тонкой прямоугольной пластины, защемленной по всему контуру. В рамках классической теории пластин прогиб w(x, у) такой пластины описывается бигармоническим уравнением
Рис. 4. Результаты решения первой тестовой задачи:
а — приближенное решение на сетке с узлами Чебышева; б — десятичный логарифм погрешности
ЛЛw = ч(х, у)/Б; Б = ЕИ3 /12(1 -V2), х е [0, а], у е [0, Ь] (27)
с соответствующими граничными условиями
dw dw
— = 0, w = 0 при х = 0, a; — = 0, w = 0 при y = 0, b, (28)
дх dy
где D — изгибная жесткость; E, v — модуль Юнга и коэффициент Пуассона соответственно; q(x, y) — нагрузка, действующая на плоскость пластины [24]. Рассмотрим случай a = b = 1, q(x, y) = q = const. С помощью представления решения в виде тригонометрического ряда и удержания в нем первых четырех членов можно найти максимальное значение прогиба пластинки с погрешностью не более 1 %: w.
2 D
= 0,00126 q / D.
Рассмотрим теперь соответствующую задачу в трехмерной постановке на основе соотношений теории упругости:
(^ + 2|д)—- + |
д2и д2и | д2и .
ах- - 2 ^ =^
^д 2v
дy2 h дг2
1 д2 w^
дxдy h дхдг
д 2v
+ (^ + 2|)
дх
д2v | д2v
дy h2 дг2
^д 2u
1 д2v ^
дxдy h дyдz
(29)
|
д2 w дх2
+ ^ д2w + (Л, + 2|) д2w _ +
2
дy 2
h2
дг2
h
д u
д2у ^
дхдг аyдz
Здесь Х = уЕ / ((1 + у)(1 - 2у)), д = Е / (2(1 + V)) — параметры Ламе; и(х, у, г), у(х, у, г), х, у, г) — перемещения точек пластинки, х е [0,а], у е [0,6], г е [-И, 0]. Граничные условия защемления имеют вид
и, V, ч = 0 при х = 0, а, у = 0, 6; (30)
на границах г = 0, -И имеем
ди / дг = -дч / дх; ^ / дг = -дч / ду при г = 0, -И;
дч д
ХИ
(
ди дv
(Х + 2 д) ^ дх ду )
при г = -И;
(31)
ХИ
(
ди дv
Л
— = q( х, у) -дг (Х + 2д) ^ дх ду)
при г = 0.
Отметим, что основные вычислительные сложности задачи (20)
-3
имеют место и в задаче (29)-(31). Пусть а = 6 = 1, q (х, у) = q = 10 . На рис. 5, а приведена относительная разностьч2П и ч3П.
0,05 0,10 0,15 А, м а
Рис. 5. Относительная разность и (а) и распределение прогиба х, у, - И) при И = 0,01 м в задаче (29)-(31) (б)
Виден характер поведения этой разности при уменьшении И: вблизи нуля она принимает значения порядка 1 %. Данный факт свидетельствует о корректности постановок задач на основе пространственной теории упругости, приведенных в настоящей работе, и эффективности описанного алгоритма. На рис. 5, б приведено распределение прогиба при г = -И, найденное при решении задачи (29)-(31). На рис. 6, а изображено распределение производной от
прогиба ж Видно, что на границах х = у = 0, х = у = 1 производные по х и по у принимают нулевые значения (жирная линия), что является еще одним подтверждением соответствия постановок задач в двумерном и трехмерном случае (см. (28)). Рис. 6, б демонстрирует распределение прогиба в сечении осью у = л/2. Видно, что деформация в направлении г близка к нулю.
Рис. 6. Результаты решения задачи (29)-(31) при И = 0,01 м: а — производная м>, принимающая нулевые значения на границах х = у = 0, х = у = 1; б — распределение ^ в сечении плоскостью у = п/2
Решение задачи об осевом сжатии сетчатой цилиндрической конструкции. Приступим к решению задачи (20) при нагрузках, соответствующих осевому сжатию конструкции. Рассмотрим конструкцию, включающую одно семейство кольцевых ребер Ф1 = 0 и два семейства
спиральных ребер ф23 =±26,6° с параметрами £1,2,3 = 140 ГПа,
Е1Т2,3 = 11 ГПа, а1,2,3 = 5,5 ГПа, О1/'3 = 5,386 ГПа, у^2'3 =0,27, V1,2,3 =
12 3
= vf ' = 0,0212. Определим сжимающие напряжения осж, возникающие в спиральных ребрах, и осевую жесткость конструкции £ при Я = 1 м, I = 4 м, Р = рг = 2 МН, Q = 0, И = 1,05 см и различных значениях ширины спиральных и кольцевых ребер 5 и , 5С и расстояний между ними а, аи, ас (рис. 7). В таблице данные значения отмечены звездочкой. Они приведены в сравнении с величинами, полученными в [10] при использовании дискретной модели (индекс «д») и упрощенной двумерной континуальной модели (индекс «к») рассматриваемой сетчатой оболочки при тех же Я, I, Р, Q. В поставленной задаче содержатся
—7
малые параметры порядка 10 при старших производных и, v, м и величины порядка 10-5 в граничных условиях уравнений (20), что обусловливает ее высокую вычислительную сложность.
Рис. 7. Основные параметры модели
Полученные решения и результаты работы [10]
Основные параметры модели Напряжения в спиральных ребрах Значения изгибной жесткости
а, мм 5И, мм 8С, мм °сж.д, МПа [10] ^сж.к? МПа [10] ^сж? МПа* 5д, МН/мм [10] МН/мм [10] МН/мм*
130,8 8,1 4,5 273,4 273 275,5 0,115 0,111 0,137
104,7 6,5 3,6 272 273 274,8 0,114 0,111 0,138
87,2 5,4 3,0 273,5 273 275,5 0,115 0,111 0,137
78,5 4,86 2,7 273,5 273 275,6 0,116 0,111 0,137
69,8 4,3 2,4 274 273 276,8 0,114 0,111 0,137
При сопоставлении данных таблицы видно, что решения, полученные при использовании принципиально различных математических и вычислительных моделей, достаточно близки. Это является косвенным подтверждением корректности полученного результата.
Заключение. Предложенные математическая и вычислительная модели позволили рассчитать НДС анизогридной цилиндрической оболочки на основании уравнений пространственной теории упругости. В случае изотропного и ортотропного материалов с линейной связью напряжений и деформаций полученные результаты близки
к решениям классической теории пластин и дискретных, и континуальных моделей сетчатых оболочек.
Возможность достижения высокой точности решений при небольшом количестве узлов интерполяции является весьма перспективной для анализа устойчивости конструкций и оптимизации их параметров, поскольку данные задачи зачастую сводятся к решению большого числа прямых задач определения НДС. Минимизируя количество узлов с сохранением высокой точности, сводим к минимуму время решения каждой такой задачи, следовательно, получаем значительное преимущество в скорости при высокой точности расчетов.
Работа выполнена при поддержке РФФИ (грант №13-01-12032 офи_м_2013).
ЛИТЕРАТУРА
[1] Васильев В.В., Барынин В.А., Разин А.Ф., Петроковский С.А., Халиманович В.И. Анизогридные композитные сетчатые конструкции — разработка и применение в космической технике. Композиты и наноструктуры, 2009, № 3, с. 38-50.
[2] Vasiliev V.V., Morozov E.V. Advanced Mechanics of Composite Materials. Ann Arbor, Elsevier, 2007, 491 p.
[3] Образцов И.Ф., Рыбаков Л.С., Мишустин И.В. О методах анализа деформирования стержневых упругих систем регулярной структуры. Механика композиционных материалов и конструкций, 1996, т. 2, № 2, с. 3-14.
[4] Бабенко К.И. Основы численного анализа. Москва; Ижевск: НИЦ «Регулярная и хаотическая динамика», 2002.
[5] Boyd J. Chebyshev and Fourier Spectral Methods. 2d edition. Ann Arbor, University of Michigan, 2000.
[6] Семисалов Б.В. Нелокальный алгоритм решения уравнения Пуассона и его приложения. Вычислительная математика и математическая физика, 2014, т. 54, № 7, с. 1110-1135.
[7] Левин М.А. Некоторые задачи о регулярных стержневых системах. Изв. вузов. Строительство и архитектура, 1965, № 9, с. 41-48.
[8] Рыбаков Л.С. О теории одной плоской регулярной упругой структуры ферменного типа. Механика твердого тела, 1995, № 5, с. 171-179.
[9] Dean D.L., Ganga Rao H.V.S. Macro approach to discrete field analysis. J. Eng. Mech. Div, ASCE, 1970, vol. 96, no. EM 4, pp. 377-394.
[10] Азаров А.В. Континуальные и дискретные модели сетчатых композитных цилиндрических оболочек. Механика композиционных материалов и конструкций, 2012, т. 18, № 1, с. 121-130.
[11] Bazant Z.P., Christensen M. Analogy between micropolar continuum and grid frameworks under initial stress. Int. J. Solids and St., 1972, vol. 8, no. 3, pp. 327-346.
[12] Бунаков В.А., Протасов В.Д. Сетчатые композитные цилиндрические оболочки. Механика композиционных материалов, 1989, № 6, с. 1046-1053.
[13] Бахвалов Н.С., Панасенко Г.П. Осреднение процессов в периодических средах. Математические задачи механики композиционных материалов. Москва, Наука, 1984, 352 с.
[14] Власов А.Н. Усреднение механических свойств структурно-неоднородных сред. Механика композиционных материалов и конструкций, 2004, т. 10, № 3, с. 424-441.
[15] Димитриенко Ю.И., Губарева Е.А., Сборщиков С.В. Асимптотическая теория конструктивно-ортотропных пластин с двухпериодической структурой. Математическое моделирование и численные методы, 2014, № 1, с. 36-57.
[16] Димитриенко Ю.И., Губарева Е.А., Сборщиков С.В. Конечно-элементное моделирование эффективных вязкоупругих свойств однонаправленных композиционных материалов. Математическое моделирование и численные методы, 2014, № 2, с. 28-48.
[17] Шешенин С.В., Скопцов К.А. Теория пластин, основанная на методе асимптотических разложений. Математическое моделирование и численные методы, 2014, № 2, с. 49-61.
[18] Алтуфов Н.А., Попов Б.Г. Континуальные модели регулярных ферменных конструкций. Механика твердого тела, 1994, № 6, с. 146-154.
[19] Митюшов Е.А. Теория армирования. Механика композиционных материалов и конструкций, 2000, т. 6, № 2, с. 151-161.
[20] Свистков А.Л., Евлампиева С.Е. Использование сглаживающего оператора осреднения для вычисления значений макроскопических параметров в структурно-неоднородных материалах. ПМТФ, 2003, т. 44, № 5, с. 151-161.
[21] Голушко С.К., Идимешев С.В., Семисалов Б.В. Методы решения краевых задач механики композитных пластин и оболочек: учеб. пособие по курсу «Прямые и обратные задачи механики композитов». Электрон., текстовые и граф. данные - КТИ ВТ СО РАН. Новосибирск, 2014, 131 с.
[22] Васильев В.В. Механика конструкций из композиционных материалов. Москва, Машиностроение, 1988, 269 с.
[23] Блохин А.М., Ибрагимова А.С., Семисалов Б.В. Конструирование вычислительного алгоритма для системы моментных уравнений, описывающих перенос заряда в полупроводниках. Математическое моделирование, 2009, т. 21, № 4, с. 15-34.
[24] Тимошенко С.П., Войновский-Кригер С. Пластинки и оболочки. Москва, Физматлит, 1966.
Статья поступила в редакцию 31.01.2015
Ссылку на эту статью просим оформлять следующим образом:
Голушко С.К., Семисалов Б.В. Численное моделирование деформирования анизогридных конструкций с применением высокоточных схем без насыщения. Математическое моделирование и численные методы, 2015, № 2(6), с. 23-45.
Голушко Сергей Кузьмич родился в 1958 г., окончил Новосибирский государственный университет в 1980 г. Д-р физ.-мат. наук, заведующий лабораторией анализа и оптимизации нелинейных систем Института вычислительных технологий СО РАН, директор Конструкторско-технологического института вычислительной техники СО РАН. Является экспертом Российского фонда фундаментальных исследований, Фонда содействия развитию малых форм предприятий в научно-технической сфере, Российского научного фонда. Автор более 170 научных публикаций в области математического моделирования, вычислительной механики, биомеханики, оптимального проектирования. е-mail: s.k.golushko@gmail.com
Семисалов Борис Владимирович родился в 1987 г., окончил Новосибирский государственный университет в 2010 г. Канд. физ.-мат. наук, научный сотрудник Конструкторско-технологического института вычислительной техники СО РАН. Победитель Всероссийского конкурса программы «У.М.Н.И.К.» (2013 г.). Автор и соавтор более 40 научных работ в области математического моделирования, вычислительных методов, нелинейных краевых задач, полупроводниковых устройств, композитов, гемодинамики. е-mail: vibis@ngs.ru
Numerical modeling of anisogrid structures deformation using schemes of high accuracy without saturation
© S.K. Golushko1' 2, B.V. Semisalov1
:Design Technological Institute of Digital Techniques, SB RAS, Novosibirsk, Russia 2Institute of Computational Technologies, SB RAS, Novosibirsk, Russia
The article describes a class of promising anisogrid structures representing mesh shell of unidirectional carbon. A brief analysis of existing approaches to modeling deformation of grid structures is presented. New mathematical and numerical models are proposed for reliable description of complex behavior of anisogrid structures under different kinds of loads. A high degree of accuracy and stability of the numerical model based on the expansions of unknown functions in Chebyshev polynomials and Fourier series is caused by the lack of saturation of such approximations. Efficiency of the proposed models and techniques is demonstrated on the example of solving test boundary-value problems and a problem of axial compression of anisogrid cylindrical shell.
Keywords: anisogrid structure, cylindrical shell, carbon, continuum model, scheme without saturation, Fourier series, Chebyshev polynomial
REFFERENCES
[1] Vasilyev V.V, Barynin V.A., Razin A.F. Petrokovskiy S.A., Halimanovich V.I.
Kompozity i nanostruktury - Composites and nanostructures, 2009, no. 3, pp. 38-50.
[2] Vasilyev V.V., Morozov E.V. Advanced Mechanics of Composite Materials. Elsevier, 2007, 491 p.
[3] Obraztsov I.F., Rybakov L.S., Mishustin I.V. Mekhanika kompozitsionnykh materialov i konstruktsiy - Mechanics of Composite Materials and Structures, 1996, vol. 2, no. 2, pp. 3-14.
[4] Babenko K.I. Osnovy chislennogo analiza [Fundamentals of Numerical Analysis]. Moscow, Izhevsk, SRC Regular and chaotic dynamics Publ., 2002.
[5] Boyd J. Chebyshev and Fourier Spectral Methods. Second edition, University of Michigan, 2000.
[6] Semisalov B.V. Zhurnal vychislitelnoy matematiki i matematicheskoi fiziki RAN - Journal of Computational Mathematics and Mathematical Physics, 2014, vol. 54, no. 7, pp. 1110-1135.
[7] Levin A. Izvestiya vuzov. Stroitelstvo i arkhitektura - Proceedings of the universities. Construction and architecture, 1965, no. 9, pp. 41-48.
[8] Rybakov L.S. Mekhanika tverdogo tela - Mechanics of Solids, 1995, no. 5, pp. 171-179.
[9] Dean D.L., Ganga Rao H.V.S. Macro approach to discrete field analysis. J. Eng. Mech. Div., ASCE, 1970, vol. 96, no. EM4, pp. 377-394.
[10] Azarov A.V. Mekhanika kompozitsionnykh materialov i konstruktsiy -Mechanics of Composite Materials and Structures, 2012, vol. 18, no. 1, pp. 121-130.
[11] Bazant Z.P., Christensen M. Analogy between micropolar continuum and grid frameworks under initial stress. Int. J. Solids and St. 1972, vol. 8, no. 3, pp. 327-346.
[12] Bunakov V.A., Protasov V.D. Mekhanika kompozitsionnykh materialov i konstruktsiy - Mechanics of Composite Materials and Structures, 1989, no. 6, pp. 1046-1053.
[13] Bakhvalov N.S., Panasenko G.P Osrednenie protsessov v peroidicheskikh sredakh. Matematicheskie zadachi mekhaniki kompozitsionnykh materialov [Averaging Processes in Periodic Media. Mathematical Problems of the Mechanics of Composite Materials]. Moscow, Nauka Publ., 1984, 352 p.
[14] Vlasov A.N. Mekhanika kompozitsionnykh materialov i konstruktsiy -Mechanics of Composite Materials and Structures, 2004, vol. 10, no. 3, pp. 424-441.
[15] Dimitrienko Yu.I., Gubareva E.A., Sborschhikov S.V. Matematicheskoe modelirovanie i chislennye menody - Mathematical modeling and Numerical Methods, 2014, no. 1, pp. 36-57.
[16] Dimitrienko Yu.I., Gubareva E.A., Sborschikov S.V. Matematicheskoe modelirovanie i chislennye menody - Mathematical modeling and Numerical Methods, 2014, no. 2, pp. 28-48.
[17] Sheshenin S.V., Skoptsov K.A. Matematicheskoe modelirovanie i chislennye menody - Mathematical modeling and Numerical Methods, 2014, no. 2, pp. 49-61.
[18] Altufov N.A., Popov B.G. Mekhanika tverdogo tela - Mechanics of Solids, 1994, no. 6, pp. 146-154.
[19] Mityushov E.A. Mekhanika kompozitsionnykh materialov i konstruktsiy - Mechanics of Composite Materials and Structures, 2000, vol. 6, no. 2, pp. 151-161.
[20] Svistkov A.L., Evlampieva S.E. Prikladnaya mekhanika i tekhnicheskaya fizika - Journal of Applied Mechanics and Technical Physics, 2003, vol. 44, no. 5, pp. 151-161.
[21] Golushko. S.K., Idimeshev S.V., Semisalov B.V. Metody resheniya kraevykh zadach mekhaniki kompozitnykh plastin i obolochek: uchebnoe posobie po kursu "Pryamye i obratnye zadachi mekhaniki kompozitov" [Methods for Solving Boundary Value Problems of Mechanics of Composite Plates and Shells: Teaching Guide on the Curse «Direct and Inverse Problems of Composite Mechanics»]. ICT SB RAS Publ., Novosibirsk, 2014, 131 p. [electronic resource]. ISBN 978-5-9905791-0-1.
[22] Vasilyev V.V. Mekhanika konstruktsiy iz kompozitsionnykh materialov [Mechanics of Composite Material Structures]. Moscow, Mashinostroenie Publ., 1988, 269 p.
[23] Blokhin A.M., Ibragimova A.S., Semisalov B.V. Matematicheskoe modelirovanie -Mathrmatical modeling, 2009, vol. 21, no. 4, pp. 15-34.
[24] Timoshenko S., Woinowsky-Krieger S. Theory of plates and shells. 2nd ed. N.Y.; Toronto; London: McGraw-Hill Book Company, Inc, 1959.
Golushko S. K. (b. 1958) graduated from Mathematical Faculty of Novosibirsk State University in 1980. Dr. Sci. (Phys.-Math.), head of the laboratory of Nonlinear Systems Analysis and Optimization at the Institute of Computational Technologies SB RAS, director of the Technological Design Institute of Computer Engineering SB RAS. He is an Expert of the Russian Foundation for Basic Research, Foundation for Assistance to Small
Innovative Enterprises in Science and Technology, Russian Science Foundation. The author and coauthor of more than 170 scientific papers, including one book, 5 Russian and 5 foreign patents. Research interests: composites, composite materials, composite structures, plates, shells, strength, fracture, mathematical modeling, computational mechanics, biomechanics, optimal design. e-mail: s.k.golushko@gmail.com.
Semisalov B.V. (b. 1987) graduated from Mathematical Faculty of Novosibirsk State University in 2010. Candidate of Physico-Mathematical Sciences, researcher at the Technological Design Institute of Computer Engineering SB RAS.Winner of the All-Russian competition of the "UMNIK" program (2013). The author and coauthor of more than 40 scientific papers, including 2 book, 3 chapters in monographs, 11 articles in journals from the list of Higher Attestation Commission, 8 publications indexed by WOS and Scopus. Research interests: mathematical modeling, numerical methods, nonlinear boundary-value problems, algorithms without saturation, best approximations, semiconductor devices, composites, hemodynamics. e-mail: vibis@ngs.ru.
