Интегро-интерполяционный метод решения плоской задачи для композита, армированного семейством криволинейных волокон Текст научной статьи по специальности «Математика»

Journal of Siberian Federal University. Engineering & Technologies 1 (2009 2), 112-120

УДК 539.3+539.4

Интегро-интерполяционный метод решения плоской задачи для композита,

армированного семейством криволинейных волокон

Н.А. Федорова*

*

Сибирский федеральный университет, 660041 Россия, Красноярск, Свободный, 791

Д

Получена 10.05.2008, окончательный вариант 17.06.2008, принята к печати 10.03.2009

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

Ключевые слова: композит, криволинейные волокна,, армирование.

1. Введение

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

2. Постановка задачи

Рассмотрим плоскую задачу упругости для среды, армированной одним семейством волокон. Пусть армирование выполнено волокнами постоянного поперечного сечения. Для описания композита используется структурная модель [4]. Введем обозначения: интенсивность армирования семейства волокон у), компонент тензора деформаций е^ (х, у), деформацию в волокнах семейства е\(х, у), напряжение в волокнах семейства - о\(х, у), осредненные напряжения обозначим через а^ (х,у), где х,у - декартовы координаты, <р(х,у) - угол армирования, индексы г, ] = 1, 2. В дальнейшем при обращении к перечисленным функциям для краткости аргументы будем опускать. Тогда структурную модель запишем в виде:

(wi/il),l + (W1Z12X2 = 0,

(1)

£11 l2i + £22/22 + 2e 12/11/12 = £?, £11,22 + £22,11 = 2£12,12-

(2) (3)

$

Corresponding author E-mail address: [email protected] © Siberian Federal University. All rights reserved

Использованы обозначения: ef = afT, ei = ei + ef, lii = cosli2 = sinaf — коэффициент линейного расширения материала семейства волокон, T = const (T — температура, заданная в единицах Си). Символы 1, 2 означают частное дифференцирование по координатам x, y соответственно. Правая часть в (2) учитывает как случай равнодеформированного волокна (ei = const, e0 = const + ef), так и случай нерастяжимого £ =0, e0 = ef) .

Заметим, что уравнение (1) - это условие постоянства сечений волокон, (2) - деформация в волокнах семейства, правая часть уравнения учитывает деформации eT, вызванные постоянным полем температур, (3) - уравнение совместности деформаций.

Осредненные напряжения aj (x,y) запишем в виде:

aij = (1 — ui)acj + ak Wi liilij. (4)

Соотношения (4) есть определение силы, действующей на слой композитов как суммы сил, создаваемых связующим материалом, и сил, создаваемых армирующими слоями. В (4) напряжения в связующем определяются по формулам:

E E

(1 _ v 2) (£ii + V£ÓÓ _ ac(1 + v)T)' aij = (1 + v) £ij' = 3 - г = 1 2)-

Здесь E, v, ac - соответственно модуль Юнга, коэффициент Пуассона и коэффициент температурного расширения связующего материала, Ei— модуль Юнга волокна. Напряжения aij должны удовлетворять уравнениям равновесия:

aii,i + ai2,2 = bi, (i = 1, 2). (5)

Правые части в (5)

bi = —((1 — ui)pc + W pi)fi

являются компонентами массовой распределенной нагрузки по направлениям прямоугольной декартовой системы координат; pc, pi— массовые плотности материалов связующего и волокона ; fi— компоненты удельной распределенной нагрузки, действующей на единицу массы.

К системе (1) - (5) присоединяются граничные условия на контуре. Уравнение контура Г задано в параметрическом виде: x = y(s), y = ^(s), s—некоторый параметр. Пусть на контуре Гр заданы статические условия с нормальными и касательными усилиями pn(s),pT(s) соответственно:

аип2 + a22n2 + 2ai2nin2 = Pn(s), (a22 — ац)пп + a^n0 — n2)= Pt(s). (6)

На другой части контура Ги заданы кинематические условия для перемещений ui,u2:

ui(r„)= u?(s), и2(Ги)= u2(s). (7)

В (6) pn(s), pT(s)— известные функции, ni = cos в, n2 = sin в, в— угол, задающий направление внешней нормали к Гр. С учетом (4) граничные условия (6) принимают вид

wiai cos2(^ — в) + (1 — Wi)[m3(eii + ve22 — LT) cos2 в +

+тз (e22 + ven — LT) sin2 в + m4ei2 sin в cos в] = Pn(s), (8)

wiai sin2(^i — в) + (1 — Wi)m3(e22 — eii + v(eii — e22)) sin2в +

+2(1 — wi)m4ei2 cos2в = 2pT(s).

a?

В (8) использованы обозначения для констант

LT = ac(1 + v)T, тз = Emi, m.4 = Em^, mi = -^, m-2 = -•

1 — v2 1 + v

Интенсивность wi задаем на той части Гш контура, где волокно входит в конструкцию:

wi (Гш) = w¡ (s) (9)

Ограничение для интенсивности армирования имеет вид

0 <wi < 0, 7•

3. Семейство равнонапряженных волокон

Тогда правая часть в (2) примет вид ei = ei + еT, напряжение в волокне ai = Const• В этом случае система (1),(2),(3) - замкнутая система восьми уравнений относительно восьми неизвестных

,a22,ai2,eii,e22,ei2^

Переформулируем названную систему в перемещениях, выразив напряжения через деформации, а деформации через перемещения, используя соотношения Коши. К системе присоединим граничные условия на контуре Ги : м(Ги) = u(s) либо граничные условия (6), где деформации выражены через перемещения.

При подстановке напряжений через перемещения в уравнения равновесия нужно вычислить следующие выражения:

(aiwi cos2 ф) i + (aiwi cos ф sin ф) 2,

(aiwi cos ф sin ф) i + (aiwi sin2 ф) При их вычислении выделяем слагаемые вида

(wi cos ф) i + (wi sin ф) 2,

которые в силу условия постоянства сечений волокон (1) равны нулю, что упрощает данные выражения, и они становятся соответственно равными:

aiwi (cos y(cos ф) i + sin y(cos ф) 2)

aiwi(cosy(sinф) i + siny(sinф) 2)

Для формулировки системы удобно ввести z = tg ф(х,у)• Тогда производные от функции, задающей угол армирования, вычисляем по формулам

z, i ; = Z, 2

ф, i = ; ф,2 =

1 + 1 + z2'

В итоге система для равнонапряженного семейства волокон относительно искомых переменных п\, щ записывается в виде:

г 1 г 2

^11 + 2 + ' 2 + ' 2 ) = О,

1 + г2 1 + г2

Щ 2 + 1 и1 1 + и 1 - гщ 2 - е0(1 + г) = О

' -у ' ' ' -у

^ТОз(иМ1 + ^2, 21) + 2шт4(и1,22 + Щ2,12) + тз(-^1д)(и1д + 2 - Ь ) +

+ 1 Ш4(-^1, 2 )(и 1, 2 + и2, 1) + ^1ш1(-(1 +г%2)2 1

(г, 1) -

(1 + г2)2

(г, 2)) = Ъъ

иш3(и2,22 + ^1,12) + 2^Ш4(и1,21 + и2,И) + Шз(-^1,2)(и2,2 + ^1,1 - Ь ) +

1 1 г

+ 2т4(-ш1,1)(и 1,2 + и2,1) + 0"1 ^1((1 + г2)2 (г,1) + (1 + г2)2 (г,2)) = Ъ2•

(10) (11)

(12) (13)

В системе использованы обозначения ш = 1 - Ш1.

Первые два уравнения системы содержат производные только первого порядка, поэтому чтобы исследовать тип системы, продифференцируем эти уравнения по одной из координат. Затем построим характеристический многочлен [5] для данной системы

Р (Л) = ¿вЬ(А11Л2 + 2АЛ + А22),

где

А11 =

1-

О О

гш1 1 + г2 О О О

00

1

- 1

г

штз 0

ШШ4

Т" )

О

А

12

Ш1

1 + г2

О

0

1 О

ш(итз + т)

О

-г

( + т4 ) О

\

,А

22

/ 0 0 0 0 \

0 0 0 0

0 0 шт4 2 0

V 0 0 0 штз /

При исследовании типа системы установили, что характеристический многочлен (полином) тождественно равен нулю, что означает вырождение типа системы. Ее решение представляет определенные трудности и требует специального подхода. В данной работе находим частные решения после введения угла армирования как заданной функции координат, что позволяет определить интенсивность армирования из уравнения (10). Уравнение (11) рассматриваем после вычисления и1, и2 как ограничение на условие равнонапряженности. Решаем уравнения (12),(13) совместно как систему относительно и1, и2. Ее характеристический полином имеет вид

Р (Л) = ¿вЬ(А11Л2 + 2АА + А22)

где

А

11 I штз 0

0 0, 5шт4 I '

2

г

0

. 12 / 0 w(vm3 + 0,5m4) \ 22 ( 0,0

у ш(ушз + 0, Бш4) 0 у ^ 0 шшз

Уравнение Р(А) = 0 после подстановки тз, Ш4 через модуль Юнга Е и коэффициент Пуассона V запишется после некоторых упрощений как биквадратное уравнение

А4 + 1=0

Находим четыре корня биквадратного уравнения:

\/-(7 - 1)(372 - 1 + 107) ± ^ЭЙ + 7872 + 6072 - 1Б + 127)

А< = -2(7-1)-, (к = 14)

Исследование выражения под радикалом (974 + 7872 + 6072 - 1Б +127) как функции от V показало, что оно имеет отрицательный знак на интервале [0; 0, 38]. Это означает, что корни Аи - комплексно сопряженные на этом интервале. Рассматриваемая система (12), (13) имеет эллиптический тип. Заметим, что для многих материалов коэффициент Пуассона принимают равным от 0, 2Б до 0, 3 [6],[7]. Следовательно, возможные значения коэффициентов Пуассона попадают в указанный интервал, где корни Аи - комплексно сопряженные. Для построения численной схемы запишем оба уравнения в дивергентном виде. Предварительно введем следующие обозначения:

У1 = ^(х, у), У2 = ¿(х, у), ¿1 = (1 - уг)тз, ¿2 = (1 - у1)т4/2. (14)

Здесь у 1, у2 — известные функции координат, 41,42- неизвестные функции. Заметим выполнение следующих соотношений для слагаемых, входящих в первое уравнение

(1 - у1)тзМ1Д1 + тз(-у1д)м1д = (г^д)д, (1 - у1)тз7М2,11 + тз(-у1д)7М2,1 = (^1^2,2)д,

(1 - у1) , 1 ( ) ( )

т4М1,22 + 2 т4(-у1,2 )и1,2 = (¿2^1,2),2,

2 -±,22 • 2 (1 - У1) , 1

2 -m4M2,12 + 2 m-4—У1,2 )и2,1 = (^2U2,l),2. В итоге оно запишется

(г1«1,1)д + (721 М2,2)д + (¿2^1,2),2 + (¿2^2,1),2 + =0, (15)

2

тт. Т , ГТ , ! у2 у2 л

Р1 = -&1 + тзу1,1ь + а1у1(-(ТГу^у2,1- у2,2).

Используем соотношения для второго уравнения

(1 - у1)тзП2,22 + тз(-у1д)м2,2 = (¿142,2),2,

(1 - у1)тз7П1 12 + тз(-у 1,2)741,1 = (7^1^1,1),2,

(1 - у1) . 1 ( ) ( )

-2-т441,12 + 2 т4(-у1,2 )41,2 = (¿241,2),1,

(1 2 у1) т442,11 + 2 т4(-у1,2 )42,1 = (¿242,1)д.

В итоге второе уравнение запишется

(¿142,2),2 + (7^1^1,1)^ + (¿241,2),1 + (¿242 д)д + ^2 =0, (16)

Р2 =-&2+тзу1,2ьТ+а1у1 ( у2,1 + (ггу^ у2^.

4. Краевая задача для прямоугольной пластинки

Пусть граничный контур Г— граница О, где О— прямоугольная пластинка со сторонами, параллельными координатным осям О = {(х 1,х2)|, 0 < хг < г = 1, 2}. В этом случае определенный выше угол в = п/2. Кинематические условия (7) не меняются. Статические условия (6) после замены компонент деформации на перемещения по формулам Коши примут вид

2

а 1г2

Рп(в) — 1 + 2

П2 2 + »411 = -+ Ьт, (17)

штз 2рт (-) + ^

+ ^ ' 2(1 + z2)

Ui 2 + VU2,1 = -.

— M

При построении разностной схемы для системы (15),(16) поставим задачу отыскания вектора u = (ui, U2) с граничными условиями на контуре в общем виде

ziUi,i + Z1VU2, 2 = —fi(xi, 0), (18)

Z2U1,2 + Z2U2,i = —fi(0, X2), Z2Ui,i + ziVU2,2 = —f2(xi, di), Z2Ui 2 + Z2U2 ,i = —f2(d2,X2).

5. Свойства дифференциального оператора поставленной задачи

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

Выпишем функционал /(и, V) = (Ьи, V) + (Е, V)

I(u, v)

д dUi д dU2, д dUi 9 dU2

я—(Zi ä— )vi + я—(ziv^— )vi + я—(z2 я— )vi + ^—(Z2^— )vi + Fivi

.dxi dxi dxi дх2 дх2 дх2 dx2 dxi

dXi dX2 +

+

д dU2 д dUi д dUi д ди2 „

Я-(zi^-)V2 + я-(Zi^-)V2 + я-(Z2^-)V2 + я-(Z2^-V + ^2^2 dxidx2,

19X2 9X2 9X2 9XI 9XI 9X2 9xi 9xi

где и = (41,42), V = (г>ьг>2), Е = (^ъ^).

Используя тождества, возникающие при нахождении интеграла по частям, например,

д ди1) д ди1)+ д41 дУ1

дх1 дх1 дх1 дх1 дх1 дх\

запишем функционал в виде ди1 дУ1

I(u, v)

9U2 9vi дщ 9vi 9U2 9vi ди2 9v2

ZiЯ—Я--+ Zi^-—---+ Z2-—---+ Z2-—---+ Zi-— ---+

9xi 9xi 9x2 9xi 9x2 9x2 9xi 9x2 9x2 9x2

д—i dv2 д—i dv2 du2 dv2 „ „

+ziv-—Я--+ Z2ä—я--+ —---+ *ivi + F2V2

dxi dx2 dx2 dxi dxi dxi

dxidx2+

(19)

+

+

fi(xi, 0)vi + fi(xi,d2)vi + /2(xi, 0)v2 + /2(xi,d2)v2 dxi+

fi(0, x2)ui + fi(di,x2)ui + /2(0, x2)U2 + /2(di,x2)u2 dx2.

Из формул (19) легко видеть, что при замене и и V /(и, V) = /(у, и), т. е. функционал симметричен.

Исследуем свойство положительной определенности оператора Ь, выпишем (Ьи, и) :

L(u, u)

( dui )2 ( du2 )2 . ( dui )2 , ( du2 )2 , du2 dui . z1(^) + + z2( Я ) + z2 + Z2-— +

dxi dx2 dx2 dxi dxi dx2

dui du2 du2 dui dui dui

+z2ä—Я--+ zi^-—---+ zi v-—-—

dx2 dxi dx2 dxi dxi dx2J

dui )2 + zi( dui )2 + 2zivdui dui + z2( dui + dui )2

dxi dx2 dxi dx2 dxi dx2

dxidx2 = dxidx2.

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

(дп\ ди2 )2 ^ о 2 дж1 дх2 > '

поскольку Х2 > 0, что следует из соотношений (14) и ограничений на значения интенсивно-стей армирования, аналогично Х1 > 0. Выражение

,dui,2 ,д—2^2 dui du2

zi( я-1 )2 + zi(^)2 + 2ziv^du2 > 0,

dxi dx2 dx1 dx2

d—i d—2 d—i d—2 (, 2 2 если dx~> 0- В противном случае, если д^< 0' тогда z^(-x) + () -

-2viduui d—2 i)>:i( (dui »2+(d—2 >2 - 2

i 2 i 2

d—i д—2 \ ( d—i д—2

dxi dx2 )= < dxi dx2

0.

Оценка следует из того, что коэффициент Пуассона 0 < V < 1.

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

d

о

о

G

G

2

6. Построение разностной схемы

Решение задачи (15), (16) с граничными условиями (18) в силу дивергентной формы записи системы, Ьи + Е = 0, где Ь— дифференциальный оператор системы, Е- правая часть, должно доставлять минимум следующему функционалу

/(и, и) = (Ьи, и) + (Е, и),

где и = (и1,и2), Е = (^1,^2). Выпишем этот функционал. Согласно методу Ритца [8] аппроксимируем пространство W2(G), в котором ищутся решения, конечномерным подпространством V и назовем приближенным решением задачи вектор и = (^1,^2),

гог € V, который минимизирует функционал /(и, и) на пространстве V. Этот функционал можно преобразовать так:

/(и, и) = Ш(и) + / / (Г1и1 + ¥2и1)в,Х1в,Х2 + ио Jо

/■ й1

+ / [Л(х1, 0)^1 + ¡1(х1,й2)и1 + /2(^1, 0)и + /2(^1, ¿2)42^x1 +

о

+ / [Л(0,Х2)и1 + /1 (¿1,Х2)и1 + /2(0, Х2)и2 + /2(¿1, Х2)и2]^Х2, (20)

о

где Ш(и)— энергия упругой деформации преобразуется к виду

Ш(и) = 1/ / Ы(иМ)2 + (и2,2)2)+ ^(К2)2 + (и2Д)2)+ (21)

+2^1^и1 1и2 2 + 2х2П2 1ги1 2]^Х1^Х2.

^ J а

Введем в области О прямоугольную равномерную сетку. Разобьем область на прямоугольные ячейки со сторонами ^1,^-2 и вершинами в узлах сетки. Каждую ячейку в свою очередь разобьем прямой, проходящей через ее противоположные вершины, на два треугольника. Обозначим левые верхние треугольники через Д+[г], а правые нижние - через Д-[г]. В качестве подпространства V пространства W2(О) возьмем пространство непрерывных в области и линейных над каждым треугольником Д^ [г] функций.

Приближенное решение задачи (15),(16) будем искать в виде и = (V!, V2) :

N1 N2

'иг ^У^, Уг3132 П3132, (22)

31=0 ¿2=0

где Пз1з2 — базис пространства V. В качестве базиса в V можно взять совокупность непрерывных в О, линейных над каждым треугольником Д^ [г] функций, каждая из которых отлична от нуля лишь в одном узле сетки. При выбранном базисе коэффициенты угз1з2 в (22) имеют значение приближенного решения в узлах сетки.

Дифференцируя /(и, и) по переменным Уг,^^ и приравнивая первые производные нулю, получим следующую систему алгебраических уравнений для определения значений решения в узлах сетки Уг31 32 :

(^1х1 X + ¿^^Х^ХТ + (¿1^2X2 )Х1 ) + + (^1X2 )Х2 + + (^2х1 )Х2 ) = Ф1(Х1,Х2 Ъ (23)

(^2х2 )Х2 +2((*1"'и1Х1)Х2 + (¿1^1x1 )Х2) + + 2 ((^2^Х2)Х1 + (^1х2 )Х1 Ж^2х1 )хГ = ^Ф2(х1,х2).

В (23) Ф1(х1, Х2), Ф2(Х1, Х2)— сеточные выражения, содержащие значения правых частей; в уравнениях введены сеточные операторы дифференцирования по следующему правилу:

_ = (иг — иг-1) = (иг+1 — иг) (24)

их1 — , , их1 — ; , (24)

П1 П1

_ = (иг — иг-1) = (иг+1 — иг)

иХ2 ; , иХ2 ; .

П2 Пл.

Заметим, что для достаточно гладких правых частей погрешности аппроксимаций имеют вид = o(h2 + [8, 9].

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

Список литературы

[1] Немировский Ю.В., Кургузов В.Д. Прочность и жесткость стеновых железобетонных панелей со сложными структурами армирования // Известия вузов. Строительство. 2003. № 2. С. 4-11.

[2] Немировский Ю.В., Федорова Н.А. Моделирование деформирования плоских авиационных конструкций, армированных семействами криволинейных волокон // Вестник Сиб. гос. аэрокосмич. ун-та. Вып. 6(13). Красноярск, 2006. С. 38-44.

[3] Федорова Н.А. Решение плоской задачи упругой среды, армированной тремя семействами волокон // Вычисл. технологии. 2005. Т. 10. С. 90-100.

[4] Немировский Ю.В., Янковский А.П. Рациональное проектирование армированных конструкций. Новосибирск: Наука, 2002. 487 с.

[5] Бицадзе А.В. Некоторые классы уравнений в частных производных. М: Наука,1981. 448 с.

[6] Тимошенко С.П., Гудьер Дж. Теория упругости. М: Наука, 1979. 550 с.

[7] Композиционные материалы. Справочник / В.В. Васильев,В.Д. Протасов и др. М: Машиностроение, 1990. 512 с.

[8] Самарский А.А. Теория разностных схем. М: Наука, 1989. 614 с.

[9] Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М: Наука, 1989. 425 с.

Integro-Interpolational Method of Plane Problem Solving for Reinforced by Family of Curvelinear Fibers Composite

Natalia A. Feodorova

Siberian Federal University 79 Svobodny, Krasnoyarsk, 660041 Russia

It is constructed the permitting system of equations of plane problem for environment reinforced by the family of curved fibers. The problem articulated in the movement for family of equal strained fibers when the angle of reinforce is the known function of coordinates. It is displayed that the system can be written in divergence form on the movement. Based on the Ritz method the numerical pattern adapted to the singularities of the task is received.

Key words: composite, curvelinear fibers, reinforcement.