Аналитическое построение функций формы и матрицы жесткости кольцевого конечного элемента в методе перемещений Текст научной статьи по специальности «Математика»

УДК УДК 539.3

Ю. И. Фомиче в, Е. А. Коровайцева

АНАЛИТИЧЕСКОЕ ПОСТРОЕНИЕ ФУНКЦИЙ ФОРМЫ И МАТРИЦЫ ЖЕСТКОСТИ КОЛЬЦЕВОГО КОНЕЧНОГО ЭЛЕМЕНТА В МЕТОДЕ ПЕРЕМЕЩЕНИЙ

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

Постановка задачи. Функции формы конечного элемента (КЭ) и соответствующие им матрицы жесткости, могут быть построены как точно [1, 2], так и приближенно в рамках интерполяционного подхода [3]. В последнем случае возникает риск отсутствия учета смещения элемента как жесткого целого. Необходимость интерполяционного подхода, несмотря на возможные риски, вызвана тем, что область применимости точного построения функций формы ограничена возможностью построения аналитических форм общего решения систем дифференциальных уравнений. При использовании численных методов расчета последнее ограничение снимается и единственной областью рационального применения интерполяционного варианта метода конечных элементов (МКЭ) являются многомерные задачи: прямые методы решения задач динамики без использования аппроксимации, например метода Фурье, и многомерные задачи статики [4]. Но это вынуждает все расчеты проводить менее эффективно, чем при использовании функций формы [2], и только в дискретной форме.

Хотя с точки зрения трудоемкости расчетов определенным компромиссом является привлечение аппарата решетчатых функций [5], что позволяет оставить форму расчета дискретной, но преимущества аналитических методов расчета будут использованы не в полной мере.

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

Рис. 1. Расчетные схемы построения функций формы при неподвижном (а) и подвижном (б, в, г) узлах рамы

антиподом методу сил и предтечей МКЭ, удается применить и для криволинейных элементов рам. Далее на примере плоского кольцевого КЭ с углом раствора 90° приведена методика построения функций формы и матрицы жесткости такого элемента конструкции, в котором пренебрегается влиянием продольных и перерезывающих сил на состояние элемента в допустимых случаях.

Расчет матрицы податливости КЭ 6 и векторов Xг амплитуд функций формы в методе сил. Рассматриваемый КЭ имеет три степени свободы в узле. Для построения функций формы и матрицы жесткости элемента (МЖЭ) методами сопротивления материалов используем для каждого, например нижнего левого, узла рам (тип 1) с неподвижными узлами одну (рис. 1, а), а для рам (тип 2) с подвижными узлами три (рис. 1, б-г) расчетные схемы.

Их использование позволяет для расчета функций формы и МЖЭ даже для рам типа 2 ограничиться (аналогично [1]) тремя элементарными единичными эпюрами изгибающих моментов (рис. 2).

Эти единичные эпюры позволяют построить матрицу податливости КЭ 8 методами сопротивления материалов. С использованием мнемоники элемент 8^ матрицы податливости 8 определяется как

% = эИ]г) * э а сама матрица имеет вид

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

S =

п R

2 EIZ

симм

п -1 2-

3п - 2

4

jR2.

EIZ R3

jR2.

Elz 1 R3

2 EIZ

п R3

4 Elz

В работе [1] показан способ получения из матрицы податливости блоков матрицы жесткости при произвольном угле раствора кольцевого элемента, применимый к расчету конструкций по МКЭ. В строительной механике не всегда необходимо пользоваться МКЭ, например при анализе состояния элементов кроме МКЭ необходимо и построение функций формы, поэтому откажемся от упомянутого способа использования матрицы податливости. Для формирования функций формы и МЖЭ с использованием матрицы 8 введем сообразно числу степеней свободы в узле три вектора В и три матрицы П^ выборки функций формы:

1 0 0

0 ; В 2 = 1 ; Вз = 0

0 0 1

0 0 1 0 1 0

1 0 ; П2 = 0 0 ; Пз = 0 1

0 1 0 1 0 0

B i =

П =

Тогда векторы неединичных узловых амплитуд функций формы внутренних силовых факторов метода сил У определяются решением трех систем линейных алгебраических уравнений:

Щ8ПгУг = -Щ8В¿; г € [1, 3] , (1)

описывающих символьное обращение матрицы податливости относительно единичных узловых неизвестных Хц = 1.

Решение систем (1) позволяет сформировать три вектора узловых амплитуд функций формы в методе сил для одного из краев элемента:

~ 2п2 - 4п - 8 1 —8п + 24 1 " 1;

X =

3п2 — 8п — 4

X2T =

X3T =

п2 - 2п - 4 п2 — 8

• R; 1

R 3п2 - 8п - 4 —2п + 8

R

п2 8

R

- 4п + 12 D -2п + 8 л • R; —2—зт • R; 1

п2 8

п2 8

Рис. 3. Расчетные схемы метода сил

Аналогичный вид имеют и еще три вектора узловых амплитуд функций формы для другого края.

Каждый из упомянутых векторов определяет решение трех задач, формально не имеющих прямого отношения к методу перемещений, но строго соответствующих решению методом сил задач для расчетных схем, приведенных на рис. 3.

Они позволяют построить и сами функции формы внутренних силовых факторов для расчетных схем (см. рис. 3). Поскольку эти функции характерны только для метода сил, откажемся от их анализа и иллюстрации, ограничившись узловой амплитудой функций. Тогда, хотя для криволинейных элементов в общем случае не избежать применения ЭВМ, для плоских кольцевых элементов при расчете МЖЭ удается вообще обойтись без привлечения каких-либо дифференциальных уравнений и обращений матрицы податливости, а для построения функций формы достаточно привлечь простейшее уравнение.

Для реализации такой возможности предлагается использование метода сил с решенными ранее типовыми задачами как встроенной процедуры метода перемещений.

Расчет МЖЭ и функций формы моментов метода перемещений. Рассмотрим в качестве примера расчет единичных функций изгибающих моментов внутренних узлов конструкции и узловых значений продольных и поперечных сил. Такая постановка допустима для длинных стержней, влиянием продольных и поперечных сил на смещения которых можно пренебречь. В то же время она не препятствует составлению условий равновесия узлов рам как в виде суммы моментов (рама типа 1), так и в виде суммы сил (рама типа 2).

Для кольцевого элемента в виде четверти кольца расчетная схема внутреннего элемента конструкции в методе перемещений приведена на рис. 4.

Согласно методу перемещений для представленного элемента могут быть построены шесть единичных функций изгибающих моментов: две соответствуют единичным поворотам узлов А и В, четыре — единичным смещениям этих узлов в вертикальном и горизонтальном направлениях.

/

Рис. 4. Расчетная схема внутреннего элемента конструкции в методе перемещений

Рис. 5. Pасчетная схема внутреннего элемента конструкции в методе сил

Рассмотрим матрицы податливости КЭ 8 и векторов Х,ь амплитуд функций формы в методе сил для расчета одной из единичных функций изгибающих моментов в методе перемещений, например функции, соответствующей единичному повороту узла А. Отказавшись при этом на время от методологии метода перемещений и соответствующего задания величины фа = 1, составим расчетную схему метода сил, адекватную схеме, приведенной на рис.5. Тогда получаем расчетную схему, уже представленную на рис. 3, а. В этой схеме внешняя нагрузка задается узловой, согласно методу перемещений, но единичной, согласно методу сил. Используя для расчетной схемы (см. рис.3,а) соответствующий вектор 11 амплитуд формы метода сил, получаем функцию суммарных изгибающих моментов в элементе с расчетной схемой:

Отметим, что для функции (2) пока нет необходимости ее анализа и построения, так как она получена для расчетной схемы метода сил, а не перемещений. Используя эту функцию, найдем элемент матрицы жесткости кольцевого элемента, соответствующий узлу А. Для расчета, например, единичной функции моментов метода перемещений, кроме функции (2), используем функцию Для рассматривае-

мого примера эти две функции позволяют найти угол поворота узла А от единичного момента в этом же узле, т.е. податливость узла А при действии единичного момента в этом узле. Применяя мнемонику, запишем

3

MS= Е, i = 1.

(2)

Рис. 6. Реакции жесткого узла Рис. 7. Функция моментов при единичном при единичном повороте повороте жесткого узла

Полученное значение является не только коэффициентом матрицы жесткости, но и коэффициентом приведения вектора амплитуд функций формы в методе сил к методу перемещений. С его использованием получаем векторы У амплитуд функций форм внутренних силовых факторов в методе перемещений для узла А:

Y = КггХг; i е [1, 3] . (3)

При вычислениях с точностью до трех значащих цифр получаем для иллюстрируемого случая расчетную схему (рис. 6), эквивалентную схеме, приведенной на рис. 5.

Отметим, что эта схема используется, во-первых, в расчетах на прочность колец малой кривизны и расчетах функций формы перемещений КЭ для построения единичной эпюры изгибающих моментов метода перемещений и, во-вторых, для расчета значений коэффициентов Kj МЖЭ. Одна из шести возможных функций моментов, а именно

EI

Mz,i(ip) = —^(14,944 - 9,485 cos ф - 12,988 sin ф), (4) R

представленная на рис.7, имеет экстремум при ipext = arctg 1,369 и EIz

Mz (ipext) = —1,145——z, а первая строка МЖЭ имеет значения R

K\% = {5,459; —9,485; 12,988; 1,955; —12,988; 9,485} . (5)

Аналогично могут быть получены остальные пять функций моментов и строк МЖЭ. Оставшиеся 12 функций поперечных и продольных сил формируются последовательным символьным дифференцированием шести функций моментов, доставляющим и этим функциям строго аналитическую форму записи. В итоге, используя не числовую (4), а функциональную форму представления матрицы жесткости, получаем

Г Mzl(0) Mz2(0) Mz3(0) Mz4(0) Mz5(0) Mz6(0) 1

N1(0) N2(0) N3(0) N4(0) N5(0) N6(0)

В работе приведены 18 функций обобщенных усилий, значения которых при ф = 0 и ф = п/2 позволяют рассчитать числовые значения всех компонент МЖЭ (5), только одна — в соотношении (3) , а графики этих 18 функций не приведены. Они в целом типовые и легко получаются по описанному алгоритму.

Отметим, что для рам типа 1 нет необходимости расчета всех шести функций моментов в отличие от подхода МКЭ, а МЖЭ — всего второго порядка. Отметим также, что даже для рам этого типа для граничных элементов конструкции, кроме расчетной схемы, представленной на рис. 5, должны быть рассмотрены и иные, например, такие как на рис. 8, т.е. необходимо табулирование функций форм изгибающих моментов и перемещений типовых внутренних и граничных элементов конструкций, как это сделано для прямолинейных элементов рам.

Расчет функций формы перемещений метода перемещений.

Обычно в точных расчетах функций форм перемещений деформируемых элементов конструкций при единичных обобщенных смещениях их границ и МЖЭ этих элементов в локальной системе координат используется однородное дифференциальное уравнение, описывающее напряженно-деформированное состояние элементов, превращенных в изолированные, согласно методу перемещений. Такой подход удобен только для прямых элементов конструкций. Отметим нерациональность использования простейшего уравнения упругой линии в некоторых курсах строительной механики: для прямого конечного элемента уравнение четвертого порядка, описывающее его состояние, более

Рис. 8. Расчетные схемы для граничных элементов конструкции

универсально и проще в использовании при получении как функций формы, так и МЖЭ. Это обусловлено полиномиальной формой всех компонентов решения дифференциального уравнения и сравнительно низким порядком уравнения. Для кольцевых элементов ситуация обратная, так как полиномиально-тригонометрическая форма компонентов решения дифференциального уравнения шестого порядка приводит к необходимости многократного решения систем линейных алгебраических уравнений такого же порядка. Аналогичная ситуация возникает и при использовании функции Грина. Поэтому откажемся от использования в расчетах функций форм перемещений полного уравнения, описывающего состояние кольцевого элемента, и рассмотрим только дифференциальное уравнение, следующее из геометрических и физических соотношений:

d3V dV R2 . ,

dV + ¥ = - ERZM

Для этого уравнения функция Грина, позволяющая конкретизировать вид интеграла Дюамеля, имеет вид

Г(ф^) = 1 - cos(^ -

а общее решение уравнения (5) представляется в форме

v

R2

V(ф) = c1 + c2 sin ф + c3 cos ф - — Mz )Г (ф,£)^. (6)

0

Оно может быть использовано для расчета любой из трех функций форм перемещений конечного кольцевого элемента при соответствующих начальных условиях:

V (п/2) = 0; V'(n/2) = 0; V"(n/2) + V (п/2) = 0.

Например, для первой из функций формы, необходимой в расчетах рам обоих типов на жесткость, получаем c\ = R; c2 = 0; c3 = —R, а функцию формы получаем в виде

v

R2

V (ф) = R(1 — cos ф) — — Mz (£ )Г (ф,£ )d£.

EI

z J

0

Аналитический вид функции (6) достаточно громоздок, но для ее иллюстрации можно использовать любой аналитико-графический построитель, например Майсад-12. Отметим также, что предлагаемый подход не только позволяет получать функции формы перемещений кольцевого КЭ в аналитическом виде, но и оказывается вне конкуренции по сравнению с методикой, использующей стандартное диффе-

ренциальное уравнение состояния кругового кольца [6]

16У Л4V 12У . , .

+2 1^4+=11

так как позволяет как при расчете МЖЭ, так и при расчете всех функций формы состояния КЭ избежать многократного численного решения систем линейных алгебраических уравнений шестого порядка, не говоря уже о трудностях построения функций Грина.

Хотя авторы считают основным назначением настоящей работы изложение методики комплексного использования возможностей методов сил и перемещений для получения аналитических форм расчета МЖЭ на конкретном примере, расширенное представление методики, безусловно, включает в себя и возможности численных расчетов. Поэтому частично приведем некоторые решения, выполненные по изложенной методике. Для этого были выбраны задачи, представленные в главе 11 справочника [6]. В качестве контрольного решения специально было использовано общеизвестное решение задачи для второго семестра курса сопротивления материалов. Поэтому далее приведены основные результаты решения одной из таких задач на примере полукольца, закрепленного жестким защемлением на левом конце и неподвижным шарниром на правом конце и нагруженного вертикальной силой Р посредине полукольца.

На фоне эталонного решения, точного в классе общеизвестных допущений теории колец, были дополнительно проведены расчеты упомянутой задачи программным комплексом А^УБ. При расчете для примера было использовано 2, 10 и 40 КЭ с равным числом КЭ по каждой из четвертей полукольца. При использовании методики вопрос о числе КЭ вообще бессмысленен, так как матрица жесткости является абсолютно точной в классе аналитических функций, поэтому скорее уместно говорить в обычном стиле МКЭ об использовании метода суперэлементов. В отличие от комплекса А^УБ для контроля возможностей изложенной в работе методики ситуация была искусственно смоделирована уменьшением числа знаков представления аналитической формы МЖЭ. Такая ситуация, например, может быть и реальной, если для ЭВМ типа ЕС или 1СЬ не использовать в расчетах двойной точности представления чисел. Минимально достаточное число суперэлементов при решении представляемой задачи по методике, изложенной в настоящей работе, равно двум. Для сопоставления с результатами расчетов, выполненных с помощью комплекса А^УБ, общее число точек дискретизации расчета МЖЭ и их разделение по суперэлементам для аналитической ее формы, предлагаемой в работе, было взято таким же, как и при расчетах с использованием А^УБ.

На рис. 9, а показана копия функции формы изгибающего момента, получаемая с помощью комплекса А^УБ для элемента Р1ре-18.

О .038686 .077373 .116059 .154746

.019343 .05803 .096716 .135402 .174089

PRINT ELEMENT TABLE ITEMS PER ELEMENT

poSTl ELEMENT TABLE LISTING *****

STAT CURRENT CURRENT

ELEM MOMZI MOMZJ

1 0. 59611E- 01 -0. 27298E- 01

2 -0. , 27298E- 01 -0. 59675E- 01

3 -0. 59675E- 01 -0. 34349E- 01

4 -0. 34349E- 01 0. 46201E- 01

5 0. 46201E- 01 0. 17409

6 0. 17409 0. 27780E- 01

7 0. 27780E- 01 -0. 69387E- 01

8 -0. 69387E- ■01 ■ -û. 10790

9 -0. 10790 -0. 83992E- 01

10 -0. 83992E- 01 0. 16653E- 15

MINIMUM VALUES

ELEM 9 8

VALUE -0.10790 -0.10790

MAXIMUM VALUES

ELEM 6 5

VALUE 0.17409 0.17409 б

Рис. 9. Экранная копия функции формы (а) и результатов расчета (б) изгибающего момента для элемента Pipe-18 комплексом ANSYS при двух (а) и десяти (б) КЭ

Очевидно, что она не имеет никакого отношения к истине. Для сравнения отметим, что методика настоящей работы дает результаты той же точности, что и эталонное решение методами сопротивления материалов, т.е. с точностью машинного нуля в погрешности расчета. Таким образом, использование предлагаемой в работе методики превращает численный для комплекса ЛК8У8 расчет в иллюстрационный аналитический расчет. Поэтому все последующее вынужденное увеличение числа КЭ при расчетах с помощью комплекса ЛК8У8 не только приводит к риску появления вычислительных сложно контролируемых ошибок расчета, но просто и экономически нецелесообразно из-за резкого и неоправданного увеличения трудоемкости решения.

На рис. 10 показаны функции изгибающих моментов для десяти КЭ; здесь по оси абсцисс отложены номера КЭ. Отметим, что кривые практически совпадают и судить об их совместном существовании можно только при цветном представлении кривых. Поэтому авторы

1 1 3 4 5 б 7 8 9 10 Рис. 10. Огибающая функции изгибающих моментов для десяти КЭ

оставили за собой право представить анализ погрешностей и не интегральный, сглаженный расчет их, а расчет по одной из наиболее жестких норм погрешности — чебышевской.

На рис. 9, б приведена копия результатов с экрана дисплея, выдаваемых комплексом Л^У8 для десяти КЭ.

На рис. 11, а, также копирующем результаты расчетов с экрана дисплея, сопоставлены результаты расчетов моментов по методике, приведенной в настоящей работе (матрица А2), и с помощью комплекса Л^У8 (матрица А3) при точности сравнения в три знака. Действительно, различия почти нет.

В то же время анализ этих же результатов при семи знаках представления говорит о совершенно другом (рис.11, б). Здесь матрица

А2 =

1

1 0.061

2 -0.027

3 -0.059

4 -0.034

5 0.046

6 0.173

7 0.027

8 -0.07

9 -0.109

10 -0.084

A3 =

1

1 0 06

2 -0.027

3 -0.06

4 -0.034

5 0.046

6 0.174

7 0.028

8 -0.069

9 -0.108

10 -0.084

А1 =

1

1 0.0607804

2 -0.0265557

3 -0.0593741

4 -0.0344623

0.0457412

6 0.1733856

7 0.0269591

8 -0.0701881

9 -0.1085465

10 -0.0843613

А2 =

1

1 0 0607836

2 -0.026553

3 -0.059372

4 -0.0344607

5 0.0457422

6 0.173386

7 0.026959

8 -0.0701884

9 •0.1085469

10 -0.0843616

A3 =

1

' 0.059611

2 -0.027298

3 -0.059675

4 -0.034349

5 0.046201

6 0.17409

7 0.02778

8 -0.069387

'3 -0.1079

10 -0.083992

Рис. 11. Сравнение функций изгибающих моментов при описании их с точностью трех (а) и семи (б) десятичных знаков по предлагаемой методике и с помощью комплекса ANSYS

0.005 0.0025

2 - dMk

! —"Г -0.0025

1 - DMk

- -0.005

-0.0075 -0.01

0.005 0.0025 0

2 - ^k

, --0.0025

1 - DMk

- -0.005

-0.0075

-0.01---------

1 5 9 13 17 21 25 29 33 37 41

б

Рис. 12. Сглаженное развитие погрешностей расчета с использованием комплекса ANSYS (кривые 2) и предлагаемй методики (кривые 1) для десяти (а) и сорока (б) КЭ

А1 соответствует решению задачи методом сопротивления материалов (эталонному), а матрица А2 получена по предлагаемой методике согласно приведенным искусственным ограничениям.

На рис. 12, а показано (сглаженное для иллюстративности) развитие погрешностей согласно рис. 11, б при использовании чебышевской нормы. На этом и последующих рисунках кривая 1 соответствует результатам, полученным по методике, приведенной в данной работе, в условиях искусственного ограничения точности представления истинных результатов, а кривая 2 соответствует расчетам, выполненным с помощью комплекса ANSYS.

Казалось бы, погрешность расчета по методике настоящей работы равна нулю. Конечно, в условиях искусственного ограничения точности это не так. Графики (см. рис. 12, а) показывают лишь наглядно, что при десяти КЭ комплекс ANSYS на несколько порядков увеличивает погрешность расчета. Безусловно, в приведенном примере дана модельная ситуация, но именно она и возможна при проведении расчетов сложных конструкций, где нет эталонной базы расчетов, а малые погрешности расчета приводят подобно неправильно заданным начальным значениям к непредсказуемым эффектам и результатам. Поэтому на рис.13 показано развитие погрешностей (см. рис.12, а) только для расчета, выполненного по методике, предлагаемой в настоящей работе.

23456789 10

а

1 2 3 4 5 6 7 S 9 10

Рис. 13. Сглаженное развитие погрешностей расчета по предлагаемой методике

Наиболее интересными, по мнению авторов, являются сравнения результатов расчетов при сорока КЭ, представленные на рис. 12, б с теми же номерами кривых.

Они, с одной стороны, прекрасно подтверждают, что к мнению классиков, в данном случае методологии КЭ [7, с. 38], необходимо прислушиваться при анализе результатов расчетов: "Действительный минимум энергии никогда не может быть достигнут ни при каком числе разбиений, так как задание (выделено нами. - Авт.) функций формы ограничивает число степеней свободы системы. Чтобы гарантировать сходимость процесса к точному решению, необходимо удовлетворить некоторым простым требованиям". С другой стороны, выявленный эффект подтверждает и мнение специалистов, использующих более высокие технологии вычислений, чем конечно-элементные [7]; хотя формально при увеличении числа точек решения задач Коши разница с аналитическим решением должна стремиться к нулю, фактически этого не происходит и существует для каждой конкретной задачи свой оптимальный шаг дискретизации задачи.

Заключение. Задачи о построении функций формы метода перемещений и МЖЭ рассматривались во многих работах. Результаты настоящей работы показывают, что для плоского кольцевого КЭ совместное использование методов сил и перемещений позволяет построить аналитические формы расчета, допускающие расчеты на уровне машинного нуля. Полученные количественные оценки свидетельствуют о дополнительных возможностях использования аналитических методов расчета при разделении системы дифференциальных уравнений на составляющие меньшей размерности за счет отделения уравнений равновесия или их части от уравнений геометрических и физических при описании состояния КЭ.

Увеличение числа КЭ в распределении локальных погрешностей расчетов кардинально ничего не изменяет, так как не изменяются их причины: при использовании МКЭ в комплексе А^УБ (иные не тестировались) по-прежнему при расчете характеристик элемента используется интерполяция перемещений в терминах работы [3], а при

аналитическом расчете МЖЭ число знаков представления, специально заниженное и фиксированное, также не изменяется.

В итоге использование аналитического расчета МЖЭ позволяет, подобно известным работам В.Е. Поповича в области вычислительной математики, выполнять расчеты с погрешностью на уровне машинного нуля. Стандартная же интерполяционная версия МКЭ дает погрешность, в принципе даже не приближающуюся к машинному нулю и сохраняет неустранимые погрешности, зависящие от способа интерполяции функций обобщенных перемещений, их степени и иных, более важных причин, но не от числа конечных элементов.

СПИСОК ЛИТЕРАТУРЫ

¡.Образцов И. Ф., Савельев Л. М., Х а з а н о в Х. С. Метод конечных элементов в задачах строительной механики летательных аппаратов. - М.: Высшая школа, 1985. - 392 с.

2. Бакулин В. Н., Каледин В. О. Конечный элемент круговой арки с конечной сдвиговой жесткостью // Механика композитных материалов. - 1988. - № 8. - С. 915-919.

3. Морозов Е. М., Никишков Г. П. Метод конечного элемента в механике разрушений. - М.: Наука, 1980. - 254 с.

4. Расчеты машиностроительных конструкций методом конечных элементов: Справ. / Под ред. В.И. Мяченкова. - М.: Машиностроение, 1989. - 520 с.

5. Тюкин В. Н. Теория управления. Ч. 2. Особые линейные и нелинейные системы. - Вологда: ВоГТУ, 2000. - 128 с.

6. Прочность, устойчивость, колебания: Справ. Т. 1. - М.: Машиностроение, 1968. - 831 с.

7. Зенкевич О. Метод конечных элементов. - М.: Мир, 1975. - 541 с.

8. Современные численные методы решения обыкновенных дифференциальных уравнений / Ред. Д. Холл, Д. Уатт. - М.: Мир, 1979. - 312 с.

Статья поступила в редакцию 18.05.2007

Юрий Иванович Фомичев родился в 1948 г., окончил Московский физико-технический институт в 1970 г. Д-р техн. наук, профессор кафедры "Космические аппараты и ракеты-носители" МГТУ им. Н.Э. Баумана. Автор 68 научных работ в области прочности космических конструкций.

Yu.I. Fomichev(b. 1948) graduated from the Moscow Physical and Technical Institute in 1970. D. Sc. (Eng.), professor of "Spacecrafts and Launchers" department of the Bauman Moscow State Technical University. Author of 68 publications in the field of strength of space constructions.

Екатерина Анатольевна Коровайцева родилась в 1986 г. Студентка МГТУ им. Н.Э. Баумана. Автор 19 научных работ в области прочности конструкций.

Ye.A. Korovaitseva (b. 1986) — student of the Bauman Moscow State Technical University. Author of 19 publications in the field of strength of constructions.