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

Вестник Сыктывкарского университета. Сер Л. Вып. 7. 2007

УДК 539.3

ОБ ОДНОМ АЛГОРИТМЕ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ С НЕЛИНЕЙНОСТЬЮ ТИПА КАРМАНА

Михайловский Е.И., Миронов В.В., Кузнецова Ю.Л.

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

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

1.1. Трансформирование краевой задачи в систему интегральных уравнений

Рассмотрим прямоугольную в плане цилиндрическую пластину (рис.1). Уравнения равновесия такой пластины по теории К.Маргера имеют вид [1]

д?

ОА2ии + Я— = К^п + А(Ф, Ч

д2ии 1

Здесь ио, Ф искомые прогиб и функция напряжений, нормальная удельная нагрузка в пересчете на срединную поверхность; И = = ЕН3/12(1 — и2); Н, II- толщина пластины и радиус срединной поверхности, Е1, у - модуль Юнга, коэффициент Пуассона материала пластины;

~ д2Ч*д2уо д2Ф д2У0 а2Ф д2У0

у ' - д(р2 - 2 д£д(р д^д(р + д(р2 '

© Михайловский Е.И., Миронов В.В., Кузнецова Ю.Л., 2007.

; , /э2 а2 \2 Д(,= 1аё + 8?)(>-

Рис.1. К расчету цилиндрической панели.

Имея в виду использование в дальнейшем квадратурной формулы Гаусса-Лежандра, произведем в уравнениях (1.1) замену переменных по формулам

х=21-1, „=?£-1, (1.2)

Со (ро

в соответствии с которыми —1<х,у<+1. Внешнюю нагрузку задаем формулой

/$.(*■¥)£ п.

Яп = < , (1-3)

I о, (х,у) ^а

где

0£ = {(ж, у) : |ж| < е, \у\ < е}. В координатах х, у система уравнений (1.1) принимает вид

д2т Eh.Il ЕН . . . ^

Д2Ф = -—гъихх--^ о го). (1.4)

4£о г&да

где

дх2 ду2 джду <9у2 <9:е2'

/ а2 д2 \2 <92Ф

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

(1.5)

К р.

да* го

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

гу = О, Муу = 0, В*п = П Ф = 0, 4 = 0. (1.6)

При этом условия (1.6)2 можно раскрыть так:

Мп = -40 ( + ) при X = ±1,

М22 = —41)

( д2ио

+ V-

ср20ду2 д2ио

при у — ±1.

Учитывая, что в силу условия (1.6)1 выполняются равенства

д2ии ^ ^ д2ш

= 0 при х — ±1, —— = 0 при у — ±1,

ду2

из (1.7) получаем д2ио

дх2

= 0 при х — ±1,

д2

IV

= 0 при у — ±1.

дх2 ' ду2

И наконец, условие (1.6)4 в координатах х, у имеет вид [1]

1

х = ±1 : = е22 =

и) Я

(а2Ф а2Ф

I —Г—7 — V-

ЕКК2 \Цдх2 1

ср20ду'

дъи\2 -5Г =0;

±1 = 4

2 Я? \ду )

_ 1 / д2Ф д2Ф 611 ~ Е/1Я2 \vldy2 ~ и£дх2

1

2\дх

(1.7)

(1.8)

Отсюда, с учетом условия (1.6)i, (1.8), (1.6)з, (1.9), будем иметь

d2w дЧ

х = ±1: w = О, -—- = О, Ф = 0, — = 0,

охг oxz

d2w дЧ

У = ± 1 : w = 0, — = 0, Ф = 0, — = 0. (1.10)

Учитывая, что левые части уравнений (1.4) записаны с использованием одного и того же оператора (1.4')2 и принимая во внимание одинаковый вид граничных условий для ги, Ф, для обращения операторов левых частей этих уравненний можно использовать единую функцию Грина

A2G(x, у; а, /3) = 8(х - а)8(у - /3); (l.ll)i

d2G

G — 0, ——- = 0 при х — ±1, ох2

d2G

G = 0, — = 0 при у = ±1. (1-11)2

Граничные условия задачи (1.11) будут выполнены, если функцию Грина искать в виде ряда

G(x, У] а, /3) = Gii cos ~Y C0S ~2~" ^ '

i,j=1,3,..

Соответственно правую часть уравнения (l.ll)i представим в виде формального (расходящегося) ряда

6(х - а)6(у = cos cos (L13)

¿,¿=1,3,..

^ /1 ю\ mux птгу

Умножив обе части равенства (1.13) на cos—-—cos —и выполнив

интегрирование по области получим

г ттта Tin [3

ошп — cos-cos-.

mn 2 2

Подставляя найденные коэффициенты в ряд (1.13), а затем, полученный ряд и ряд (1.12) в уравнение (l.ll)i, после несложных преобразований придем к формуле

ma jiif3 тх jiry

-I а cos-cos-cos-cos —

lb ^ ^ 2 2 2 2

G(*,.«.й^ЕЕ —2 J -P-14)

i,i=l,3v ( _ +

£ &

Используя функцию Грина (1.14) краевые задачи (1.4), (1.10) можно привести к следующей системе интегро-дифференциальных уравнений:

-К4 С

у) = 77^ / 2/5 Р)Яп(а, /3)с1а(1(3—

1

i? С

/ у; а, /3)Фжж(а, f3)dadf3+ Jnx

С№D Jni

G(x, у; а, /3)[ФЖЖ^ - 24*xywxy + /3)dad/3,

EhR

Ф(х, у) = - I G(x, y\a, /3)wxx(a, /3)dad/3-

Eh

/ - (1.15)

Систему (1.15) можно решать итерационными методами с использованием конечно-разностной аппроксимации производных и формул численного интегрирования. Учитывая, что вблизи криволинейных кромок реализуется краевой эффект, и что при малых значениях £ в центральной части пластины напряженное состояние имеет большой показатель изменяемости [2], решение задач с применением формул численного дифференцирования требует использования густой равномерной сетки, что неизбежно связано с увеличением погрешности. Однако, в силу того, что в правых частях уравнений (1.15) фигурируют лишь вторые производные от искомых функций, эти интегро-дифференциальные уравнения можно трансформировать в интегральные. Действительно, дифференцируя дважды каждое из уравнений соответственно по „хх",

ихУи 1 „ууи, получим

R4 f R С

wxx = / Gxxqdadf3 - —— / Gxx4>xxdadf3+ Jn 4D& Jn

- 2Фxywxy + Фyywxx\dad/3,

KoVo-U Jn

R4 Г R f

Wxy = TTTn / Gxyqdadf3 - / Gxy^xxdadfi+

Jn Jn

XX 1GD

1 г

2I ^xyxx^yy 2ФxywXy 4fyywxx]dad^ ¥оU Jn

gvlD

i?4 f R

Wyv = 16D Jn GyyqdadP ~ ¡^2 J Gyyyxxdadf3+

+ с2 \ п [ Gyy[4fxxWyy - 24fxywxy + 4fyywxx]dadß,

^xx = [ Gxxqdadß - f Gxx[wxxwyy - w2xy]dadß,

44o Jn Чо^о Jn

EhR f Gxyqdadß - f Gxy[wxxwyy - w2xy]dadß, Jn So ^Pо Jn

vT; =

^ XII

EHR С Eh С

Фyy = / Gyyqdadß - —-- / Gyy[wxxwyy - w2xy]dadß. (1.16)

4So Jn ^o^o Jn

1.2. Переход к системе алгебраических уравнений

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

/ Guvqda dß = wquv(x, у), Jn

(под uv следует понимать „хх", „ху" или »yyL\ впредь для сокращения записи будем использовать цифровые обозначения: u = х => и = 1, и = у и = 2, v = x=>v = l, v = y=>v = 2) полученной системы (1.16), воспользуемся квадратурной формулой Гаусса-Лежандра [3]

/ /(х, y)dxdy = / / /(х, y)dxdy = J-iJ-l

n m

^ ^ ^ ^ ^4-г/п B^Lm f (KXl/n, У ¡лт ) ; (^'l)

г/=1 ¿1=1

где -1 < < X2n < ... < ЖПп < 1, -1 < Ulm < У2т < ••• < Утт < 1 корни полиномов Лежандра 1/п(ж), Lm(y) (узлы квадратурной формулы Гаусса-Лежандра)

4 2(1 - xln) 2(1 - У;та)

^-z/n __9/7" \\9 5 -ЕЗцт

п2 (Ьп_1 (х^))2' га2 (х^ш))2

коэффициенты квадратурной формулы Гаусса-Лежандра. Придавая переменным х и у значения в узлах квадратуры Х{ — — ХгП (¿61: п), У1 — Ууш 0 ^1: т) получим систему из 6тп уравнений, относительно 6тп неизвестных гиоиу{х11у^)1 Фиу{р^иУз)ч (^ 6 1 : п, ] £ 1 : ш, и = х, у, г; = х, у)

д4

xx^yy

¿,¿=1,3,-

Фад^ (^г? 2/j) ^ ^ ^ Al/nBjJjrnGuv (yXii yj ? OIvï (^i/J

z,j= 1,3,..

Eh

'Ш

y^ y^ AvnBpmGuvixiiyjiav, P^WxxWyy - w2xy](o^,/ЗД (2.2)

¿,¿=1,3,..

Значения для коэффициентов при заданной нагрузке q

считаются известными. В самом деле, учитывая формулу (1.3), имеем

q / ч _ 4Д4 Qq ул

7Г2D 4^ / / 2 /2 4 2

kulv kiïXi U7T fl^Vi vtt.t cos(--1--) cos(--1--

k,l=1,3,.. [ k I l I 72 +

Ço ^o

kira Г^ Itï/3 1 kira ^

£C0S(—)da y^cos(—= ^çysin(—) ^

2sin(*fs)sin(*f)

kirs lire

2 2

1 . Лтг/З,

sm

(2.3)

(2.4)

1.3. Построение итерационной схемы стационарного метода Ричардсона (СМР)

Введем обозначения ТУ = [гиц(1,1), ..,1/;11(п,77г); 11712(1,1), ••, «^(п, ш); «^2(1,1), »,^22(п,т)],

УУ9 = [^11С1^ !). ••> МгЛп, «'У1,1), ••> ^\2{п, т); ыд22(1,1),.., тд22(щ т)],

Ф = [Ф и(1,1),.., Фц(п, т); Ф12(1 ? 1)^ Фгг^? т); Ф22(1 ^ 1)^ Фгг^? га)];

(3.1)

^ = ^(И7, Ф), С = (^(И7) - нелинейные правые части в системе (2.2)

Система (2.2) приводится с учетом обозначений (3.1) к следующему виду:

IV = }¥д + Р(\¥, Ф),

Ф = о + G(W).

(3.2)

Положив ]Уо — Фо — 0 (формулы нулевого приближения), запишем итерационную схему СМР

ФЛ+1 = (1-г)ФЛ + гС(ВД, А; = 0,1,..

ИЪ+1 = + (1 - т)ИЪ + тГ(\¥к, ФЛ+1), к = 0,1,.. (3.2)

2. Численный эксперимент

2.1. Об учете пологости панели

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

Л4 Я

^ хх 1

л 2 Т

Л Ф = -ТРГ^-

Подействовав на первое уравнение оператором Ад, и произведя вычитание из него второго уравнения, придем к классическому уравнению типа В.З.Власова [4]

д4 12(1-//2)Д2а4^ Л4 л2 А" ги Н-------=-А а

Граничные условия остаются прежними (см. форм. (1.10)).

Для решения полученной линейной задачи уместно применить метод двойных тригонометрических рядов, согласно которому, для удовлетворения граничным условиям, необходимо принять следующие представления для прогиба и нагрузки:

Ш \ ^ ^ ( Ыкт \ ъкх тгту

Подставив эти представления в последнее уравнение, получим следующую формулу для вычисления коэффициентов Фурье:

Я4 /я-2\2 (к2 т2Х 2 +

ш =_шр\22) ^_

кт + /тг^Л2 12(1 -

22) ^о ^о/ \ 22 / 16£04/г2 \ 22

(неучет подчеркнутого слагаемого соответствует решению для задачи об изгибе плоской пластины).

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

Из приведенных рисунков видно, что рассматриваемые прогибы удовлетворительно согласуются лишь для очень пологих (<р0 < 0.5 ^ 30°) и относительно толстостенных (h/R > 0.1) панелей.

Рис.2. Графики прогиба цилиндрической панели (ии1) и соответствующей плоской пластины (ш11) при I = Я = 50 см, (2о — Ю кг, £ = 1.

2.2. Об учете краевого эффекта в цилиндрической панели

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

Рис.За. Графики изгибающих моментов при I = Я = 50 см, (¿о = 100 кг, £ —

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

Рис.36. К оценке ПКЭ по нелинейной теории (графики изгибающих моментов при I = Я = 50 см, (¿о = 100 кг, г = 1, сро = 7г/8, Н = 2 см).

2.3. Об учете нелинейности

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

Рис.4а. Графики для прогиба и изгибающих моментов при I = Я = 50 см, г = 0.5, (ро = 7г/8, Н = 1 см.

Рис.46. Графики для прогиба и изгибающего момента Мц при I = Я = 50 см, г = 0.3, сро = 7г/4, Н = 3 см.

Графики на рис. 4 представлены для значений нагрузки фо? обеспечивающей сходимость итерационного процесса, т.е. при больших (меньших) значениях (^о итерационный процесс расходится. При этом, индексом I отмечен график для величины, рассчитанной по линейной теории.

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

Вычисления проводились с использованием 20 членов по каждой переменной в рядах для функции Грина и ее производных и при 30 узлах квадратурной формулы Гаусса-Лежандра в каждом направлении на неравномерной сетке. Использовалась итерационная схема стационарного метода Ричардсона, в которой параметр т каждый раз подбирался так, чтобы обеспечить сходимость итерационного процесса.

•6000-

300

Рис.4в. График изгибающего момента М22 при I = Л = 50 см, г = 0.3, сро = 7г/4, Н = 3 см.

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

Для изображения изгиба цилиндрической панели при различных значениях нагрузки, приведем соответствующие поверхности при 0 = 100 кг и д0 = -Ю0 кг.

Рис.5а. Поверхность прогиба при 1 = В = 50 см, г = 0.5, сро = 7г/4, Н = 1 см, (¿о = 100 кг.

Рис.56. Поверхность прогиба при I = R = 50 см, г = 0.5, ipo = 7г/4, h = 1 см, Q0 = —100 кг.

Литература

1. Михайловский Е.И. Математические модели механики упругих тел. Сыктывкар: Изд-во Сыктывкарского ун-та, 2007. 516 с.

2. Михайловский Е.И., Черных К.Ф. Развитие механики оболочек в трудах школы академика В.В.Новожилова // Успехи механики. 2003. Т.2. №3. С.87-126.

3. Абрамовиц М., Стиган И. Справочник по специальным функциям. М.: Наука, 1979. 832 с.

4. Власов В.З. Избранные труды. T.l. М.: Изд-во АН СССР, 1962. 528 с.

Summary

Mikhailovskii E.I., Mironov V.V., Kuznetsova J.L. About one solution algorithm for solving of nonlinear boundary problem Karman-type.

On the basis of Green's well-known function, an algorithm of edge problems with Karman-type non-linearity to the corresponding system of algebraic equations. The algorithm is realized on tne examples of simply supported opened cylindrical shell.

Сыктывкарский университет

Поступила 31.10.2007