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

УДК 539.3

ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ СИСТЕМЫ POLE ДЛЯ ИССЛЕДОВАНИЯ ЗАКРИТИЧЕСКО-ГО ПОВЕДЕНИЯ ПЛАСТИН

ЛИННИК А.Б.

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

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

1. Постановка и метод решения задач о закритическом поведении пластин

Для решения задач о закритическом поведении пластин нужно исходить из общей теории гибких пластинок, учитывающей одновременно напряжения в срединной поверхности и напряжения изгиба. Уравнения равновесия в этом случае имеют вид [1]

D V4w = ст h

x

32w 32w , 32w

—2 + сту—- + 2t——

dx2 У dy2 dxdy ’

AU = l(w),

(1)

(2)

причем операторы A,l определяются следующим образом:

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

Nn = Pi0 Sn = P20,

(3)

где Pi0 ,p2 — внешние нормальные и касательные к кромке пластины усилия.

Линеаризация исходной нелинейной системы уравнений базируется на основе итерационного метода. Систему (1), (2) запишем в виде

4h V w‘+1 = D СТ

d2wi+1

xi

+ а

dxz Уі ду

AUi+1 = ^wi+i)

32wi+1 d2wi+1

----+ 2ті-------i±1 (4)

dxdy , ( )

(5)

здесь і — номер итерации.

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

Можно показать, что вариационная формулировка (4) сводится к нахождению минимума функционала

l(w) - J{(Aw)2 - 2i _v)(w'11w'22 - ^2)+

Q

-(niw2 + 2Swiw>2 + ^w^j/DjdQ .

Предполагается, что усилия NbS,N известны из предыдущего шага итерации (для простоты записи индекс і опускается).

Рассмотрим вариационную постановку задачи (5), (3). Пусть функция w(x,y) известна. Принимая во внимание (3) и связь между напряжениями и перемещениями, краевые условия запишем в следующей форме:

Eh2 (шп1 + v'nm + v[v'Xl- шхт]) = Fi0, x є 3Q ,

1 -v2

2^Eh^ (lv'n - mu'n + U41 + v'im) = F2> x є сП. . (6)

Здесь

AU = - -2 [(i + v)grad divU + (i -v)au];

ґ

!(w) =

1A

2 dx

2

dw1 dx

+ v

dw

~\2

dy

1 -v d ! dw dw

2 dy ( dx dy

1 -v d f dw dw 1 d

(rdw^2

2 dx (9x dy ) 2 dx

dy

+ v

dw

dx

-,C\

dj

Fi0 =-

E

f2 =

E

2(1 + v)

(w2n +VW?X)+ pi0,

W'nW4 + p2 .

Так как краевые условия являются неоднородными, решение задачи (5), (6) будем искать в виде

U = U0 + Ui, (7)

Уравнения (1), (2) дополняются определенными граничными условиями. Условия для прогиба пластинки имеют такой же вид, что и для жестких пластинок [1]. Будем считать, что пластинка под-

где U0 = (u0,v^ — любая достаточное число раз дифференцируемая в q вектор-функция, удовлетворяющая неоднородным краевым условиям:

РИ, 2000, № 2

123

Eh

2

(uQ'nl + V0'nm + v[vq41 - Uq4^ = if, x єдП

1 -v

Eh 0 2^ + vj (1v0' n - mu0' n + u0' x1 + v0' xm) = F20, x ;

Uj = (u^vj) — вектор-функция, удовлетворяющая однородным краевым условиям:

Eh

(ujnl + vj' nm + v[vj' xl - uj' xm]) = 0, x edD.,

2

!-v Eh

r(lvr n - mur n + ur xl + Vj xm) = 0, x є 3Q..

+v)v

Подставляя (7) в уравнение (5), получаем

AUj = l(w) - AU0 . (8)

Можно показать, что оператор А задачи (5), (6) является положительным и, следовательно, задача (8) эквивалентна вариационной:

I2(U) = (AUbuO - 2[l(w) - AU0,U0]. (9)

После некоторых преобразований функционал (9) будет иметь вид:

= J{uU + V°2 + 2vu1,1V1,2 4 2—(uj,2 + +

Q

+ w'l(ul,l + УП,2) + w (V j, 2 +vuu) +

+ (j _ ^(uj,2 + +

+ + "0,2М,2 + v(V0,2uU + и0ДМ,2) +

2—(u0,2 + "0дХии + Vu)l}d^ +

+ J {[(w2 +VW22 ^-v)w jw'2 +

w2; +vw2]в'2 + ^ -V)w'lW'2Ю'lJvl -

5Q

+

- 2

—E~ (ulnFl0 + VjnF20)}i

dS,

где ujn = ujl + Vjm, Vjn = -ujm + Vjl.

При этом учтено, что направляющие косинусы l и m на границе 5Q области Q соответственно равны l = -ю j; m = -ю '2 , где ю(x) = 0 — нормализованное уравнение границы области [2].

Дискретизация функционалов осуществляется на множестве функций, удовлетворяющих по крайней мере кинематическим краевым условиям. Структуры решений для функции прогиба в зависимости от способов закрепления краев пластины выведены в работе [2]. Остановимся более подробно на выводе структуры решения для перемещений u и v, удовлетворяющей краевым условиям (6). Согласно методу R-функций представим искомый вектор перемещений U в форме

u - Фц +юФ2j , V = Ф^ +юф22 ,

(10)

где Фij (i, j = ^2) — неопределенные компоненты структурных формул. Используя методику, пред-

ложенную в [2], продолжим условия (6) внутрь области q :

E

(

\-V

да да

— Dju н---Djv + v

dx dy

да да

-т-ТГ" TjU

dx dy

= F

E (да D да D | да Т | да Т ^ F (11)

^7Djv“ТГ^и + ^_TjU + ДГТ1V l = F2.

2(i + v) ^ ox dy dx dy j

Здесь Fj = EC^F® ) F2 = EC^F®) — продолжение фун-

кций Fj0 и F-2 [2]. Дифференциальные операторы имеют следующий смысл:

df

0

_ , да df да df

Sx dx dy dy

Т f - I + da df

1 Sy dx dx dy

5Q

дП

9n

dr

Подставляя в (11) выражения (10) и используя свойства дифференциальных операторов

DjO = j, x є dQ Tjro = 0, x є dQ

получаем систему уравнений, которую можно разрешить относительно Ф2j, Ф22. Тогда (10) можно представить следующим образом:

Л .2 ^ „л \ ^ А

u = а

j-V2 да 2А + v)_ да

----Fj---^F2 —

E j Sx E 2 dy

Эю^| (дал2

+ фц +

+ ю{-D1Ф11 +

5y )

dx

ТlФl2 +

/ \ да da m ^ „

+ (j + v)— -т-^Фп} +

dx dy

+ Ю

V = ю

^j_v2 da 2A + v) da

-----Ф 3----*----^Ф 4

E Sx E dy

j-v2f da 2^ + v)F da

Л

f2

E dx

+ ©{-DjO^ +

Fj----h

E 1 dy - —J

Эю^| ( 9ю^2'

+ ®j2 +

V

dy

dx I

Т^п -

/

/ \ da da m ^ „

_ ^ + ДГТ1Ф12>+

dx dy

+ ©

j _v 2 da 2A +v) da

E 9y E dx

Или U = U0 + Uj, где

j-v2 da 26 +v) da U0 — ю -----Fj------------Fo —

0 E j dx E 2 5y

V : J

^\-v2f da 2^ + v)F da^

V0 =ro

-Fj-----V

E j 5y

F2

E dx

124

РИ, 2000, № 2

f

U1 = Фц +ю{-01Ф11 +

V

Эю

dy

2

-V

Зю

Эх

2''

I Т1Ф12 +

/, \ Эю Эю m ^ „

+11 + vh“^-ТіФп},

ox dy

v1 - ф12 +ro{_D1®12 +

/ \ Эю Эю „ ,

- і1 + vj——— Т1Ф12}.

Эх dy

ґ fдюл2 ЛФ.

Эю

Эх

2''

І Т1Ф11 -

Для нахождения неопределенных компонент в приведенных структурных формулах разложим их

в ряд по некоторой полной системе функций {ф; }, например, по полиномам Чебышева, тригонометрическим, степенным или другим. Тогда искомые решения для вектора перемещений примут вид

N

w = Е aiBo№> raj,

i=1

N1

u = E сА(фь ra),

i—1

N2

v = E ciB^9i, ®), i=N1+1

где Bo, B1, B2 — операторы соответствующих струк-

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

3I^w)

3ai

= 0, i = 1,N;

= 0, i = 1iN2.

Эс; 2

(12)

Уравнения (12) можно представить в матричном

виде AX = B . Решая матричное уравнение, находим неопределенные коэффициенты ai, с; , а следовательно, и функции перемещений.

2. Описание алгоритма решения задач геометрически нелинейной теории пластин

Реализация изложенного выше алгоритма осуществлена в условиях эксплуатации программирующей системы POLE-RL [3]. Эта система содержит специальную библиотеку программ, которые осуществляют организацию данных сложной структуры, реализуют динамическое распределение памяти для них, а так же выполняют операции кортежной алгебры (арифметические, вычисления элементарных функций, R-функций, дифференциальных операторов и др.).

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

Рабочая программа состоит из блоков формирования матриц ФМ1, ФМ2, решений алгебраических систем РС1, РС2, организации итерационного процесса ITER. Исходными данными для ФМ1 являются нулевые коэффициенты С0= [c0i].

Блок формирования матриц ФМ1 осуществляет вызов модуля для интегрирования функций ащ и bii, который, в свою очередь, обращается к модулям вычисления структуры решения СТР U для функций u,v; преобразования геометрической информации для функции ю(х) в аналитическую; определения функциональных компонент краевых задач ФК. Работа блока формирования матриц завершается созданием матрицы A1=[a1ij] и столбца свободных членов B1=[b1i].

В блоке РС1 в качестве исходных данных используются матрицы A1 и B1. В результате решения алгебраической системы (она может быть не только системой Ритца) формируются коэффициенты С1=[с1і], которые являются решением плоской задачи теории упругости.

Следующие блоки ФМ2, РС2 работают аналогично, с той лишь разницей, что в них решается задача изгиба пластины с учетом напряжений. С1 являются исходными данными для блока ФМ2. Результатом работы этих блоков будут коэффициенты С2=[с2і].

Рис. 1. Блок-схема для решения задачи о закритическом поведении пластин

РИ, 2000, № 2

125

Для эффективного использования памяти вычислительной машины массивы A1, B1, A2, B2 уничтожаются после получения соответствующих коэффициентов С1 и С2. В свою очередь коэффициенты СО, С1 и С2 сохраняются программой до следующего шага итерации.

После решения второй задачи блок ITER вызывает для сравнения коэффициенты СО и С2. Сравнение осуществляется по формуле

Z (c0i _ c2i)2

|5Cj| =-i-----------------100% <є

Е c02i i

(13)

где |8C^ — относительная погрешность в процентах; є — достаточно малая положительная величина.

Если условие (13) не выполняется, то в блоке ITER выполняются следующие операции:

1) коэффициенты С2 записываются на место СО (прежние СО уничтожаются);

2) управление передается на начало ФМ1 (метка М). Весь процесс вычислений повторяется.

Если условие (13) выполняется, тогда управление передается в блок выдачи результатов ВР. И после этого программа завершает свою работу.

Приведем текст программы на языке RL для рассматриваемой задачи.

declare

bandy f1;bandx f2; sfpl,ql,rl;pol0 lm3,e1,e2; sf0 p01,p02;

omega

w1=f1;w2=f2;

w3=w1&w2;

function

w2=f2;w1=f1;

w4=w1&w2;

function

om=w4;p=rl;

u3=om*p;

function

u4=sum(1, 0 ,u3); function

om=w4;wx=dx(om);wy=dy(om); wxy=wx*wy; r=u4;

r3=d 1 (r,om);r4=t 1 (r,om);

r1=-(r3*r3+e1*r4*r4)/2-4.5*(1-e1*e1)*w2/(w1+w2);

r2=-r3*r4;

u0=om*(wx*r1-wy*r2);v0=om*(wy*r1+wx*r2);

p=pl;q=ql;

u1=(p+om*(lm3*wxy*t1(p,om)-d1(p,om))#

om*(1.0-lm3*wx*wx)*t1(q,om))+0.0;

v1=om*(lm3*wy*wy-1.0)*t1(p,om)#

q-om*(lm3*wxy*t1(q,om)+d1(q,om));

function

vr=sum(2,v0,v1);

ur=sum(2,u0,u1);

u=f1;

end

function

tr10=u; u31=u3; function

u5=sum(3,0,u3);

program

rko;

nel;

ama;

gauss(s1,fa,fb);

gil(s11,fbk);

sis(aa,bb);

ama;

gauss(s1,fa1,fb1);sis(aa,bb); йЙє(‘Проги6 в центре пластины’); tabp(ss1,fs2); time;

iter(s29,FITE);

FITE=0.001;

fs2(u5)=u5;

fa(v1)=u1(i,2)*u1(j,2)+v1(i,3)*v1(j,3)+e1*(u1(i,2)*

v1(j,3)+v1(i,3)*u1(j,2))+(1-e1)*

(u1(i,3)+v1(i,2))*(u1(j,3)+v1(j,2))/2;

fb(v1)=-(u0(2)*u1(i,2)+v0(3)*v1(i,3)+e1*

(u0(2)*v1(i,3)+v0(3)*u1(i,2))+(1-e1)*

(u0(3)+v0(2))*(u1(i,3)+v1(i,2))/2+(r(2)*r(2)*

(u1(i,2)+e1*v1(i,3))+r(3)*r(3)*(v1(i,3)+e1*u1(i,2))+

(1-e1)*r(2)*r(3)*(u1(i,3)+v1(i,2)))/2);

fbk(v1)=-((1-e1)*r2*(v1(i,1)*wx-u1(i,1)*wy)/2+

r1*(u1(i,1)*wx+v1(i, 1)*wy)+

(u1(i,1)*(wx*(r(2)*r(2)+e1*r(3)*r(3))+

(1-e1)*r(2)*r(3)*wy)+

v1(i,1)*((1-e1)*r(2)*r(3)*wx+

wy*(r(3)*r(3)+e1*r(2)*r(2))))/2);

aa=ai(1);bb=bi( 1);

126

РИ, 2000, № 2

p1=(ur(2)+e1*vr(3)+(r(2)*r(2)+e1*r(3)*r(3))/2)/ 3. Численные результаты

(1-e1*e1);

p3=(ur(3)+vr(2)+r(2)*r(3))/(2+2*e1);

p2=(vr(3)+e1*ur(2)+(r(3)*r(3)+e1*r(2)*r(2))/2)/

(1-e1*e1);

fa1(u31)=(u31(i,4)+u31(i,6))*(u31(j,4)+u31(j,6))-

(1-e1)*(u31(j,6)*u31(i,4)+u31(j,4)*u31(i,6)

-2.0*u31(i,5)*u31(j,5));

fb1(u31)=- e2*(p1*u31(i,2)*r(2)+

p3*(u31(i,2)*r(3)+r(2)*u31(i,3))+

p2*u31(i,3)*r(3))+

0.001*u31(i,1)*e2;

end

value

const=3,3,1,4,80, 1,1; tabl=nu,2,1,2,2,nv,1,2,2,2,nw,1,1,2,2,

0,3,1,0,0, 0,1,2,0,0, 0,3,1,0,0;

Пример 1. Исследуем закритическое поведение свободно опертой квадратной пластины, сжимаемой усилиями, направленными вдоль оси ОХ (рис. 2).

А

а

V

Y / к

О

^ Х

Рис. 2. Пластина под воздействием силы

На рис.3 результаты, полученные авторами, сравниваются с результатами работы [4]. Значение приложенной нагрузки вычисляется по формуле

\2

, а безразмерный параметр прогиба в

Pi0 =а-ЕІ h

центре пластины определяется так:

:: wo = w(0,0)/h .

nu=12; nv=12; nw=11; pl=1,pr; ql=2,pr; rl=3,pr; pr=am,bm,a,b; p01=1; p02=2; am=-.5; a=.5; bm=-.5; b=.5; f1=0,a; f2=0,b; s11=k1,a,0,a,b; k=3; s1=k,0,0,a,b;

----метод R- функций ...[4]

k1=5;

s29=1,3,2,2,0;

e1=0.3; e2=10.92;

ss1=0,0;

lm3=1.3;

end

В разделе declare описываются переменные, используемые в дальнейшем, например, полосы, полиномы, кортежи и т.д. Задание геометрической информации и структур решения происходит в function. Далее начинается выполнение тела программы. Считываются начальные коэффициенты rko. Интегрирование по области и по контуру производится с помощью процедур gauss и gil. Полученные алгебраические системы решаются в процедуре sis. Блок итераций заканчивается командой iter. В разделе value присваиваются значения переменным, а также используются стандартные идентификаторы const и tabl для задания информации о неопределенных компонентах (их количестве, типе симметрии), о матрицах (их структуре, типе, количестве), о максимальном порядке дифференциальных кортежей и т.д.

В табл. 1 изучена практическая сходимость значений прогиба в центре пластины при фиксированной нагрузке а =4,5 при изменении количества координатных функций Ni для u, v и N — для w.

Таблица 1

N 6 10 15 21

6 1,06 1,18 1,19 1,19

10 1,07 1,20 1,21 1,21

15 1,08 1,21 1,21 1,21

21 1,09 1,21 1,21 1,22

28 1,09 1,21 1,21 1,22

Проанализировав результаты табл. 1, можно сделать вывод, что, начиная с 10 координатных функций для u, v и w, стабилизация значений прогиба идет в третьем знаке.

Пример 2. Рассмотрим свободно опертую квадратную пластину с врезами на двух параллельных сторонах, сжатую усилиями, направленными вдоль оси ОХ (рис. 4).

Исследуем соотношение между прогибом и нагрузкой при различной глубине врезов: R/a = 0,1 (табл. 2); R/a = 0,01 (табл. 3).

РИ, 2000, № 2

127

У- a ^.

—=> p О JT*\ є с . P .

—э —э і і і і і і it 4 X е £

2R <—>

Рис. 4. Пластина с врезами

Таблица 2

а

6,1 6,3 6,5 7 7,5 8 8,5

wo 0 0,523 0,781 1,25 1,62 1,96 2,27

Таблица 3

а

4,9 5 5,5 6 6,5 7 7,5

wo 0 0,120 0,889 1,29 1,63 1,94 2,23

Из приведенных выше результатов следует, что при увеличении радиусов врезов пластина становится более устойчивой.

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

Литература: 1.ВольмирА.С. Гибкие пластинки и оболочки. М.: Гостехиздат, 1956. 420 с. 2.РваневВ.Л. Теория R-функций и некоторые ее приложения. К.: Наук. думка, 1982. 552 с. 3.Шевненко А.Н. Программирование структур решения краевых задач на языке RL-1. Харьков: Ин-т пробл. машиностр., 1979. 23 с. 4. Корнишин М.С., ФайзулинаМ.А. Геометрически нелинейные задачи изгиба и устойчивости пластин сложной формы // Исследования по теории пластин и оболочек. 1984. Вып. 17, ч. 1. С. 20-31.

Поступила в редколлегию 29.02.2000

Рецензент: д-р техн. наук, проф. Курпа Л.В.

Линник Анна Борисовна, аспирант кафедры прикладной математики ХГПУ. Адрес: Украина, 61142, Харьков, ул.Блюхера, 38а, кв.249, тел.40-09-41. E-mail: linik@kpi.kharkov.ua

УДК 519.713

О РЕШЕНИИ ЗАДАЧ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ В ЦЕЛЫХ ЧИСЛАХ

ЛАРИОНОВ Ю.И., ХАЖМУРАДОВ М.А.________

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

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

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

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

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

Число симплекс-итераций при решении задач ЛП зависит в основном от числа ограничений, а не от количества переменных. В грубом приближении число необходимых итераций лежит в пределах от 1,5 до 2m, где m — количество ограничений — уравнений в начальной симплекс-таблице.

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

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

128

РИ, 2000, № 2