Редуцирование расчетных моделей в задачах прочности

Том XXII

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 1991

М1

УДК 629.7.015.4 + 629.7.015.4 : 533.6.013.42

РЕДУЦИРОВАНИЕ РАСЧЕТНЫХ МОДЕЛЕЙ В ЗАДАЧАХ ПРОЧНОСТИ

Д. Д. Евсеев, Е. К. Липин, В. В. Чедрик

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

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

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

и комплексов программ решения многих сложных технических задач. Принимая во внимание различную структуру исходных данных континуальных, дискретноконтинуальных и дискретных моделей метода конечных элементов, возникает проблема формирования единых исходных данных о геометрических и жесткост-ных параметрах конструкции. Так как полные исходные данные о конструкции требуются для расчета ее НДС, то очевидно, что геометрическая и жесткостная модель метода конечных элементов может быть принята в качестве основной, из которой может быть сформирована соответствующая континуальная и дискретно-континуальная модели. Причем приведение матрицы жесткости диск-кретной модели МКЭ к соответствующим матрицам континуальной и дискретноконтинуальной моделей позволит на единой расчетной модели осуществить решение как задач прочности, так и аэроупругости [1]. Для решения данной проблемы предлагается метод редуцирования матрицы жесткости дискретной модели конструкции к матрицам жесткости континуальной и дискретно-континуальной моделей с помощью матриц перехода, полученных при формировании вектора узловых перемещений по заданным формам перемещений континуальной или дискретно-континуальной моделей.

1. ОСНОВНЫЕ СООТНОШЕНИЯ РАСЧЕТНЫХ МОДЕЛЕЙ

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

Г = {и, V, ш}Т = {«1, VI, , Щ, 1)1, и)С, . . . , и*, Уд,, 1Юн)Т,

где «'—индекс узловой точки; N—количество узлов дискретной модели; и, V,-, до,-—компоненты вектора перемещений г,- 1-й узловой точки в направлении осей х,' у, г глобальной системы координат. При известных соотношениях между силами и перемещениями для каждого конечного элемента с помощью МКЭ можно описать упругие свойства и исследовать поведение конструкции в целом. При этом через узловые перемещения по заданным формам внутри каждого конечного элемента определяется их НДС.

Вектор узловых перемещений г определяется из решения системы уравнения равновесия

Кг = я.

С математической точки зрения метод конечных элементов имеет тесную связь с методом Рэлея — Ритца, поэтому уравнения равновесия (1) могут быть получены как следствие условия стационарности полной потенциальной энергии дискретной модели конструкции

N N N

Э - и - т=4Е Е -Е ^ - я’г .

;=1 /=1 ;=|

где К— матрица жесткости дискретной модели; —--элементы матрицы К; Я—

вектор эквивалентных узловых сил.

Конструктивно-анизотропная пластина. Конструкция с помощью данной концепции представляется тонкой двухслойной пластиной, несущие слои которой могут быть подкреплены сосредоточенными элементами, создающими конструктивную: анизотропию. В рамках гипотезы «прямой нормали» перемещения любой точки пластины в направлениях осей х, у, г определяются через поперечный прогиб V (х, г) в следующем виде:

и (х, У, 2) = — у — , и (х, У, 2) = и (х, 2), Ш(х, у, 2) = — у ^~.

Для получения решения о деформированном состоянии пластины используется метод Ритца, согласно которому поперечный прогиб V (х, 2) представляется как линейная комбинация заданных форм перемещений:

П

V

к= I

где С = {Сь С2, ••• ,СП} т — вектор обобщенных координат; фк— заданные координатные функции, удовлетворяющие кинематическим граничным условиям. В качестве координатных функций могут выбираться степенные функции от безразмерных координат х и 2 вдоль направлений осей х, 2 [3]:

фк

где рк, дк — показатели степеней, значения которых определяются из условия

удовлетворения кинематических граничных условий. Вектор обобщенных координат С определяется

энергии конструктивно-анизотропной пластины:

п п п .

Э = и - Т = + 2 2 ЕктСкСт- 2 <?к С = + Сс GС — дтС

к=1 т=1 к=]

как решение системы уравнений равновесия обобщенных сил

GC=Q,

где G — обобщенная матрица жесткости пластины; я^. — элементы матрицы О; <3 — вектор обобщенных сил.

В результате определяются перемещения любой точки конструкции, представляемой конструктивно-анизотропной пластиной:

- п _ П

/ , дх V"1 ^ да дх

и (х, у, 2)=—у лг! Ск -%= -уж

к= 1 к= 1

п п

V (х, у, 2) = I Скфк = I СКХРк;/\ к= 1 к= 1

»

к= 1 к= 1

Тонкостенная балка. Конструкция представляется тонкостенной балкой, наделенной эффективными жесткостями на изгиб £7„зг и кручение G^кp [4]. Деформированное состояние балки без учета влияния поперечного сдвига представляется двумя независимыми функциями: то (2) — поперечный прогиб, 0 (2) — угол закручивания поперечных сечений относительно оси жесткости г, через которые определяются перемещения любой точки балки по следующим формулам:

v (х, 2) = И»о (2) + хв(г), и (х, у, 2) = — у — = — у8 "

V (Х, у, 2) = —у = —у. (Шо + хв').

п _(Р«

ЧРк*

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

р я

шо(г) = £ (г), 9(2) = £ А,1],(г).

=1 ;=1

Поперечный прогиб V (х, г) при этом может быть определен в обобщенной форме:

П

и(х, г) = £ Скфк(х, г),

*=1

где принимаются следующие обозначения: для к =1,2, .. . ,р С1 = АЬ ер| = г|5ь

С2 = А2, ер2 = '\'2; ... , Ср = Ар, фр = ,|5р; для к = р+1, р + 2. п Ср+| = ВЬ

Фр+1 = х11ь Ср+2=В2, ерр+2 = хч2; ..., С„ = Вя, ф„ = х1]я; А, В — вектора обобщенных координат; % (г), г^ (г) — заданные координатные функции, удовлетворяющие кинематическим граничным условиям в задачах изгиба и кручения.

В качестве координатных функций , 1] могут выбираться степенные функции от безразмерной координаты г вдоль оси жесткости. Тогда поперечный прогиб v в обобщенной форме через соответствующие координатные функции может быть представлен выражением (3), где для к = 1, 2, . .., р показатели степени рк равны О, а для к = р + 1, р + 2 , ..., п показатели степени рк принимаются равными 1; 9К — показатель степени, значения которого определяются кинематическими граничными условиями для функций прогиба шо и угла закручивания 0. Вектор обобщенных координат С определяется из условий стационарности полной потенциальной энергии тонкостенной балки:

L L

Э = и — Т = + (£/изА2 + С/кр9'2) С1г — (£о^о + ткр0) ^г =

о о

п п п

= + I £ £ктСкСт — £ РкСк = + СТСС-РТС

к=1 т=1 к=1

как решение системы уравнений равновесия (4) , где О —обобщенная матрица

жесткости балки; £о(г)— погонная поперечная нагрузка; ткр (г) — распределенный крутящий момент. В результате определяются перемещения любой точки конструкции, моделируемой тонкостенной балкой:

л Я Р я л Я

и (х, у, г) = — у £ С<г “, V (х, у, г)= £ Скг “ + £ Скхг к,

к=р+1 * к=1 к=р+1

к>(х, у, 2) = — СК(<7К — 1)5(?к }+ ^ Ск^(^к—1)2<'к V

\А=1 *=Р+1 /

2. РЕДУЦИРОВАНИЕ РАСЧЕТНЫХ МОДЕЛЕЙ НА ЗАДАННЫХ ФОРМАХ ПЕРЕМЕЩЕНИЙ

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

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

Соответствие между дискретной и континуальной моделями конструкции наибо

точках. Из отмеченного различия в аппроксимации .

что в качестве общих точек наиболее удобно выбрать узлы конечно-элементной расчетной схемы. Таким образом, получим связь между ее узловыми перемещениями г и обобщенными координатами С при выборе координат континуальной модели равными координатам узлов х = х,-. у = у|, 2 = 2,-:

где П — матрица перехода от перемещений континуальной модели к узловым перемещениям дискретной модели:

Замену дискретной модели конструкции континуальной расчетной моделью при установленной связи узловых перемещений г и обобщенных координат С в форме (5) можно произвести из условий равенства их потенциальных энергий деформации и и внешних сил Т. Данные"условия позволят установить связи между матрицами жесткости и обобщенными силами дискретной и континуальной расчетными моделями. Подставим соотношения (5) в выражения для и и Т дискретной модели и получим равенства

V (х„ у„ г,) = £ Скх‘г-' = V,,

или в матричной форме

г={и, V, до}т= ПС,

/ дг »(*|—1 )\ / дг _р2 _(^2 1X

{-У^Ях* ^zx ) {-у11ГЯ2Х1г1 )

)

и = тгТКг = т СТПТКПС, Г = ЯТГ = ЯТПС.

Отсюда' из условий равенства потенциальных энергий дискретной' и  континуальной моделей . ?

уГтКг = 1стСС, Ятг = фтС

следуют соотношения между матрицами жесткости и обобшенными силами ввиде

0° = Пт/Ш, д°=ПтУ?.

В результате полученных соотношений (6) дискретная расчетная модель конструкции' приводится к континуальной'расчетной схеме конструкции, деформированное состояние которой описывается линейной комбинацией заданных форм перемешений на всей области конструкции. Неизвестные обобшенные координаты С определяются из решения уравнений равновесия обобшенных сил

О°С°=0°, С0=С0_1<30,

где ф°, 0 — вектор обобшенных сил и обобшенная матрица жесткости расчетной модели, полученной в результате операции редуцирования (6). При этом следует отметить, что для дискретной модели может быть сформирована сингулярная матрица жесткостИ К, при получении которой в уравнениях равновесия не обязательно исключать степени свободы, обеспечиваюшие отсутствие смешения конструкции как твердого тела, и учитывать кинематические граничные условия. Это становится возможным, так как заданные формы перемешений континуальной модели выбираются из условий удовлетворения заданным кинематическим граничным условиям. В этом случае операция редуцирования в форме (6) обеспечивает единственность решения уравнений равновесия обобшенных сил для континуальной расчетной модели, так как при формировании обобщенной матрицы С° редуцированием на заданных формах исключаются перемешения конструкции как твердого тела и учитываются кинематические граничные условия. Получив решение для континуальной модели конструкции в форме (7) и используя условия связи обобшенных координат С с вектором узловых перемещений г дискретной модели, можно по узловым перемешениям

г° = ПС°

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

о° = ОВг°,

где D — матрица упругости, В— матрица дифференциальных операторов, D • В —матрица напряжений.

3. ФОРМИРОВАНИЕ ДИСКРЕТНО-КОНТИНУАЛЬНОЙ МОДЕЛИ КОНСТРУКЦИИ НА ОСНОВЕ ДИСКРЕТНЫХ МОДЕЛЕЙ МКЭ .

Сложные силовые конструкции планера ЛА, а также нерегулярные конструкции его агрегатов в. задачах аэроупругости заменяются набором й упругих подконструкций, представляемых в расчетной схеме континуальными моделями, связанными в -общих узловых точках [5]. В этих точках упругость связей между подконструкциями моделируется набором пружин с соответствующими жесткостями. Подконструкции в'этой модели могут представляться как тонкостенной балкой, так и конструктивно-анизотропной пластиной. Деформированное состояние каждой

описывается через поперечный прогиб  V,- (х, ?) в форме (2)."

V, (х, г) = £ Ск/фк/(Я г) ,

где j — номер подконструкции. Связь между /-й подконструкцией и ^-ми подконструкциями, имеющими с ней смежные границы, устанавливается по оотносительным перемещениям общих узловых точек, которые согласно (5) определяются в следующем виде:

Дг/Ч, = Г; — Г^ = П/С; — П^с,.

Здесь П, и П,., — матрицы перехода в представлении перемещений и, V, Ф /-й и ">-ых подконструкций по заданным формам степенных функций (3) к узловым перемещениям на общих границах рассматриваемых подконструкций. Вектора обобщенных координат С определяются из условий стационарности полной потенциальной энергии дискретно-континуальной модели, включающей энергии й упругих подконструкций О, и энергию деформации пружин связи Ди/Ч

й й й

э = и Т = I (и; + ди*) -I П = -4-1 (C]GiCi + Ar],.,K,.,Mi1j» -/=1 /= /=

й й

-I о?С = 4-1 (с;/,+с;п; к^п,с, - с;п;к=,п,.,с, +

/=1 /=1

й

+ с;п;к;^п,.,с,., - с;п;к", п;.с;.)-£ о?с,

/='

где С/, 0/— обобщенная матрица жесткости и вектор обобщенных сил /-й подконструкции; —диагональные матрицы жесткости пружин связей, связывающих в узлах общей границы ' относительными перемещениями Д/= = { Дм;Ч, Ди;,." , Ддо;Ч^ т, являющимися согласно (9) функциями обобщенных координат, /-ю подконструкцию с ф-ми подконструкциями. При соответствующем задании в матрице диагональных коэффициентов жесткости пружин связей в расчетных моделях как дискретно-континуальной, так и соответствующей дискретной модели метода конечных элементов могут быть . реализованы любые связи (шарнир, упругое или жесткое закрепление). В энергию деформации для /-й подконструкции О, включаются энергии деформации пружин связи Ди,,." , установленных на общих границах с ^-ми подконструкциями, которые затем уже не суммируются в энергии деформации этих смежных подконструкций.

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

где Дг,'Ф = г, — г'Ф — вектор относительных перемещений общих узловых точек на смежных границах /-й и 'Р-ых подконструкций; К/, Я/— матрица жесткости и вектор эквивалентных узловых сил для /-й подконструкции. Вектора узловых перемещений Г; для подконструкций /=1,2, ... , (й— 1), й, образующих в рамках МКЭ расчетную схему конструкции, определяются из решения системы уравнений равновесия, полученной как условия стационарности полной.потенциальной энергии (11).

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

Г/=П/С/, Дг/Ч|,— П/С/ — Пу,Сф,

где П/,  П'Ф — матрицы перехода от перемещений континуальных моделей у-и и 'Р-ми подконструкциями к узловым перемещениям их дискретных моделей, может быть осуществлен из условий равенства потенциальных энергий деформаций и и внешних сил Т. Подставив соотношения (12) в выражение (11), получим связь между матрицами жесткости и обобщенными силами дискретной и дискретно-континуальной расчетными моделями как для подконструкций, так и всей конструкции:

й а

и = +2 (г? / + Аг^А^) = + 2 [ с; "? К/Н/С/ +

, +(С/СП,С- С£Г1$) /С/С,(П,С/- II,С,)], т = 2 Rfг,-= 2 Rfuicj.

• 1  /=! 1=1

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

с; = п; к,п/, о>°

= п;к ^П/,

,  = пфк^п,, се;

В результате дискретная модель конструкции приводится к дискретно-континуальной расчетной схеме, деформированное состояние которой описывается набором сумм роздействия заданных форм перемещений на каждой подконструкции. Неизвестные вектора обобщенных координат С] являются решением уравнений . равновесия обобщенных сил Q,( полученных из условий стационарности полной потенциальной энергии редуцированной дискретной модели конструкции. Матрицы жесткости К, в дискретной модели могут быть сформированы как для свободных подконструкций, а отсутствие их смещений как твердого тела и кинематические граничные условия учитываются соответствующим выбором форм перемещений в подконструкциях или при формировании относительных перемещений упругих связей. Используя условия связи обобщенных координат С'( с вектором узловых перемещений г°

г'(=П/С'(,

по известным соотношениям метода конечных элементов (8) определяются напряжения в элементах, составляющих подконструкции.

= П/Я/;

= ПТК); П * ; - ПТкС; п г

4. ПРИМЕРЫ РЕДУЦИРОВАНИЯ РАСЧЕТНЫХ МОДЕЛЕЙ

КОНСТРУКЦИИ

Многостеночная регулярная конструкция крыла (рис. 1). В дискретной модели МКЭ расчетная схема крыла содержала 120 узлов с 324-мя неизвестными перемещениями Иі, У;, W;. Обшивка М0делир0валась изопарамет-рическими конечными элементами, работающими в сложном напряженном с0стоянии, а продольные и поперечные стенки представлялись двумя типами элементов — изопараметрическими КЭ и элементами, работающими тольк-о на сдвиг. Континуальная модель консоли крыла представлялась в виде пластш^ поперечный прогиб которой описывался следующим рядом координатных функций:

П

. и(х, г) = 2 С.х к г *,

где показатели степени рк и дк с учетом жесткого закрепления корневого сечения принимались равными

Рк = О, 1, 2, 3, О, 1, 2, 3, О, 1, 2, 3, О, 1, 2; дк = 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5.

НДС консоли крыла в дискретной континуальной и редуцированной моделях определяются от действия сосредоточенной силы Я = 104 Н в концевом сечении. На рис. 2—4 представлены результаты расчета поперечных прогибов

и, компонент напряжений ах,, ау„, тх,у, и аэка в обшивке для рассмотренных моделей консоли крыла.

Напряжения приведены для двух сечений крыла 2= 0,05 и 0,95 (соответственно корень и конец крыла). В корневом сечении все компоненты напряжений и эквивалентные напряжения практически совпадают в МКЭ и модели пластины, полученной редуцированием конечно-элементной матрицы жесткости. Стоит также отметить, что даже в концевом сечении (2= 0,95), где напряжения на порядок ниже различия в напряжениях не существенны.

Из сравнения результатов следует, что редуцирование дискретной модели МКЭ к континуальной модели при существенном уменьшении разрешающей системы уравнений (324Х Х324 — МКЭ и 15Х 15 — пластина) не снизило точности определения НДС регулярной конструкции крыла.

Многостеночная нерегулярная конструкция крыла (рис. 5). Дискретная модель МКЭ данного крыла,' содержащего протяженный вырез, составлялась из трех регулярных: под-конструкций. Полная расчетная схема крыла со-

6 пластина

} сОНогоНь/е стенки о редуцированиеJ

—НКЗ

X

Рис. 2

г=0,05

6,і,

кг/пп1

20

,О

кт/им11

20

10

И

Рис. 3

Рис. 4

1 2 3 ¥ 5 б 7 8 9 10

Рис. 5

держала 120 узлов с 324 неизвестными перемещениями. Обшивка и стенки моделировались конечными элементами, используемыми. в расчегной схеме регулярной конструкции_крыла. Дискретно-континуальная модель сосгавлялась гакже и3 треХ подконструкций в виде пластин, где поперечные Прогибы описывались одинаковыми рядами координатных функций

V; (*- С*,*ГЧ '' І = 1'2'3;

*=і

Рк, = О, 1, О, 1, О, О, О, 1, 1, 2, 2, 3, 3; ?„, = О, О, 1, 1, 2, 3, 4, 2, 3, 2, О, О, 1.

, пластина

• ПКЭ

— МКЭ

— МКЗ

х редуцирование

Рис. 6

При данном наборе координатных функций граничные условия удовлетворяются соответствующим выбором жесткостей пружин связей,  а именно, в корневом сечении они задаются в узлах связи по направлениям перемещений и, и, и бесконечно большими. Разрешающая система уравнений в МКЭ имела размерность (324Х324), а в дискретно-континуальной модели составляла (ЗХ 13) X (ЗХ 13) = 39X39. Результаты расчетов дискретной, дискретноконтинуальной и редуцированной моделей в виде поперечных прогибов для нескольких поточных сечений представлены на рис. 6. Сравнение и анализ полученных результатов показывает, что разработанная процедура редуцирования и на нерегулярных конструкциях обеспечивает достоверные результаты по моделированию их упругих свойств.

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

ЛИТЕРАТУРА

1. И л ь ич е в В. Д. Матричиые методы сиитеза дииамических и упругих характеристик лииейной неконсервативной коиструкции. — М.: Ученые записки ЦАГИ, 1975, т. 6, № 2.

2. 3 е и к е в и ч О. Метод коиечных элементов • в технике. — М.: Мир, 1975.

3. Б у н ь к о в В. Г. Особенности свободной схемы летательного аппарата при решении задач аэроупругости. —М.: Труды ЦАГИ, 1969, вып. 1166.

4. Л и п и н Е. К., Т е и я е в а В. Е. Определение жесткости на кручение силового кессона многонервюрной конструкции несущей поверхности летательного аппарата. — М.: Ученые запискн ЦАГИ, 1988, т. 19, №, 4.

5. Е в с е е в Д. Д., Рыб а к о в А. А. Алгорнтм расчета матриц податливости конструкций летательных аппаратов методом подконструкций применительно к задачам аэроупругости. — М.: Ученые записки ЦАГИ, 1981, т. 12, № 5.

Рукопись поступила 18/Х 1989 г.