Применение метода дискретных операторов к расчету анизотропных пластин и оболочек из композиционных материалов Текст научной статьи по специальности «Физика»

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

________ 1 Q 87 М 2

УДК 539.3 : 678

ПРИМЕНЕНИЕ МЕТОДА ДИСКРЕТНЫХ ОПЕРАТОРОВ К РАСЧЕТУ АНИЗОТРОПНЫХ ПЛАСТИН И ОБОЛОЧЕК ИЗ композиционных МАТЕРИАЛОВ

А. А. Постников

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

Технология создания композиционных материалов и конструкций из них тесно связана с решением задач механики деформируемого твердого тела. Линейная теория упругости анизотропного тела [1]— один из ее разделов, имеющий широкий спектр приложений. На основе соотношений теории упругости с помощью гипотез и допущений построены различные по сложности модели [1—3, 13], позволяющие анализировать специфические особенности поведения композитов в конструкциях. Поскольку точные решения возникающих при этом задач удается получить лишь в частных случаях анизотропии материалов и геометрии исследуемых объектов, особую актуальность приобретает разработка приближенных методов расчета. На практике наибольшее распространение получили методы конечных разностей [4] и конечных элементов [5].

Метод конечных разностей (МКР) зарекомендовал себя при расчете пластин и пологих оболочек более экономичным, чем вариационные методы и метод коллокаций [6]. Однако за видимой простотой процедуры МКР стоит сложность практической реализации метода на нерегулярных сетках для расчета неканонических по форме в плане пластин и оболочек. Зная физический смысл исходных дифференциальных уравнений, разностные уравнения можно получить, используя ин-тегро-интерполяционный метод или метод балансов [4]. Если исходить из математического представления задачи, то для получения конечно-разностного аналога удобно использовать метод контурной аппроксимации [7, 8]. Процедура этого метода использует формулы Грина для установления связи между интегралами по области и интегралами по границе области. Данные подходы применялись при построении конеч-но-разностных схем на регулярных и нерегулярных сетках для дифференциальных уравнений в частных производных второго и четвертого

порядков, соответствующих задачам плоской теории упругости, изгиба и колебаний пластин [9, 10].

Метод конечных элементов (МКЭ) построен на вариационной основе и потому свободен от характерных для МКР недостатков. Однако расчет МКЭ пластин и оболочек из композиционных материалов связан с разработкой сложных конечных элементов [11, 13], что приводит в итоге к завышению порядка разрешающей системы алгебраических уравнений. В этой связи понятен интерес к методам, сочетающим простоту и экономичность МКР с универсальностью МКЭ. Известно, что основные соотношения МКЭ можно получить непосредственно из дифференциальных уравнений, минуя этап отыскания функционального эквивалента исходной краевой задачи. Преимущества такого подхода обсуждались в работе [5], а примеры реализации для задач плоской теории упругости и изгиба пластин содержатся в работах [14—18]. В настоящей статье для построения основных соотношений МКЭ непосредственно из дифференциальных уравнений используется метод контурной аппроксимации. В качестве аппроксимирующих функций выбираются функции формы конечных элементов. По аналогии с МКР вводится понятие операторной области, которая строится вокруг каждого узла конечно-элементной сетки. Дискретные схемы для дифференциальных операторов получаются с помощью теоремы о среднем и формулы Грина. Представленную процедуру можно условно назвать методом дискретных операторов (МДО) [15]. Имея много общего с МКР и МКЭ, алгоритм МДО выгодно отличается от традиционной формы МКР простотой постановки и реализации граничной задачи, а по сравнению с МКЭ имеет пониженный порядок разрешающей системы уравнений.

1. Постановка задачи. В работах последних лет установлены пределы применимости различных вариантов уравнений слоистых анизотропных пластин и оболочек для расчета элементов конструкций из композиционных материалов. Не обсуждая этот вопрос в данной статье, рассмотрим в качестве исходных уравнения весьма пологих анизотропных оболочек в смешанной форме [3], основанные на гипотезах Кирхгофа — Лява для всего пакета слоев:

1Х (А*)® + у*? = 0; |

М^ы)? —= 0, } ^

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

жесткости и податливости пакета по известным формулам [3]; q — ин-

тенсивность распределенной нагрузки; Ьи Уд — линейные дифференциальные операторы следующего вида:

I. Ф») = Ли тт + 4018 71Т + 2 (£>12 + 20бв) -Т71Т +

дх1 дх3 ду дх2 ду"

, .д4 . г-. д4

дх ду3 ду4 ’

‘■ЛЛ„)-А,г-~-2А,,-^- + (гАп + } (2)

_М, * +Д, * ■

’ дх ду3 ду*

_ 1 г____1 с*2

Ул ~ Я22 "т“ 7?,, <ЭуЗ ’

где х, у и /?1ь /?22—координаты и радиусы кривизны срединной поверхности оболочки соответственно.

Выбор уравнений (1) обусловлен их простой структурой, удобной для наглядного представления процедуры построения дискретного аналога исходной дифференциальной задачи. Если в уравнениях (1) положить V* =0. то система распадется на два независимых уравнения, относящихся к задачам изгиба пластин и плоской теории упругости. Поставим в соответствие срединной поверхности оболочки некоторую двумерную, регулярную область Й, ограниченную контуром I В качестве искомых выбираем функции до, ф, относительно которых формулируются граничные условия. Задача сводится к отысканию неизвестных функций до, ф, удовлетворяющих в области й уравнениям равновесия (1) и граничным условиям на контуре Ь.

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

Ф + + (3)

где Ф, Л'о, А^г — векторы следующего вида: Ф = {да, <р}т; Л^0 = = {<*<>, «з}т; = <*4)т; Л/2 = {а2, а5р-

Если обозначить значения вектора функций Ф в узловых точках элемента через Фь Ф2, Фз, то векторные коэффициенты N0, N1, Ы2, входящие в выражение (3), представляются в следующем виде:

Л'о = [(Х1У2 — У1 х2) Ф0 + (х, у0 —у2 х0) + (х0Уг — Уо ) Ф2];

А'1 = -£-=-^-(^12Ф0+^20Ф1+Л1Ф2); (4)

N2 = (-«12 ®0 + -*20 Ф1 + -*01 Фг).

где х1) = х1 — ху, у,7 =у1 —у/, Д1—площадь треугольного элемента; 2Д] = х10 у20 х20 Ум.

Далее, по аналогии с МКР [4], вводится понятие операторной области 50, которая состоит из частей элементов, окружающих узел А0, и ограничена контуром /, соединяющим центры элементов и середины их общих сторон

й г=1

Узлы конечно-элементной сетки Аг (XI, уi) классифицируются на внутренние и граничные (рис. 1).

6

Рис. 1

Дискретные схемы для дифференциальных операторов второго порядка строятся путем аппроксимации последних в операторной области 50 с помощью теоремы о среднем и использования затем формулы

После подстановки выражения для из (4) в правую часть ин-

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

В отличие от выражений (6) дискретные схемы для дифференциальных операторов во внутреннем узле не содержат интегральных слагаемых. Идексация в формулах (7) также зависит от типа узла. Так, при обходе внутреннего узла Л0 против часовой стрелки г—1 = гс, когда г=1, и 1+1 = 1, когда / = п. В то же время для граничного узла г—1 = 0, когда /=1, и /+1=0, когда г = п.

Для того, чтобы получить дискретные выражения для дифферен-

М ф

Грина. Для-------, например, получим

дх2

(5)

(6)

где

а, о (х) =

450

(х*=ьу);

(7)

д2Ф

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

= ЧГ. Тогда, например, для во внутреннем узле (рис. 2) можно

получить следующую схему:

где

тп

Аналогичным путем строятся дискретные выражения и для других дифференциальных операторов четвертого порядка. Приведем их окончательный вид для граничного узла Ло (см. рис. 2):

п Г т

[^]о = 2[а*о(^2а;/И(фу-ф|) -

— |^2 а1 о (*>] 2! а10 (*) (Ф/ - фо) +

+ е*1МСЮ-,1у+°*,1х) Г ^-йу-

^ J дх ' ^ 5А2 ] дх *

2 а1 О (*) а а

50

[^]. = ,| ["' “МФ<> —^2 °1 о (*)| 2 а1 о (у) (ф/ - фо) — — £*_!_(£) ^йх +

$н\ J ду ^2 J

+

2а‘о(*) а

<-1 С дФ

д3 Ф

Г—Лх -)—— Г ■

50 J ду .] дхдуз

йу, (х+=*у);

г а1 ф 1 _ ^

[ дх3ду ]0 ~ £

а/оИ2ал(*У) (фу-ф«)

У=1

2а*о(*)

г=1

Ё<*ю(*У) (Ф^~Ф0) +

12=1

5*1 ду 6’й2 Л ду

дз Ф

2 аг о (х) в

50 л ду Яо .1 д*2ду

6 ь

йу, (х ч=ь у).

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

3. Построение дискретного аналога исходных уравнений. Используем полученные дискретные схемы (8) для аппроксимации дифференциальных операторов (2). Интегральные члены, входящие в выражения (8), позволяют учесть условия закрепления контура оболочки. Для ортотропной слоистой оболочки, положив в уравнениях (2) .016=1)26=0 и Л16=Л26=0, получим

-^вн +

ац I (*) 5*

*0+

ак 1 (У)

1{ки~д^йу~кп^лх)

50

С (*« 4“ аУ — -Т7 аХ) +

дх

ду

ь

д3 ® I ьг Ф А12

дх3

дх ду1

•)«*у -

(9)

где Лвн — дискретный аналог уравнений во внутреннем узле А0\

п т

Л.„ = «о Фо + 2 Р< Ф< + 2 Юф;! “о. Р/. V ~ матрицы:

/=1

/=1

п п

ао!=Ца10(х)-В1 + 2а/о(у)-Я2 +

1-1 1=1

( Г т п 1 Г т п

Р/ = — |51 2 а)1 (*) + 2 а1 о (*) + В2 2 аП (У) +2 а1 О (У)

I Ь«1 1=1 ] Ь/=1

п п

Ту = В12 ап (х) а«о И + 52 2 (У) аю(у)>

г-1 /ш1

/I /I

#1 = *112 °/ о (*) + £ а1 о Су):

1=1

п

1=1

п

В2 — К12 2 а1 о (X) ~Ь *22 2 °г 0 (у)’

Кп =

оп о

О А 22 „

1=1 (=1

012 + 2066

*12 = *22 ==

L

О22 О

о _

О — (2Л13 + Ава)

/?

О

- 2 а, о (■*) + 2 о (У)

*М1 /=1 *'22 /—1

1 П 1 П

р“ 2 а/ о С*) + -р-2 а * о Су) 1^1 ^22 я

Р'0— вектор, /?0 = {А,о. °}т; Яй = -^-\)Я^хйу.

В уравнении (9) учтены условия, определяющие способы закрепления и нагружения боковых поверхностей прямоугольной в плане оболочки. Для учета граничных условий, которым должны удовлетворять функции прогибов ш и напряжений ф на криволинейном контуре оболочки, выберем за систему координат направление внешней нормали п и касательной t в граничном узле Ад. С учетом преобразований

дФ

дх

йу —— с!х = Г-^5- соз(л, х) +сое (га, у)

[ дх

дФ

ду

дФ

ду

<Эф Л, = —-Ш] дп

дх3

дхду2

дх3 ду

ду3

уравнение (9) принимает следующий вид:

Лвн + Г йЬ + ^ Г — <Н---------------------------^-------Г — Л = Г,, + Я0. (10)'

5Л1 ] ди ^ Я** ,) дп Б0 } дп о . о V /

а <1 Ь

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

Авн + Лог + Лг == Ро ) (11)

где

-^ог Л

и

п п

Е о м

І = 1 1=1

А=1 I 1=1

Лг = 2 0 2 [*11 аі о (я) + *12 «/ о (0] (ф< - Ф0); <30= р0+ Ро'>

' Я|о(л) '

а1о (і)

сое2 (га, л:)зіп2(га, л:)8іп(га, х)со8(п, х)

сіп2/и оло2 / и *уЛ сіп / и ■ъЛ пг\ с ( п

*10 (х)

а*о(У)

аіо(х)+аго(У)-^

Влияние граничных условий учитывается с помощью дискретных операторов во внутренних околограничных узлах Лог и граничных узлах Лг. Если, например, рассматриваемый узел является внутренним, то в уравнении (11) следует положить Лог=Лг = 0 и /5о=0. Для внутренних околограничных узлов следует положить равными нулю Лг и Р0. Дискретная схема для Лог сохраняется целиком только в случае свободного края и не учитывается в случае заделки.

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

Естественным оказывается обобщение методики на решение задач о колебаниях пластин и оболочек. При этом массовая плотность материала, как и давление, считается постоянной в каждом конечном элементе.

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

6 «Ученые записки» 2 81

Рис. 3

чи, характеризующее объем требуемых вычислений. Для МДО оно равно количеству узлов Л/, а для использованных высокоточных треугольных конечных элементов — 6УУ2+12А^ + 6. Следует отметить, что решение на основе прямоугольного изгибного конечного элемента пластины практически совпадает с полученным МДО, однако при равенстве количества узлов общее число степеней свободы для МКЗ в три раза больше. Эта разница становится наиболее ощутимой при решении нестационарных и нелинейных задач. Еще большего выигрыша можно ожидать при расчете оболочек.

В качестве примера применения МДО к расчету оболочек рассмотрим задачу об изгибе пологой сферической панели при действии равномерного внутреннего давления <7. Панель переменной толщины, характеризуемая параметрами /?=1,17 м, 0 = 27°, изготовлена из армированных волокон со свойствами £1 = 0,16-106 МПа, £2 = 0,5-103 МПа, 612 = 0,45 • 103 МПа, Vl2 = 0,28 путем звездообразного армирования по

схеме фг=ф(г—1), ?= —, /= 1, 2,..., 8. Коэффициенты упругости в урав-

8

нениях (1) вычисляются по известным формулам теории армирования [3]. На рис. 4 показано изменение меридиональных напряжений ое в сечении по образующей в зависимости от координаты х/а. Максимальных значений напряжения достигают вблизи заделки на внутренней поверхности панели (кривая 1), при этом параметр толщины панели Н!к изменяется линейно от #//1 = 50 в опорном сечении х/а= 1 до 11111= 100 в сечении х/а = 0,9 и далее остается постоянным. Точность

численного решения данной задачи можно частично проконтролировать с помощью известных аналитических решений [6, 19].

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

ЛИТЕРАТУРА

1. Лехницкий С. Г. Теория упругости анизотропного тела.—

М.: Наука, 1977.

2. Мал м ей стер А. К., Та муж В. П., Тетере Г. С. Сопротивление полимерных и композиционных материалов. — Рига: Зинатне, 1980.

3. Амбарцумян С. А. Общая теория анизотропных оболочек. —

М.: Наука, 1974.

4. Самарский А. А., Андреев В. Б. Разностные методы для эллиптических уравнений. — М.: Наука, 1976.

5. Зенкевич О. Метод конечных элементов в технике. — М.: Мир,

1975.

6. Справочник по теории упругости. /Под ред. Варвака П. М. и Рябова А. Ф. — Киев: Стройиздат, 1971.

7. М а с (Neal R. Н. An asymmetrical finite difference network. —

Quarterly of Applied Mathematics, 1951, vol. 11, N 3.

8. Бухман В. E. Метод контурной аппроксимации дифференциальных уравнений в частных производных системой конечно-разностных уравнений. •—В кн.: Об эффективности аналоговых методов решения краевых задач.—М.: Наука, 1969.

9. Приказчиков В. Г. Интегро-интерполяционный метод построения разностных уравнений в задаче колебаний пластины. — Ученые записки ЦАГИ, 1973, т. 4, № 4.

10. Lis zk а Т., Orkisz J. The finite difference method at arbitrary irregular grids and its application in applied mechanics. — Computers

and Structures, 1980, vol. 11, N 1—2.

11. Reddy J. iN., Chao W. C. A comparison of closed-form and

finite-element solutions of thick laminated anisotropic rectangular plates.— Nuclear Engineering and Design, 1981, vol. 64, N 2.

12. Пелех Б. JI., Марчук М. В. Метод конечных элементов при

решении краевых задач для анизотропных пластин из композитных материалов. — Механика композиционных материалов, 1983, № 1.

13. Алфутов Н. А., Зиновьев П. А., Попов Б. Г. Расчет многослойных пластин и оболочек из композиционных материалов. — М.: Машиностроение, 1984.

14. Райсман X. Ф. Метод конечных разностей как вариант метода конечных элементов. — В кн.: Труды Ленинградского кораблестроительного института, вып. 85. — М.: Стройиздат, 1973.

15. Wu Chang-chum. A new numerical method for analysing bent elastic thin plates-—the discrete operator method. — Scientia sinica, 1975, vol. 18, N 5.

16. Reissman C., Huhn H.—J. Die Anwendung des finiten Bi-lanzverfahrens die orthogonal—anisotropen Potentialfeldern. — Wissen-schaftliche zeitschrift der universitat Rostoc. Math. — naturwiss. R„ 1974, vol. 23, N 3.

17. Снигирев В. Ф. К расчету плоских панелей методом конечных элементов. — Изв. вузов «Авиационная техника», 1977, № 4.

18. Постников А. А. Расчет элементов конструкций из композитных материалов методом дискретных операторов. — В кн.: Материалы второй всесоюзной научно-технической конференции. — Прочность, жесткость и технологичность изделий из композиционных материалов. — Ереван, 1984, т. 3.

19. Тимошенко С. П., Войновски й-К р и г е р С. Пластинки и оболочки. — М.: Наука, 1966.

Рукопись поступила 22/У 1985 г.