CC BY
51
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
Том XXIII 1992 № 3
УДК 629.7.015.4.023. + 539.4.013.3 ■
ОПТИМИЗАЦИЯ ФОРМЫ ПРОСТРАНСТВЕННЫХ ЭЛЕМЕНТОВ КОНСТРУКЦИЙ С КОНЦЕНТРАТОРАМИ НАПРЯЖЕНИЙ
Ф. В. Рыбаков
Решена задача определения формы нерегулярных зон пространственных элементов конструкций, минимизирующей концентрацию напряжений, в упругой постановке. Методика решения задачи основана на применении МКЭ в перемещениях. Задача свободной минимизации решается градиентными методами. Дан анализ эффективности применения метода переменной метрики Давидона—Флетчера—Пауэлла в сравнении с методом наискорейшего спуска. Предложен метод автоматического перестроения конечно-элементной модели в процессе оптимизации. Приведена оптимальная форма подкрепления пластины с отверстием, позволяющая значительно снизить концентрацию напряжений.
Наличие отверстий, вырезов, выступов, переходов и других геометрических нерегулярностей создает сложное распределение напряжений в конструкции. При этом возникают зоны местного возрастания напряжений, которые, как правило, ограничивают прочность конструкции.
Выработанные на протяжении десятков лет формы различных конструктивных элементов достаточно эффективны, и дальнейшее их улучшение возможно при точном анализе напряженно-деформированного состояния, например с помощью метода конечных элементов (МКЭ), и применении численных методов оптимизации. Несмотря на то что как те, так и другие достаточно хорошо изучены и описаны в целом (ряде монографий, соединение их в единый аппарат для решения задачи оптимизации формы требует, во-первых, создания достаточно эффективных быстродействующих алгоритмов; во-вторых, обеспечения высокой точности вычисления целевой функции и ее производных, что является необходимым условием сходимости алгоритма минимизации.
В данной статье применительно к объемным задачам развивается методика, предложенная ранее автором и описанная в [1], с использованием метода переменной метрики. Предлагается эффективный метод определения градиента целевой функции. Приводятся результаты решения задачи об оптимальной форме подкрепления отверстия в пластине в объемной постановке.
1. Постановка задачи оптимизации и метод ее решения. В качестве целевой функции в задаче оптимизации формы конструкции примем коэффициент концентрации напряжений, т. е. функцию
где а, — интенсивность напряжений по Мизесу в точке области /?, а — номинальное напряжение; 1 — вектор конструктивных параметров, определяющий форму границы области /?.
На область допустимых изменений конструктивных параметров накладываются линейные ограничения типа неравенств
*„(1X0, т=» 1,2,....М. (2)
Как правило, данные ограничения не являются активными, минимум функции находится внутри области допустимых значений достаточно далеко от ее границ. Ограничения необходимы для контроля за процессом оптимизации, чтобы не допустить образования физически невозможных форм. Поэтому для учета ограничений наиболее удобным является метод штрафных функций.
Введем новую целевую функцию
2(1) = Ф(1) + 0[*,(ВД,
где О —некоторая функция штрафа. Подробно ознакомиться с реализацией метода штрафных функций Читатель может в работе [1]. Здесь же рассмотрим задачу минимизации функции 1 без ограничений методом наискорейшего спуска и методом Давидона—Флетчера—Пауэлла [2].
Пусть при переходе от Л-й итерации к Л+1-й вектор параметров вычисляется по формуле
1*+1 = Ь — Л,* Н* VZ (и), (3)
где матрица Н*+1 вычисляется по формуле
мкм1 н>д(уг>)А(уг,)гн,
мткь(чгк) д(уг*)Чд(у2*)
и и I 4<~к) /лл
“*+1 — "* "I а*га/^г~ ГГЗГГТгТГТТЗГП ’ ' '
'VZ(tk) — градиент целевой функции на к-м шаге, Д^. и А (^Я*)— разность величин на &-м и £ — 1-м шаге.
Параметр шага А.* вычисляется из условия
г(\к- чЩк))
*
с помощью аппроксимации функции г (Ь — Я.Н* ЧЪь) полиномом третьей степени по к.
Если в формуле (3) в качестве матрицы Н* принять единичную матрицу, то получается метод наискорейшего спуска. Изменение матрицы Н* в процессе оптимизации по формуле (4) соответствует методу переменной метрики. При многократном преобразовании.матрицы по формуле (4) возможно накопление ошибки. Поэтому для обеспечения сходимости процедуры оптимизации нужно периодически восстанавливать первоначальное значение матрицы. В данной работе восстановление осуществлялось на каждом п-м шаге, где п — количество оптимизируемых параметров.
2. Перестроение конечно-элементной модели при изменении параметров оптимизации. При решении задачи оптимизации формы конструкции возникает необходимость изменения конечно-элейентной модели при различных значениях проектных параметров. Это можно осуществить с помощью специальных программ автоматического генерирования модели. Однако если варьируемая поверхность представима в виде гладкой функции, то задачу можно существенно упростить с помощью процедуры перестроения исходной конечно-элементной сети.
Пусть контур концентратора задан в некоторой местной системе координат в неявной форме уравнением вида
ф (х, I) = О,
где х = (х1, х2, х3) — координаты точки в местной системе координат. Переход из местной системы координат к исходной осуществляется по формуле:
X = А х + Ь, (6)
где X — координата точки в исходной системе координат; А — матрица преобразований координат; Ь — вектор начала местной системы координат в исходной системе координат.
Пусть Ь — единичный вектор, определяющий направление трансляции контура. В местной системе координат это будет вектор I = А-1 Ь. Рассмотрим некоторую точку контура с координатами хо в местной системе координат при значении вектора параметров = (*®, .... *2). Возьмем приращение па-
раметра Л/,-, тогда точка хо перейдет в точку
«
х, = Хо 4- (7)
где с1 — неизвестная скалярная величина, определяющая расстояние между точками X) й Хо.
Точка XI лежит на новой границе тела, соответствующей значению параметра 1 = 1° + т. е. координаты этой точки удовлетворяют уравнению
ф(хо + <*1, 1° + Д1) = 0.
Линеаризируя функцию ф, получим
ф(х0, 1) + {'Фх<рУШ = 0.
Учитывая, что ф(х„, /) Ф 0, получаем
=тг£' <8>
(^ф)1
где вектор(у,ф)'=(-£,
Формулы (7) и (8) определяют приближенное значение новых координат произвольной точки границы. Повторяя в цикле указанную процедуру, можно получить новые координаты с заданной точностью. Как показывает опыт, формулы (7) и (8) очень хорошо определяют новые значения координат и обычно достаточно трех—пяти итераций, чтобы получить новые координаты с точностью до седьмого знака.
В тех случаях, когда изменяемую границу нельзя описать гладкой функцией, задаваемой выражением {5), необходимо создавать автоматические генераторы сетей.
3. Вычисление градиента целевой функции. Для пространственных конструкций в соответствии с критерием Мизеса удобно воспользоваться интенсивностью напряжений в виде
~ 2(1+V) (°*е* "Ь ауеу ~ь~ агег + УхуХху + + УгхХгх) — 2'(1 + у) (°х °У '
являющейся функцией энергии формоизменения. С помощью закона Гука для линейно-упругого тела эта формула может быть записана в более простом виде:
4 [(о*- ст*)2 + (°# - аУ + (»* — Ох)2 + 6 (т2ху + х% + т2,)]. (9)
Дифференцируя соотношение (9), легко получить выражение для производной doi/dtj через производные от напряжений и деформаций. Дифференцируя известные из МКЭ связи между напряжениями и деформациями, деформациями и перемещениями, а также основную систему уравнений МКЭ в перемещениях, получаем все необходимые соотношения для вычисления производной doi/dtj. Эти соотношения аналогичны исходным формулам, по которым определяется сама величина, с той лишь разницей, что в них входят производные от координат узлов конечно-элементной модели. Вывод этих соотношений имеется в работах [3, 4]. При вычислении градиента целевой функции по этой методике для заданной конечно-элементной модели и машинной точности точность вычисления градиента определяется точностью производных от узловых координат. Методика определения производных от координат узлов для форм варьируемых границ, задаваемых уравнением (5), подробно изложена в работе [4]. В случае если варьируемую границу нельзя описать в виде гладкой функции, то для вычисления производных от координат узлов дается малое приращение параметра Atj, перестраивается конечноэлементная модель и производные вычисляются через отношение приращений АХ/А tj. В этом случае необходимо исследовать влияние величины приращения параметра A tj на точность вычисления производных.
4. Сравнение двух методов оптимизации. Рассмотрим круглый стержень со ступенчато изменяющимся сечением (рис. 1). Относительные размеры следующие: D/d = 2, l/d = 2, L/d — Ъ. В качестве параметров оптимизации будем рассматривать полуоси эллипса, в форме которого выполнен галтельный переход, поделенные на диаметр меньшего сечения, т. е. величины a/d и b/d (см. рис. 1). За исходную принята форма радиусного перехода с радиусом 0,2 d. Ей соответствует коэффициент концентрации напряжений, равный 1,48 (точка А). Коэффициент концентрации вычисляется по формуле К = шах о,/а, где о, — интенсивность напряжений по Мизесу, о — номинальное напряжение на участке малого диаметра. Конечно-элементная модель состояла из 235 квадратичных 8-узловых конечных элементов, что соответствует приблизительно 1600 степеням свободы. Оптимизация проводилась как методом Давидона—Флетчера—Пауэлла, так и методом наискорейшего спуска. Результаты расчета приведены на рис. 1. Ломаная АВ соответствует трем итерациям метода переменной метрики, ломаная АВ' — десяти итерациям метода наиско-рейшего спуска. Из рисунка видно, что в методе наискорейшего спуска идет колебательное движение с одного «берега оврага» на другой. Метод перемен-
/
ной метрики обеспечивает гораздо более высокую скорость сходимости. Для наглядности на рисунке приведены линии уровня целевой функции, т. е. max Oi/ст. Они получены путем оптимизации параметра эллиптической формы b/d при различных фиксированных значениях параметра a/d.
5. Исследование точности вычисления.градиента целевой функции. Рассмотрим пластину с отверстием, подкрепленным кольцевой симметричной окантовкой при одноосном растяжении (рис. 2). Относительные размеры следующие: отношение диаметра отверстия к ширине пластины d/W —0,2, отношение диаметра отверстия к толщине пластины d/т = 10, отношение радиуса сопряжения окантовки с пластиной (окрестность точки В на рис. 2) к толщине пластины г/х — 0,25. Объем материала окантовки Va равен объему отверстия V0. В качестве параметра оптимизации возьмем отношение внешнего диаметра окантовки к диаметру отверстия D/d. Высота окантовки для данного значения D/d при фиксированном объеме материала определяется по формуле
{D/df + VJV. - I п/% =-----------------.
(D/df-l
В данном примере варьируемая граница тела не является гладкой, поэтому для перестроения конечно-элементной модели в процессе оптимизации используется программа-генератор сети и производные от координат вычисляются через отношение приращений.
Исследуем на этом примере влияние величины приращения Д< на точность вычисления градиента целевой функции VZ. Определим относительную ошибку (в процентах) в вычислении градиента по формуле
На рис. 2 приведена полученная в результате расчетов зависимость ошибки е от величины приращения Л/ при значении параметра I = 1,2. Наилучшая точность достигается, если приращение А/ лежит в интервале значений от 10_3 до 10-4. Ошибка в вычислении градиента в этом случае составляет порядка 5%.
Рис. 3
Конечно-элементная модель, используемая в расчетах, состояла из 20-узло-вых квадратичных изопараметрических объемных элементов и имела около 5,3 тыс. степеней свободы.
6. Оптимизация формы подкрепления пластины с отверстием. Проведены расчеты по оптимизации параметров кольцевого подкрепления в пластине с отверстием (относительные размеры описаны в п. 5) для пяти значений отношения объема накладки к объему отверстия VJV0 = 1; 2; 3; 4; 6. На рис. 3 приведены полученные в процессе оптимизации максимальные значения а,/о в зависимости от отношения внешнего диаметра окантовки к диаметру отверстия D/d. Коэффициент концентрации напряжений для неподкрепленного отверстия в этом случае равен 3,14 (см. [5]). Максимум интенсивности напряжений реализуется всегда на внутренней поверхности отверстия в точке А (см. рис. 2). Из рис. 3 видно, что при фиксированном объеме окантовки минимум коэффициента концентрации напряжений, который можно достичь, оптимизируя соотношение внешнего диаметра окантовки и ее высоты, изменяется от 2,35 при VJV0 = 1 до 1,51 при VJV0 = 6. На рис. 3 также показаны оптимальные сечения окантовки при VJV0= 1; 2; 3; 4; 6. Параметры этих сечений и соответствующие им максимальные значения о,/о приведены в табл. 1.
Таблица!
Параметры оптимального подкрепления
VJV0 D/d *А . Oi/O
1 1,25, 2,79 2,35
2 1,36 3,36 2,04
3 1,43 3,38 1,83
4 1,50 4,20 1,71
6 1,56 5,0 1,51
Подкрепление отверстия в виде наплыва. Рассмотрим ту же пластину, что и в предыдущем пункте. Для снижения концентрации напряжений введем симметричное подкрепление в виде утолщения в области максимальных напряжений (рис. 4). Из условия симметрии на рисунке показана 1/8 часть пластины в окрестности отверстия. Размеры утолщения определяются тремя параметрами t\, ti и <з (см. рис. 4) Параметр t\ определяет максимальный размер наплыва в поперечном сечении, равный t\d/2\ t2 — угол начала наплыва на контуре отверстия, равный hn/2, ta— толщину наплыва, равную /зт/2. Форма поверхности BCD наплыва определяется следующим образом. Поверхность наплыва гладкая, по дуге CD плавно переходит в плоскость внешней поверхности пластины. В вершине В касательная плоскость параллельна пластине. Кривая ВС, лежащая на цилиндрической поверхности отверстия, задана в цилиндрических координатах уравнением
2* = /3T/2sin2--~^”/2 , (10)
где а —угол, изменяющийся в интервале от hn/2 до я/2.
Дуга CD, лежащая в плоскости внешней поверхности пластины, задана уравнением
г* = d/2 (l + f 1 Sin ° x j”/2 ). (11)
Кривая BD и другие линии пересечения внешней поверхности наплыва с произвольной поверхностью ф= const задаются полиномом третьей степени
z** = 2* (а) + А (г - d/2У + В (г- d/2)3, (12)
где коэффициенты А и В определяются из граничных условий и равны:
Л— 3 z*(a)~T/2 • В — 2 z*(g)~T/2
(/■*(<*)-d/if ' (г* (а) -d/2f
Здесь и далее запись z* (а) и г* (а) означает, что соответствующие величины являются функциями угла а.
Объем материала, необходимого для формирования такого наплыва:
л/2 г* г**
Vu = ^ rdrdadz = ^ da ^ rdr ^ dz,
V a,, d/2 т/2
где ao = t2n/2 при <2^0 И ССо = 0 при /2<0.
После подстановки соотношений (10) — (12) и выполнения интегрирования получаем:
при /2^0
Кн = ^/2(й/2)%( 1 -12) + 4-),
при /2 < О
(13)
+«■££(• ■
(14)
Используемая при расчетах конечно-элементная модель (рис. 4) состояла из 20-узловых объемных конечных элементов с квадратичной аппроксимацией поля перемещений и имела около двух тысяч степеней свободы.
Коэффициент концентрации напряжений, определенный на данной модели, при отсутствии наплыва равен 3,30, что на 4,8% выше точного значения.
Оптимизация проводилась по всем трем параметрам при ограничениях t\ ^ 0, h 1, *з ^ — 1 ■ Начальные значения параметров t\ = 0,7, U — 0, t3 — 3,0. Ограничений на объем материала не накладывалось. В результате расчетов получено, что оптимальной форме наплыва соответствуют параметры t\ = 1,15, /2 = 0,06, /з=3,09. При этом минимум интенсивности напряжений по Мизесу равен 1,62, объем материала наплыва Уя= 1,84 V0. Как видим, при решении задачи минимизации концентрации напряжений в пластине с отверстием при одноосном растяжении в объемной постановке процедура оптимизации не приводит к бесконечному росту толщины. Существует оптимальное значение объема подкрепляющего наплыва, равное 1,84 объема отверстия. Это происходит вследствие того, что с увеличением объема материала наплыва увеличит вается кривизна поверхности в области сопряжения наплыва с плоскостью пластины и в этой зоне возникает локальная концентрация напряжений.
Оптимизация формы утолщения при фиксированном объеме материала.
Рассмотрим предыдущую задачу в предположении,, что объем материала в процессе оптимизации не меняется, т. е. наложим ограничение
VJVQ= const (15)
и исследуем зависимость оптимальной формы от отношения VJVQ. Пусть форма наплыва определяется соотношениями (10) — (12). Условие (15) позволяет выразить один из параметров через два остальные. Используя соотношения (13), (14) и то, что V0 = nd2/32, из (15) получаем:
при *2 ^ 0
/з=Уи/У0я/[4/|(1-/2)(-^Л+4-)].
, при /2 < 0
{-я^, [т» + (1 - ^ - 4<1 - <>) +
, *2Я/2 /, 1 2 12Л/2\\ V
+ с05—('--со5 —)} >'
Далее оптимизация проводится по двум параметрам и /г при 'значениях отношения Уя/У0, изменяющихся от 0,4 до 4,0. Начальные значения параметров
Таблица 2
Оптимальные параметры наплыва прн фиксированном объеме материала
VJV0 /Со *ор. h 1 + /3
0,4 2,31 1,99 0,72 0,43 3,23
0,8 1,94 1,84 0,69 0,22 4,40
1,2 1,93 1,66 1,11 0,16 3,94
1,6 2,02 1,59 1,05 0,06 4,70
2,0 2,09 1,61 0,97 —0,06 4,87
2,8 2,19 1,62 0,87 -0,11 7,16
4,0 2,27 1,70 0,79 —0,12 10,82
t\ и /2 соответственно равны 0,7 и 0, В табл. 2 приведены основные результаты расчета, т. ё. оптимальные значения параметров /2, 1 + Н, соответствующий им коэффициент концентрации напряжений Кор{ и коэффициент концентрации напряжений при исходных значениях параметров Ко-
Из таблицы видно, что при увеличении объема наплыва Ун > 2У0 минимум коэффициента концентрации напряжений начинает расти, что подтверждает верность описанных выше расчетов при оптимизации по трем параметрам.
На рис. 5 для сравнения показаны оптимальное кольцевое подкрепление при Ун/К0=6 и оптимальный наплыв при Ун/У0 = 1,84. Соответствующие йм коэффициенты концентрации напряжений, равные 1,51 и 1,62, различаются лишь на 3% значения коэффициента концентрации напряжений для непод-крепленного отверстия. В то же время затраты материала на наплыв в три раза меньше, чем на кольцевое подкрепление.
ЛИТЕРАТУРА
1. Гришин В. И.. Рыбаков Ф. В. Оптимизация формы элементов авиационных конструкции с концентраторами напряжений//Ученые записки ЦАГИ. — 1985. Т. 16, № 6
2. Пол а к Э. Численные методы оптимизации. — М.: Мир, 1974.
3. Francavill a A., RamakrishnanC. V., Zienkie-w icz О. С. Optimization of shape to minimize stress concentration // J: Strain Analysis. — 1975. Vol. 10, N 2.
4. Гришин В. И., Рыбаков. Ф. В. Определение градиента целевой функции в задаче о минимуме концентрации напряжений на основе метода конечных элементов // Ученые записки ЦАГИ. — 1990. Т. 21, № 3.
5. Петерсон Р. Е. Коэффициенты концентрации напряжений. — М.: Мир, 1977.
Рукопись поступила 29/VIII 1990 г.
