Определение градиента целевой функции в задаче о минимуме концентрации напряжении на основе метода конечных элементов Текст научной статьи по специальности «Физика»

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

Том XXI 1990

№ 3

УДК 629.7.015.4.023

ОПРЕДЕЛЕНИЕ ГРАДИЕНТА ЦЕЛЕВОЙ ФУНКЦИИ В ЗАДАЧЕ О МИНИМУМЕ КОНЦЕНТРАЦИИ НАПРЯЖЕНИЙ НА ОСНОВЕ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ

В. И. Гришин, Ф. В. Рыбаков

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

Одним из наиболее эффективных методов обеспечения прочности элементов авиационных конструкций является поиск форм деталей, минимизирующих коэффициент концентрации напряжений ^(¿)=тахст/р, где а — максимальное главное напряжение

для плоских или интенсивность напряжений по Мизесу для осесимметричных элементов конструкций; Г — контур варьируемой поверхности; р — номинальное напряжение; < — вектор проектных переменных.

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

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

В данной работе для анализа чувствительности применяется конечно-элементный подход с использованием 8-узлового изопараметрического конечного элемента с квадратичной аппроксимацией поля перемещений. Вместо обычно используемого вычисления производных от координат путем простого численного дифференцирования предложен аналитический метод. Приводится исследование точности вычисления градиента целевой функции в задаче оптимизации формы тела при варьировании:

а) только граничных узлов конечно-элементной модели;

б) граничных узлов и непосредственно примыкающих к ним промежуточных узлов в граничном слое элементов;

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

1. Конечно-элементный подход к определению градиента целевой функции.

Рассмотрим систему уравнений равновесия конструкции

ки=и,

тде К — матрица жесткости конструкции; £/ — вектор узловых перемещений; Н — вектор внешних сил.

Дифференцируя по параметру оптимизации tj, получаем

(1.1)

Если варьируемая граница тела не нагружена внешними силами, то из (1.1) получаем

ди

ді

£;и)- <1Л>

Для каждого изопараметрического конечного элемента, используя стандартную методику определения матрицы жесткости [6], получаем для плоского элемента:

і і

дК1 ГС Г дВт дВ д | У П

А-= I —— | У | + В7 — | УI + Вт ОВ

і .) J і ді) 1 діі ді) ]

ді

-ї -і

или для осесимметричного элемента: і і

ді

1 Г С [дВ дВ а ІУ І

г2’ і і кг>в|'|,'+а’0^ІЛг+г,І0Я^Г+

-1 -1

+ ВТ IУ1 ~1 ¿5 Л) ,

™)

где В — матрица деформаций,

/5 — матрица упругости,

]У|—определитель Якобиана преобразования координат,

г — координата текущей точки элемента вдоль радиуса (для осесимметричного тела).

дВт дВ д |УI дг Производные ------ ; -— ; ——; -— легко могут быть выражены через проору д{] (И^ а2у

изводные от узловых координат. Их значения для двумерных тел приводятся в [7], а для осесимметричных получаются аналогично.

Определив все вышеизложенные величины и производные из уравнения (1.2), по-

ди

лучаем вектор -— и, далее, дифференцируя известные соотношения МКЭ, получаем д1)

производные от деформаций е и от напряжений а

дв дВ . ди1 да де.

— =----- и1 + В ■— , — = О — .

д1] дtj дЬ)

дв да

Через величины г— и — легко получить выражение для производных любой це-дЬ) дЬ]

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

да1 1 I дах дау ______________________________________________________________________________________

ді) 2 \д£у д^ V (ах + ау)3 — 4 ( ах ау — ^ху)

и производной от интенсивности напряжений 01 в случае осесимметричном задачи

дн ] Г ЪЕ д 1 — 2^ ' д 7

^7 = 2^ [2(1 + .) й, (°г £г + °г £г + £0 + 1гг) ~ 2(1+.) М) + °г + °в)1'

где

/ЪЕ 1 — 2ч

2 (1 -)- м) (аг Ег +Зг Е* + а8 Е6 + Тег) — 2 (1 + м) (аг + аг + ае) ’

Е — модуль упругости материала, V — коэффициент Пуассона.

Величина (Гг получена через энергию формоизменения, являющуюся разностью потенциальной энергии и энергии изменения объема. В данной работе в качестве целевой функции задачи оптимизации использовались главное напряжение О! для плоских задач и интенсивность напряжений Яг для осесимметричных задач. И в том и в другом случае величины Ст1 и а, нормировались на некоторое номинальное напряжение р. Особенности применения метода подконструкций к вычислению производных от напряжений излагаются в работе [1].

2. Метод варьирования координат. Пусть форма контура концентратора напряжений в некоторой местной системе координат (ОХ^ ) задана уравнением

ф (*М, 9 =0 ’ <2Л>

где хм = [х1м, х1) — координата точки,

¿ = (<1, ¿2. ••■> *п) — вектор конструктивных параметров, и пусть переход из местной системы координат к исходной осуществляется по формуле

х = Ахм + Ь ,

где х — координата точки в исходной системе координат, А — матрица преобразования системы координат, Ь — вектор начала местной системы координат в исходной системе координат.

Пусть / — единичный вектор, определяющий направление трансляции контура. В местной системе координат это будет вектор 1Ж=А~Ч. Рассмотрим некоторую точку контура с координатами в местной системе координат при значении вектора

параметров *= ¿л) и возьмем приращение параметра Ыр тогда точка

п л дХуд

хи перейдет в точку хм = Так как направление трансляции за-

дх1

дается вектором /м> то вектор —— параллелен вектору /м, т. е.

д(]

дхм

^М = М. (2.2>

Для определения коэффициента пропорциональности (1 рассмотрим уравнение контура при значении параметра = $ + Ы] :

Ф (+ 1и? ы*’ * ■" $ + М'.. *”) = ° ’

Разложим функцию Ф в ряд Тейлора, ограничившись лишь линейными членами

*(<<•)+(ч*.-57)

дФ I

Учитывая (2.1) и (2.2), получаем Л= — — /(у-Км'Мм).

I

дФ

,

т е' (Ч,ф-'-) "•

Или в исходной системе координат

дх

=-Л

дФ діу

дЬ " (Ч,Ф’Ч

/м .

(2.3)

Если в качестве вектора /м взять вектор у* нормальный к контуру, то получаем

дФ

Ж)

V*.ф-

(2.4)

Выражения (2.3) или (2.4) определяют производную от координат произвольной точки варьируемого контура. При варьировании внутренних узлов подконструкции необходимо также определить производные от координат этих узлов.

Пусть конечно-элементная модель подконструкции, описывающей область концентратора, является топологическим прямоугольником (см. рис. 1). Рассмотрим произвольный внутренний узел І, который лежит на линии 1—/г, где к — узел на варьируемой границе Г, 1 — узел на противоположной неподвижной границе подконструкции.

дхь

Производную от координат узла к — найдем по выражению (2.4). Представим вектор

01)

дхк дхк

^5—в виде = /а , где і — единичный вектор направления трансляции узла, ори-

01]'. ОС)

вотированный внутрь тела, а — скалярная величина, определяющая длину вектора про-

дхі

изводной и его направление по отношению к вектору ]. Тогда производную —— найдем по выражению

дх-,

ді)

хі-1 - ■*

І+1

(¿-1)

(£ - 1) .

(2.5)

Таким образом, формулы (2.4) — (2.5) определяют производные от координат любого узла подконструкции, конечно-элементная модель которой меняется при изменении параметров оптимизации.

В случае если варьировать только граничные и промежуточные узлы в слое элементов, расположенных вдоль границы, в качестве производной координат промежуточ-

200

р— л

" X

а)

Рис. 1

Рис. 2

ного узла (наприме, узла к—1 на рис. 1), возьмем половину производной от координат соседнего граничного узла, т. е. для узла й—1 (рис. 1) принимаем

__ дхи

2

3. Численный анализ влияния учета перемещений внутренних узлов на точность вычисления градиента. Для оценки точности в вычислении градиента проведен численный эксперимент на примере пластины с отверстием (рис. 2, а). Из условия симметрии рассматривалась четверть пластины. Расчеты проводились для двух конечно-элементных моделей (КЭМ) № 1 (рис. 2,6) со 150 степенями свободы и № 2 (рис. 2, в) с 366 степенями свободы и трех размеров отверстия диаметром а!=20 мм, 30 мм, 40мм.

Форма одной четверти отверстия аппроксимировалась кривой второго порядка, уравнение которой имеет вид:

у2 + Аху Сх2 + йх + Еу + 0 = 0

с фиксированными точками М и N (рис. 2, а) на концах контура, с горизонтальной касательной в точке М и вертикальной касательной в точке N. В качестве параметра оптимизации ^ выбрана ордината точки контура с абсциссой а( = 0,71а, поделенная на а, где а — радиус исходного отверстия. Рассматривались формы отверстия с параметром ¿т = 0,710-^0,860, с шагом Д£ = 0,025.

Градиент целевой функции 1ЯХ1Р’ где —-максимальное главное напря-

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

/ дРлт дР*т

1 обозначим соответствующие величины производных —— , —— , —— 11 а также

V дt дЬ д1

методом прямого дифференцирования, то есть по формуле

дРт

дt

т + 1 '

т—1

2М

Ошибки

= 2, 3,..., 6 в вычислении градиента определялись по от-

ношению к величине, полученной методом прямого дифференцирования, т. е., например:

дРя

=

т

дЬ

д_1т

дt

100%

Итоговые усредненные по пяти значениям параметра ¿т = 0,735; 0,760; 0,785; 0,810;

0,835 ошибки е“р, £®р, £®р, полученные на двух конечно-элементных моделях, для трех размеров диаметра приводятся в табл. 1. Из таблицы видно, что наибольшая точность достигается при учете перемещений внутренних узлов. В этом случае точность вычисления градиента в 5—10 раз выше, чем при учете перемещений только граничных узлов. В табл. 2 для каждого способа варьирования координат приводится среднее время процессора ЭВМ ЕС-1055, необходимого для вычисления значения целевой функции и градиента один раз (величины ^а, тб, тв). Из табл. 1 и 2 видно, что первый способ вычисления производных, т. е. с учетом перемещения только граничных узлов не пригоден, так как по сравнению со вторым способом при равных временных затратах дает в 2—3 раза большую ошибку. Что же касается оставшихся двух способов, то учет перемещений всех внутренних узлов подконструкций для более грубых моделей дает заметное увеличение точности (до 5 раз), при увеличении же густоты разбиения конечно-элементной модели точность вычисления градиента этими двумя способами сближается. Предпочтение следует отдать третьему способу, учитывающему перемещение всех внутренних узлов, так как его совокупные характеристики лучше, чем у первых.

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

Таблица 1

КЭМ № 1 К ЭМ № 2

с1 *ср 4> *ср *СР 4 £ср

20 6,08 3,55 0,74 2,88 1,33 0,36

30 6,44 2,92 0,64 1,49 0,41 0,40

' 40 4,48 1,14 0,29 1,03 0,48 0,55

Таблица 2

КЭМ № 1 К Э М № 2

а Та Т® Тв Та Т® тв

20 1'02" 59" Г 59" 2'49" 2'50" 5'49"

30 59" Г02" 1'57" 2'47" 2'47" 5'33"

40 59" Г 00" Г55" 2'50" 2'50" 5'39"

в миллиметрах. На трех участках заданы распределенные внешние нагрузки Ри р%= = 0,128ри р3=1,25р1 (см. рис. 3).

Рассматривались эллиптические формы перехода (участок МЫ, рис. 3) с фиксированной величиной одной полуоси а=3,5 мм и переменной другой Ь = 1а, где t — безразмерный параметр оптимизации. Точка М фиксирована, положение точки N меняет-

3/ **• «у

ся при изменении параметра t. Целевая функция /'=—5125- представляет собой отноше-

Р\

ние максимальной интенсивности напряжений ч шах на контуре МЫ к номинальному напряжению р4 (рис. 3). Рассматривались формы перехода, соответствующие значениям

параметра <=1,00-5-1,30 с шагом Д<=0,05. Ошибки в вычислении градиента Ет’ >

Рис. 3

£т' т=2. 3,..6 усреднялись по пяти значениям параметра /то = 1,05; 1,10; 1,15; 1,20; 1,25; 1,30. В итоге получены средние относительные ошибки е®р = 31,7%; =

=4,69%, £ср =3,15%. Наилучшие результаты дает третий способ. Соответствующее время процессора, требуемое для вычисления один раз целевой функции и градиента в этих случаях, составляет та = 9'56", тб = 9'43", т® = 14'34" .

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

В заключение отметим, что в приведенных выше примерах по сравнению с традиционными исходными радиусными отверстиями и переходом оптимизация формы концентратора напряжений позволяет снизить концентрацию напряжений на 21—28%.

ЛИТЕРАТУРА

1. Гришин В. И., Рыбаков Ф. В. Оптимизация формы элементов авиационных конструкций с концентраторами напряжений. — Ученые записки ЦАГИ, 1985, т. 16, № 6.

2. В otkin М. Е. Shape optimization of plape and shell structures.— AIiAA Paper, april 1981, 81-0553.

3. Ramakrishnan С. V. and Francavilla A. Structural shape optimization using penalty functions. — J. of Structural Mechanics, 1975, vol. 3, N 4.

4. H о u J. W., Chen J. L. and Sheen J. S. Computational method for optimization of structural shapes. — AIAA Journal, June 1986, vol. 24, N 6.

5. К r i s t e n s e n E. S., Madsen N. F. On the optimum in—plane loading cases. — Int. J. for Numerical Methods in Engineering, 1976, vol. 10, N5.

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

1975.

7. .Francavilla A., Ramakrishnan С. V., Z i е n k i е-w i с z О. С. Optimization of shape to minimize stress concentration. — J. Strain Analysis, 1975, vol. 10, N 2.

Рукопись поступила 13/V 1988 г.