ТЕОРИЯ И АНАЛИЗ МЕТОДОВ ТОПОЛОГИЧЕСКОЙ ОПТИМИЗАЦИИ Текст научной статьи по специальности «Физика»

УДК 624.04

DOI: 10.18698/2308-6033-2023-4-2264

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

© П. А. Косых1, А.В. Азаров1,2

1 МГТУ им. Н. Э. Баумана, Москва, 105005, Россия 2 АО Центральный научно-исследовательский институт специального машиностроения, Хотьково, 141371, Россия

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

Ключевые слова: топологическая оптимизация, оптимальное проектирование, тело переменной плотности, SIMP-метод, BESO-метод, метод конечных элементов, аддитивные технологии

Введение. В последнее время наблюдается стремительное развитие аддитивных технологий. Современные методы 3D-печати позволяют получать изделия уже не только из обычного пластика, но и из металлов и композитов. Существует много разных способов печати изделий: из пластиков — методами послойного наплавления полимерной нити FDM (Fused Deposition Modeling) и стереолитографии SLA (stereolithography), из металлических порошков — методами селективного лазерного спекания SLS (Selective Laser Sintering) и сплавления SLM (Selective Laser Melting), из композитов — методами печати с использованием непрерывного армирующего волокна (CFF — Continuous Fiber Fabrication, CFC — Continuous Fiber Coextrusion).

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

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

Один из путей решения задачи оптимального проектирования — топологическая оптимизация (ТО). В настоящее время появилось достаточно много методов ТО, из которых наиболее широкое распространение получили SIMP-метод, BESO-метод, Level-set-методы, пропорциональная ТО [3-9]. В разных изданиях представлены способы оптимизации структур, напечатанных на 3Б-принтере и армированных непрерывным волокном. При таком подходе оптимизируется не только топология изделия, но и ориентация волокон [10-13].

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

Постановка задачи ТО. Рассмотрим постановку задачи ТО для упругого тела. Пусть в рабочей области D, для которой определены места креплений и приложения нагрузки, ведется поиск оптимальной по некоторому критерию конструкции Q с границей Г. Данная конструкция находится под воздействием объемных сил f(х) и поверхностных сил t(x), где х — координата точки. Под действием приложенных сил возникает поле перемещений и(х). Решение ищется в виде распределения материала в области D.

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

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

Выражение для податливости l(u), которая определяется через работу внешних сил над упругим телом, можно представить следующим образом:

l(u) = J fudQ + J tudr, (1)

Q

а выражение для потенциальной энергии деформации а(и, и) будет иметь вид

a(u, и) = J Ет (х)г tJ (и)г и (u)d Q (2)

du, du.

Q

f д,. Л

где Eijki (x) — тензор упругих постоянных; Sj (u) = —

j

dx: dx, V J i J

линейные деформации тела.

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

С учетом всех введенных обозначений и ограничений постановка рассматриваемой задачи ТО будет иметь следующий вид:

min l (u );

a(u, v) = l (v); (3)

JdQ<V,

Q

где v — кинематически допустимые вариации поля перемещений; V — максимальный объем, допустимый для конструкции.

Для решения сформулированной задачи традиционно используется метод конечных элементов (МКЭ). В рамках представленной работы решение задачи ТО выполняется в пакете MATLAB. Для того чтобы избежать трудностей, связанных с построением конечно-элементной сетки (КЭС), нумерацией узлов и элементов, а также с составлением глобальной матрицы жесткости, во всех рассматриваемых задачах рабочая область D имеет простую прямоугольную форму, приняты простые граничные условия (ГУ), а сама КЭС состоит из квадратных четырехузловых элементов.

Приведем формулировку задачи ТО, используя принятые в КЭ-анализе обозначения:

min fти;

Ku = f; (4)

ZV < V, Ve e Q,

где f — вектор приложенных сил; u — вектор перемещений; K — глобальная матрица жесткости; Ve — объем элемента e.

Г

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

В текущей работе рассмотрены два варианта оптимизируемых конструкций, которые стали классическими для задач ТО. В первой задаче — консольно закрепленная балка, нагруженная на свободном конце, во второй задаче — балка, которая опирается концами на шарниры и нагружена по центру. Во второй задаче часто выполняется оптимизация только половины балки вследствие ее симметричности, поэтому граничными условиями для такой балки с одной стороны станет шарнир, а с другой — запрет на горизонтальные перемещения. Первую из них будем называть задачей о консольной балке, а вторую — о шарнирной. Схемы этих задач приведены на рис. 1.

81МР-метод. Этот метод является одним из наиболее распространенных подходов к решению задачи ТО. Его развивали авторы М.П. Бендсё, Н. Кикучи, О. Зигмунд, а более подробное описание 81МР-метода приведено в монографии [3]. Широкое распространение получил код на языке МЛТЬЛБ, реализующий данный метод с использованием всего 99 строк [14].

Определим для материала два состояния — наличие или отсутствие. Тогда проектную область В можно разбить на две части: ^ — область, где материал есть, ЭЮ — область, где его нет. При таком варианте задачи получается черно-белое решение, т. е. при этом область ^ можно представить черным цветом, а область ЭЮ — белым. Однако чаще всего дискретное поведение материала будет неудобным для нахождения решения, поэтому стоит ввести промежуточные состояния материала, которые можно условно назвать серыми. При таких состояниях тензор жесткости Ещ (х) содержит лишь часть

свойств от полного тензора Еук1, а решение имеет наряду с черными

и белыми областями еще и серые. Для описания промежуточных состояний материала следует ввести непрерывную величину, характеризующую наличие и свойства материала в точке. Эту величину принято называть плотностью. Примеры полученных 81МР-методом решений с различной степенью серости приведены на рис. 2.

Рис. 1. Схемы задач топологической оптимизации о консольной балке (а) и о шарнирной балке (б)

Рис. 2. Демонстрация полученных SIMP-методом решений с различной степенью серости

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

№(x) = Р(x)pew , p >1

E

О < р(x) < 1

x е Q,

(5)

где p(x) — плотность материала, изменяющаяся в пределах от О (полное отсутствие) до 1 (полное наличие); p — коэффициент штрафа.

Как видно из представленной модели, при значениях p > 1 вклад промежуточного значения плотности в жесткость элемента уменьшается. При таком штрафовании промежуточных значений их невыгодно использовать.

Описанная в (5) модель материала, которая служит основой SIMP-метода, объясняет его название — Solid Isotropic Material with Penalization (пенализация для твердого изотропного тела). Используя ее, задачу ТО можно свести к задаче поиска оптимального распределения плотностей p(x). Представим эту задачу, дополнив постановку (3) моделью (5):

min l (u);

a(u, v) = l (v);

EljM (x) = p( x) РЦЫ.

Jp( x)d Q

< V,

p >1; о <Р,

(6)

<p( x) < 1,

где рш1п — нижняя граница плотности, которая необходима для того, чтобы обеспечить разрешимость задачи равновесия.

Приняв, что в силу условия равновесия а(и, и) = 1(и) и что объем оптимизированной конструкции равен максимально допустимому, введем множитель Лагранжа Л, а также функцию Ь, минимум которой соответствует минимальной податливости:

Г \

L = a(u, u) + Л

Jp( x)d П- V

VQ

(7)

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

в j (и )в kl (u) = Л. (8)

ф

Для промежуточных значений плотности (pm;n < p(x) < 1) выражение (8) можно записать так:

pp(x)p-1 j у (и)в к! (и) = Л. (9)

Отсюда видно, что левая часть (9) равна значению удельной энергии деформации, домноженной на коэффициент pp(x)p-1. Если предположить, что высокая энергия деформации элемента соответствует его низкой жесткости, схему обновления плотностей на i-й итерации можно представить следующим образом:

Р/+1 =

max {(1 - С)Р/, Pmin}, если PiB? < max К1 -С) Р/, Pmin}; min{(1+ С)Р/,1}, если min{(1+ С)Р/,1}<рВ; (10)

pB — в остальных случаях,

где Б1 определяется по формуле

Б, х)Р"1 Е»к1 г у (и)г « (и). (11)

Здесь Б, используется как коэффициент для оценки жесткости элемента, а Л, устанавливает уровень энергии деформации на текущей итерации. Ясно, что оптимальное значение Б, = 1 для каждого элемента. Согласно схеме (10), плотность элемента увеличивается при недостаточной жесткости элемента (при Б, > 1) и уменьшается при завышенной жесткости. При этом не допускается выход плотности из установленных границ. Параметр £ контролирует максимальное приращение плотности, которое может произойти за одну итерацию,

а параметр п определяет степень влияния В, на плотность. Стандартные значения для £ и п равны 0,5 и 0,2 соответственно.

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

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

Зависимость решения от конечно-элементной сетки. Замечено, что при использовании разного разбиения проектной области на конечные элементы решения могут различаться. Вопреки ожиданиям того, что более мелкая КЭС уточнит решение и сделает границы более четкими, в действительности может приводить к кардинально другому решению. При уменьшении размера конечного элемента наблюдается решение с большим количеством мелких элементов, более «тонкой» структурой и большим количеством отверстий. Такой эффект обусловлен тем, что при неизменном объеме увеличение количества отверстий в конструкции и измельчение ее структурных элементов повышают эффективность этой конструкции (рис. 4).

Для того чтобы решить отмеченные проблемы, можно применить несколько подходов. Среди основных — ограничение длины периметра и усреднение значения плотности элемента или его чувствительности с соседними элементами. Рассмотрим метод усреднения чувствитель-

б

Рис. 3. Эффект шахматной доски в задаче о консольной (а) и шарнирной (б) балке

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

д!

дРе

= -РРр 1и тК еие,

(12)

где ие — вектор перемещений элемента; Ке — матрица жесткости элемента.

Рис. 4. Изменение решения при измельчении сетки консольной балки размером 60x30 мм с количеством элементов 60x30 (а), 120x60 (б), 240x120 (в) и шарнирной балки 80x20 мм с количеством элементов 80x20 (г), 160x40 (д), 320x80 (е)

Усреднение чувствительности элемента к с окружающими его элементами е происходит по следующей схеме:

д!

дРк

1 N

-N-2 Нере

Рк 2 не е=1

е=1

д!

дре

(13)

где N — общее количество элементов в сетке; Не — независимый от сетки весовой коэффициент,

Не = Гшп -^(к,е), ё^к,е) < гШ1П, к,е = 1,...,N.

(14)

Здесь ё1в1:(к,е) определяется как расстояние между центром текущего элемента к и центром соседнего элемента е, с которым происходит усреднение.

Согласно (14), весовой коэффициент Не линейно убывает при удалении от центра текущего элемента и равен нулю за пределами окружности радиуса гт1п.

Влияние параметров 81МР-метода на решение задачи ТО. Для

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

В качестве тестовых рассмотрим задачи о консольной и шарнирной балках толщиной 5 мм. Консольная балка имеет длину 60 мм и высоту 30 мм, шарнирная — 80 мм и 30 мм соответственно. Для консольной балки размер КЭС составляет 180x90 конечных элементов, для шарнирной — 240x90. Обе балки изготовлены из материала, имеющего модуль упругости Е = 200 ГПа, коэффициент Пуассона V = 0,3. В тестовых задачах варьировались два параметра. В то время, когда изменялся один из них, другой сохранял стандартное значение. При стандартной настройке радиус фильтрации гШп = 1,5 мм, коэффициент штрафа р = 3, параметры £ = 0,2 и п = 0,5. Во всех вариантах расчета было установлено ограничение объема, составляющее половину от начального. Для каждого варианта было принято максимальное количество итераций, равное 100.

Результаты расчетов для консольной и шарнирной балок при варьировании двух параметров представлены в табл. 1, 2. В них каждому параметру отведено по три строки: в первой указано значение параметра, во второй — значение целевой функции I, в третьей — количество проведенных итераций N. Получившиеся в результате расчетов конструкции приведены на рис. 5 и 6.

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

Таблица 1

Результаты расчетов для консольной балки при варьировании параметров гт1„ и р

Варьируемый параметр Значение параметра / Значение целевой функции 1, Н мм / Количество итераций N

Радиус фильтрации Гшт, мм 0,5 0,7 1,0 1,5 2,0 5,0

63,4 64,5 65,5 67,4 69,4 98,9

39 32 26 27 34 32

Коэффициент штрафа р 1 2 3 5 8 15

56,5 69,5 67,4 69,6 70,7 80,5

16 47 27 26 44 63

Таблица 2

Результаты расчетов для шарнирной балки при варьировании параметров гт1„ и р

Варьируемый параметр Значение параметра / Значение целевой функции 1, Н мм / Количество итераций N

Радиус фильтрации Гтш, мм 0,5 0,7 1,0 1,5 2,0 5,0

145,4 146,5 150,2 153,8 156,8 192,1

51 50 44 42 34 45

Коэффициент штрафа р 1 2 3 5 8 15

131,0 156,5 153,8 156,3 164,6 175,2

17 51 42 46 100 100

б

Рис. 5. Варианты конструкции при варьировании параметров гтш, мм (а)

и р (б) (см. табл. 1)

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

б

Рис. 6. Варианты конструкции при варьировании параметров rmin, мм (а)

и p (б) (см. табл. 2)

BESO-метод. Для решения задачи ТО также применяется относительно простой BESO-метод (Bi-directional Evolutionary Structural Optimization — двунаправленная эволюционная структурная оптимизация), подробное описание, развитие и применение которого приведено в [4], где представлен и реализующий его код в MATLAB.

Основная идея данного метода заключается в оценке нагружен-ности каждого элемента и отбраковке недостаточно нагруженных элементов. Критерием нагруженности могут выступать, например, напряжения по Мизесу или энергия деформации элемента, используемая в данной работе. Результатом решения задачи ТО BESO-методом является черно-белое распределение материала. Черным областям конструкции соответствуют элементы, включенные в решение, белым — отбракованные в процессе решения.

Для того чтобы оценить значимость элемента для текущей конструкции, вводится число чувствительности элемента ае, равное энергии его деформации:

«е = "т K eu е, (15)

где ue — вектор перемещений элемента; Ke — матрица жесткости элемента.

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

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

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

— а/ +а/

где / — номер итерации.

Во время процесса оптимизации для каждой итерации устанавливается целевой объем, т. е. тот объем, которого необходимо достигнуть. Поскольку требуемый объем V может быть больше или меньше объема первоначального предположения, целевой объем на каждой итерации постепенно либо уменьшается, либо увеличивается до тех пор, пока не будет достигнуто требуемое значение. Изменение целевого объема в процессе оптимизации можно выразить так:

V1+1 = Г (1 ± ЕЯ), (17)

где ЕЯ — скорость эволюции.

Как только требуемый объем будет достигнут, он станет поддерживаться постоянным для оставшихся итераций.

Ранее было отмечено, что БЕБО-метод позволяет получать черно-белое решение вследствие наличия в рабочей области у всех элементов лишь двух дискретных состояний. При этом не происходит полного удаления из расчета отбракованных элементов, вместо этого их жесткость сильно снижается — до такого уровня, при котором они не влияют на конечно-элементный анализ. Снизить жесткость можно введением двух уровней плотности элемента и зависимости жесткости от плотности, аналогичной зависимости в 8ГМР-методе:

Еуы(*) = р(Х)рЩы, р > 1; (18)

Р(х) = Ртш или 1.

Постановка задачи ТО с целевой функцией податливости и ограничением по объему имеет вид

min l (u ); a (u, v) = l (v);

Eijkl (x) = p(x)pE0kl, p >1; (19)

J p(x)dQ < V, p(x) = Pmin или 1

Q

Определим производную целевой функции по плотности. Поскольку модель материала в BESO-методе задается подобно тому, как задавалась в SIMP-методе, производную можно записать таким же образом:

я/

Я =-ppf "1ит K eu e. (20)

Ф

Учитывая введенную модель материала, можно определить число чувствительности элементов с различной плотностью:

1 dl к еи е, если ре = 1;

=--d_ = 1 p 1 т (21)

P dPe [p^UтКeUe, если Ре = 0.

Схема обновления плотностей элементов на i-й итерации будет иметь вид

i 11, если а\ > Л;

ре=1 е (22)

[0, если ае < Л,

где Л — предельное число чувствительности.

Элементы с числом чувствительности меньшим Л отбраковываются, с большим — остаются в конструкции. Значение Л выбирается во внутреннем цикле отдельно для каждой итерации таким образом, чтобы удовлетворить ограничению по объему. Это можно сделать, например, методом бисекции.

Влияние параметров BESO-метода на решение задачи ТО. Для BESO-метода параметрами, влияющими на решения, являются радиус фильтрации rmin и скорость эволюции ER. Оценим их влияние, как и в случае применения SIMP-метода, использовав задачи о консольной и шарнирной балках с теми же размерами, параметрами КЭС и свойствами материала. Примем следующие стандартные значения соответствующих параметров: rmm = 1,5 мм, p = 3, ER = 0,02, оставив без изменений критерий завершения и максимальное количество итераций.

Результаты расчетов для консольной балки представлены в табл. 3 и на рис. 7.

Таблица 3

Результаты расчетов для консольной балки при варьировании параметров гт1„ и ЕЯ

Варьируемый параметр Значение параметра / Значение целевой функции 1, Н мм / Количество итераций N

Радиус фильтрации Гтш, мм 0,5 0,7 1,0 1,5 2,0 5,0

61,4 61,5 61,5 62,0 61,6 72,6

43 44 50 41 50 49

Скорость эволюции ЕЯ 0,005 0,01 0,02 0,05 0,07 0,1

52,7 61,7 62,0 62,1 62,4 539,8

100 79 41 33 29 19

а б

Рис. 7. Варианты конструкции при варьировании параметров гтш, мм (а)

и ЕЯ (б) (см. табл. 3)

Результаты расчетов для шарнирной балки представлены в табл. 4 и на рис. 8.

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

вается при значении гт1п = 1,5 мм, а отклонение от этого значения ухудшает сходимость. В случае консольной балки четкой зависимости сходимости от радиуса фильтрации не наблюдается.

Таблица 4

Результаты расчетов для шарнирной балки при варьировании параметров гт1„ и ЕЯ

Варьируемый параметр Значение параметра / Значение целевой функции 1, И-мм / Количество итераций N

Радиус фильтрации Гтт, мм 0,5 0,7 1,0 1,5 2,0 5,0

140,3 140,3 140,0 141,6 143,9 147,5

52 51 52 45 52 60

Скорость эволюции ЕЯ 0,005 0,01 0,02 0,05 0,1 0,25

122,6 140,0 141,6 141,9 151,1 143,7

100 77 45 40 100 100

б

Рис. 8. Варианты конструкции при варьировании параметров гт1п, мм (а)

и ЕЯ (б) (см. табл. 4)

Малое значение скорости эволюции приводит к малой скорости решения. Так, в случае и консольной, и шарнирной балки значение ЕЯ = 0,005 не позволило процессу оптимизации закончится за 100 итераций, в связи с чем необходимое количество материала не было удалено. Этим объясняется большая жесткость. При установке этого параметра на слишком высоком уровне можно наблюдать сильное искажение конструкции, уменьшение жесткости и ухудшение сходимости. Наилучшие решения обеспечивают значения ЕЯ в диапазоне от 0,01 до 0,05.

Сравнение SIMP- и BESO-методов. Наиболее очевидное различие между представленными подходами заключается в том, что SIMP-метод рассматривает плотность элемента как непрерывную величину, а BESO-метод — как дискретную. Следовательно, в первом случае решением является серая картина, которую не всегда можно интерпретировать как целостную конструкцию, а во втором — решение всегда будет черно-белым, что гораздо удобнее при проектировании.

Анализ изображений на рис. 5, а, 6, а и 7, а, 8, а показал, что BESO-метод позволяет путем варьирования радиуса фильтрации получать конструкции с разной степенью измельчения структуры, что дает возможность регулировать жесткость. В SIMP-методе радиус фильтрации оказывает большее влияние на степень серости решения, чем на степень его измельчения.

По результатам расчетов (см. табл. 1-4) можно сделать вывод, что BESO-метод в основном дает более жесткие конструкции. Это обусловлено определением более удачных структурных схем, а также тем, что решение SIMP-метода содержит элементы с промежуточной плотностью, жесткость которых занижена.

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

Следует отметить, что SIMP-метод имеет более строгое математическое обоснование. Это позволяет решать с помощью аналогичного подхода задачу ТО, например, для ортотропного материала, но необходимо будет ввести в модель материала (5) зависимость тензора Eijkl от угла армирования.

Заключение. В данной работе рассмотрены проблемы, связанные с задачей топологической оптимизации. Были представлены постановка задачи топологической оптимизации и ее особенности. Приведено описание подходов к ее решению: с помощью SIMP-метода и BESO-метода. Для этих методов рассмотрено построение решения, а также изучено влияние различных настроек на результат. Кроме того, описаны основные проблемы, возникающие при решении данной задачи, и предложены методы их устранения.

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

ЛИТЕРАТУРА

[1] Azarov A.V., Antonov F.K., Golubev M.V., Khaziev A.R., Ushanov S.A. Composite 3D printing for the small size unmanned aerial vehicle structure. Composites PartB: Engineering, 2019, vol. 169, pp. 157-163. https://doi.org/10.1016/j.compositesb.2019.03.073

[2] Blakey-Milner B., Gradl P., Snedden G. Metal additive manufacturing in aerospace: A review. Materials & Design, 2021, vol. 209. https://doi.org/10.1016/j.matdes.2021.110008

[3] Bendsoe M.P., Sigmund O. Topology Optimization: Theory, Methods and Applications. New York, Springer Verlag, 2003, 271 p.

[4] Huang X., Xie Y.M. Evolutionary topology optimization of continuum structures: methods and applications. Chichester, John Wiley & Sons, 2010, 228 p.

[5] Allaire G., Gournay F., Jouve F., Toader A.-M. Structural optimization using topological and shape sensitivity via a level set method. Control and Cybernetics, 2005, vol. 34 (1), pp. 59-80.

[6] Challis V.J. A discrete level-set topology optimization code written in MATLAB. StructMultidisc Optim, 2010, vol. 41, pp. 453-464. https://doi.org/10.1007/s00158-009-0430-0

[7] Biyikli E., To A.C. Proportional Topology Optimization: A New Non-Sensitivity Method for Solving Stress Constrained and Minimum Compliance Problems and Its Implementation in MATLAB. PLoS ONE, 2015, vol. 10 (12). https://doi.org/10.1371/journal.pone.0145041

[8] Комаров В. А. Проектирование силовых аддитивных конструкций: теоретические основы. Онтология проектирования, 2017, № 2 (24), с. 191-206. https://doi.org/10.18287/2223-9537-2017-7-2-191-206

[9] Кишов Е.А., Комаров В.А. Топологическая оптимизация силовых конструкций методом выпуклой линеаризации. Вестник Самарского университета. Аэрокосмическая техника, технологии и машиностроение, 2018, № 1, с. 137-149. https://doi.org/10.18287/2541-7533-2018-17-1-137-149

[10] Papapetrou V.S., Patel C., Tamijani A.Y. Stiffness-based optimization framework for the topology and fiber paths of continuous fiber composites. Composites PartB, 2020. https://doi.org/10.1016/j.compositesb.2019.107681

[11] Gandhi Y., Minak G. A Review on Topology Optimization Strategies for Addi-tively Manufactured Continuous Fiber-Reinforced Composite Structures. Appl. Sci., 2022, vol. 12. https://doi.org/10.3390/app122111211

[12] Федулов Б.Н., Федоренко А.Н., Антонов Ф.К., Ломакин Е.В. Алгоритм топологической оптимизации конструкции, выполненной из анизотропного материала с учетом параметров ориентации армирования. Вестник Пермского национального исследовательского политехнического университета. Механика, 2021, № 3, с. 182-189. https://doi.org/10.15593/perm.mech/202L3.17

[13] Azarov A.V., Latysheva T.A., Khaziev A.R. Optimal design of advanced 3D printed composite parts of rocket and space structures. IOP Conference Series: Materials Science and Engineering, 2020. https://doi.org/0.1088/1757-899X/934/1/012062

[14] Sigmund O. A 99 line topology optimization code written in matlab. Structural andMultidisciplinary Optimization, 2001, vol. 21, no. 2, pp. 120-127.

Статья поступила в редакцию 31.03.2023

Ссылку на эту статью просим оформлять следующим образом: Косых П. А., Азаров А.В. Теория и анализ методов топологической оптимизации. Инженерный журнал: наука и инновации, 2023, вып. 4. http://dx.doi.org/10.18698/2308-6033-2023-4-2264

Косых Павел Андреевич — студент кафедры «Ракетно-космические композитные конструкции» МГТУ им. Н.Э. Баумана. e-mail: kpa18m334@student.bmstu.ru

Азаров Андрей Валерьевич — д-р техн. наук, доцент кафедры «Ракетно-космические композитные конструкции» МГТУ им. Н.Э. Баумана, главный научный сотрудник АО Центральный научно-исследовательский институт специального машиностроения. e-mail: azarova@bmstu.ru

Theory and analysis of the topology optimization methods

© PA. Kosykh1, A.V. Azarov1,2

:Bauman Moscow State Technical University, Moscow, 105005, Russia 2Central Research Institute of Special Machinery Joint Stock Company, Khotkovo, 141371, Russia

As of today, additive technologies are making it possible to create products being close to the optimal shape. Topological optimization is widely used to design such products. The paper considers two common approaches to solving this problem: SIMP and BESO methods. Essence of the topological optimization problem, its formulation in a general form and typical examples to demonstrate this problem are described. Theoretical foundations and implementation features are presenting each method; algorithms sensitivity to the initial settings is analyzed. The problems that arise within the solution, such as the chessboard problem and dependence on the finite element mesh, are analyzed; options for solving these problems are provided. Two approaches were compared. It is concluded that the BESO method is offering more efficient and more design-friendly solutions.

Keywords: topological optimization, optimal design, variable density body, SIMP, BESO, additive technologies, finite element method

REFERENCES

[1] Azarov A.V., Antonov F.K., Golubev M.V., Khaziev A.R., Ushanov S.A. Composite 3D printing for the small size unmanned aerial vehicle structure. Composites PartB: Engineering, 2019, vol. 169, pp. 157-163. https://doi.org/10.10167j.compositesb.2019.03.073

[2] Blakey-Milner B., Gradl P., Snedden G. Metal additive manufacturing in aerospace: A review. Materials & Design, 2021, vol. 209. https://doi.org/10.1016/j.matdes.2021.110008

[3] Bendsoe M.P., Sigmund O. Topology Optimization: Theory, Methods and Applications. New York, Springer Verlag, 2003, 271 p.

[4] Huang X., Xie Y.M. Evolutionary topology optimization of continuum structures: methods and applications. Chichester, John Wiley & Sons, 2010, 228 p.

[5] Allaire G., Gournay F., Jouve F., Toader A.-M. Structural optimization using topo-logical and shape sensitivity via a level set method. Control and Cybernetics, 2005, vol. 34 (1), pp. 59-80.

[6] Challis V.J. A discrete level-set topology optimization code written in MATLAB. StructMultidisc Optim, 2010, vol. 41, pp. 453-464. https://doi.org/10.1007/s00158-009-0430-0

[7] Biyikli E., To A.C. Proportional Topology Optimization: A New Non-Sensitivity Method for Solving Stress Constrained and Minimum Compliance Problems and Its Implementation in MATLAB. PLoS ONE, 2015, vol. 10 (12). https://doi.org/10.1371/journal.pone.0145041

[8] Komarov V.A. Proektirovanie silovykh additivnykh konstruktsiy: teoreticheskie osnovy [Theoretical basis for design of load-bearing structures produced using additive technologies]. Ontologiya proektirovaniya — Ontology of designing, 2017, no. 2 (24), pp. 191-206. https://doi.org/10.18287/2223-9537-2017-7-2-191-206

[9] Kishov E.A., Komarov V.A. Topologiya optimizatsii silovykh konstruktsiy metodom vypukloy linearizatsii [Topology optimization of a load-bearing struc-

ture via the method of convex linearization], Vestnik Samarskogo universiteta. Aerokosmicheskaya tekhnika, tekhnologii i mashinostroenie — Vestnik of Samara University. Aerospace and Mechanical Engineering, 2018, no, 1, pp, 137-149, https://doi.org/10.18287/2541-7533-2018-17-1-137-149

[10] Papapetrou V.S., Patel C., Tamijani A,Y, Stiffness-based optimization framework for the topology and fiber paths of continuous fiber composites, Composites PartB, 2020, https://doi,org/10,1016/j,compositesb,2019,107681

[11] Gandhi Y,, Minak G, A Review on topology optimization strategies for additive-ly manufactured continuous fiber-reinforced composite structures, Appl. Sci., 2022, vol, 12, https://doi,org/10,3390/app122111211

[12] Fedulov B,N,, Fedorenko A,N,, Antonov F,K,, Lomakin E,V, Algoritm topolog-icheskoy optimizatsii konstruktsii, vypolnennoy iz anizitropnogo materiala s uchetom parametrov orientatsii armirovaniya [Algorithm for Topology Optimization of a Structure Made of Anisotropic Material with Consideration of the Reinforcement Orientation Parameters], Vestnik Permskogo natsionalnogo is-sledovatelskogo politekhnicheskogo universiteta. Mekhanika — PNRPU Mechanics Bulletin, 2021, no, 3, pp, 182-189,

https://doi,org/10,155 93/perm, mech/2021,3,17

[13] Azarov A,V,, Latysheva T,A,, Khaziev A,R, Optimal design of advanced 3D printed composite parts of rocket and space structures, IOP Conference Series: Materials Science and Engineering, 2020, https://doi,org/0,1088/1757-899X/934/1/012062

[14] Sigmund O, A 99 line topology optimization code written in matlab, Structural andMultidisciplinary Optimization, 2001, vol, 21, no, 2, pp, 120-127,

Kosykh P.A., Student, Department of Space and Rocket Composite Structures, Bauman Moscow State Technical University, e-mail: kpa18m334@student,bmstu,ru

Azarov A. V., Dr, Sc, (Eng,), Associate Professor, Department of Space and Rocket Composite Structures, Bauman Moscow State Technical University; Principal Research Scientist, Central Research Institute of Special Machinery, e-mail: azarov@anisoprint,ru