Full Text

Введение Разработка методов построения форм тел, оптимальных по сопротивлению и/или глубине проникания в общей постановке, возможна лишь при наличии упрощающих предположений о характере взаимодействия тела и среды [1-3], например гипотезы локальности [4-9]. В основе получивших широкое распространение [9-14] моделей локального взаимодействия (МЛВ) лежит предположение, что каждый элемент контактной поверхности взаимодействует со средой независимо от других участков тела. Применение различных модификаций МЛВ позволило выделить [9, 15-18] класс абсолютно оптимальных тел (АОТ), содержащий общее решение достаточно широкого круга задач оптимизации формы тела, проникающего в грунтовые среды. В пространственной постановке нормаль к поверхности таких тел составляет с направлением движения постоянный оптимальный угол [15]. Среди тел вращения АОТ является круговой конус с тем же оптимальным углом раствора, определяемым скоростью тела и прочностными характеристиками среды. Для грунтов малой прочности осесимметричные АОТ обладают значительным удлинением. Для поиска формы тела вращения минимального сопротивления внедрению при заданной длине тела, меньшей длины абсолютно оптимального конуса, которая, как известно, содержит торец [16], необходимо применять численные методы [19]. Известные методы решения задач поиска образующей тел с притуплением, внедряющихся на максимальную глубину, основаны на использовании МЛВ. Постановка и численный метод решения вариационных задач на основе генетического алгоритма изложен в работе [19], разработке модификаций метода локальных вариаций посвящены работы [20, 21], где рассматривались задачи проникания оптимальных тел с торцом в среды типа металла или бетона. Для грунтовых сред экспериментально и теоретически - сравнением с результатами численных расчетов в осесимметричной постановке на основе модели грунтовой среды Григоряна - была показана применимость МЛВ к описанию проникания острых конусов [22-23]. Погрешность МЛВ на базе решения задачи о расширении сферической полости в определении сил сопротивления внедрению затупленных тел исследована в работах [14, 23]. Показано (на примере задачи поиска оптимального тела вращения заданной длины и площади основания) [6], что решение с учетом нелинейных эффектов обтекания в двумерной постановке позволяет существенно уточнить как форму, так и силовые, и кинематические характеристики оптимальных затупленных тел при проникании в грунтовые среды. В данной статье аналогично [24, 25] сила сопротивления внедрению в грунт - целевая функция в задаче оптимизации формы тела, определяется на базе упрощенной модели контактного взаимодействия при различных параметризациях образующей [26-28]. Далее полученное решение принимается в качестве начального приближения в итерационном процессе построения оптимальной формы ударника при движении в грунтовых средах на основе численных расчетов в осесимметричной постановке. 1. Постановка задачи и методика численного решения Рассматривается стадия развитого проникания ударника при постоянной скорости внедрения V0, причем обтекание искомой головной части ударника грунтом плотностью ρ0 происходит безотрывно. Введем цилиндрическую систему координат r0z, 0z - ось симметрии. Ударник - тело вращения заданной длины L с радиусом основания R, образующая боковой поверхности тела описывается кривой r(z) и включает передний торец радиусом R0. Действующее контактное давление определяется только нормальной компонентой вектора скорости в соответствии с МЛВ общего вида где η - угол, образованный направлением вектора скорости и касательной к искомой кривой r(z). Полное осевое сопротивление F проникающего тела вращения [19, 27, 28] с плоским передним торцом имеет вид . Положим , - площадь миделя; - радиус притупления; - удлинение тела, тогда . (1) Таким образом, при известной функции имеем простейшую задачу вариационного исчисления, в которой необходимо определить такую выпуклую кривую , удовлетворяющую условию , что полное сопротивление (1) имеет минимум. Приближенное решение этой задачи на основе известного метода локальных вариаций [29] сводится к отысканию чисел, минимизирующих функцию N переменных , где - искомые точки кусочно-линейной аппроксимации кривой при известном соответствующем распределении координат (см. также [20, 21]). Ранее была показана [28] сходимость последовательных приближений при кусочно-линейном представлении образующей боковой поверхности оптимального тела и убывание погрешности определения силы сопротивления с ростом числа узлов образующей. В задачах газовой и аэродинамики доказана эффективность аппроксимации образующей в классе полиномов или степенных функций, содержащих небольшое число искомых параметров оптимизации [26]. Отмечено [27], что аппроксимация квадратичным полиномом в форме Безье является удачной параметризацией образующей оптимального тела и в случае проникания его в грунт. Модифицированный алгоритм [27, 28] позволяет единообразно проводить вычисления как при кусочно-линейном представлении образующей, так и при функциональном, когда компоненты вектора являются параметрами полинома. Для применения алгоритма необходимо конкретизировать способ расчета интегральной силы сопротивления внедрению D в соотношении (1) - на основе МЛВ или одной из известных моделей грунтовой среды в постановке механики сплошных сред [14, 23]. Вектор неизвестных параметров r = , является одним из локальных экстремумов (1), для определения которых составляется система уравнений D'(r) = 0, . (2) Система нелинейных уравнений (2) решается итерационным методом Ньютона: D''(rk) (rk+1 - rk) = -D'(rk), , , (3) где k - номер итерации; D''(rk) - симметричная матрица коэффициентов (в случае кусочно-линейного представления образующей имеет трехдиагональный вид). Частные производные первого и второго порядка могут вычисляться точно или методом конечных разностей с шагом варьирования [27, 28]. Проведение вычислений на каждой итерации (3) эквивалентно нахождению минимума приращения сопротивления , представленного в виде квадратичной формы разложением в ряд Тейлора в окрестности rk: . (4) Другие полиномиальные интерполяции целевой функции в области изменения параметров радиусом могут быть получены методами планирования эксперимента [30]. Наиболее полная модель, основанная на полном ортогональном плане второго порядка, включает все промежуточные взаимодействия факторов, но требует для реализации проведения вычислений целевой функции и при больших N является, очевидно, избыточной. В ротатабельных планах второго порядка число расчетов не меньше чем . Для вычисления коэффициентов (4) на основе центральных разностей второго порядка аппроксимации также может быть составлена матрица планирования, количество операций равно . Существенно уменьшить количество параметров N, как было отмечено выше, позволяет представление образующей оптимального тела кривой Безье [31]. Параметрический полином в форме Безье описывается выражением . (5) Рис. 1. Представление образующей параметрическим полиномом Безье С учетом поставленных ограничений - заданных безразмерных величин радиуса основания и длины искомого оптимального тела , z0 = 0, z2 = l, вектор параметров r содержит три параметра (N = 3), представляющих неизвестные координаты вершин , характеристического треугольника полинома Безье (рис. 1), и - углы при соответствующих вершинах характеристического треугольника. Таким образом, на каждой итерации метода Ньютона (3) целевую функцию необходимо вычислять 15 раз. Это не критично при аналитическом вычислении функции, но существенно увеличивает трудоемкость при численном расчете обтекания ударника даже в осесимметричной постановке. В данной работе аналогично [24] предварительное исследование точности и сходимости методов оптимизации формы ударника, обоснование выбора приемлемой аппроксимации образующей проводятся в упрощенной постановке на основе МЛВ. Как будет показано далее, решение задачи параметрической оптимизации методом циклического покоординатного спуска с изменением шага варьирования позволит уменьшить число вычислений. 2. Результаты предварительных расчетов В качестве тестовой рассмотрим задачу проникания в грунт, поведение которого описывается моделью нелинейно-сжимаемой жидкости. Как было показано ранее [14], в этом случае реализуется наибольшая ошибка МЛВ при определении силы сопротивления внедрению. Сжимаемость грунтовой среды характеризуется ударной адиабатой: , (6) представленной линейной зависимостью скорости ударной волны с от массовой скорости u за ее фронтом. Постоянная А характеризует скорость распространения волны сжатия в грунте при малых давлениях, l - предельную сжимаемость грунта. Численные расчеты проводятся при следующих значениях параметров: А = 460 м/с; l = 2,3; r0 = 1700 кг/м3, размеры тела R = 0,01 м; L/R = скорость проникания = 400 м/с. Поведение нелинейно-сжимаемой жидкости описывается одночленной МЛВ . Среднее значение константы С, определенной из решения одномерной задачи о расширении сферической полости [3, 14, 32-34], , , , при примем равным С 1,1. На рис. 2 приведены результаты расчетов оптимальных форм: кривые 1 есть численные решения системы уравнений (3) на основе МЛВ, кривые 2 получены в осесиметричных численных расчетах (см. далее). Светлые квадраты соответствуют кусочно-линейному представлению образующей (N = 16), сплошная и штриховая линии 1 - представлению образующей в форме полинома Безье при N = 3 и 2. Видно, что решение, определяемое тремя параметрами, практически точно соответствует кусочно-линейному представлению (фрагмент в окрестности вершины тела см. на рис. 2). Пунктирная линия также является решением системы (3) при N = 3 и получена в сходящемся итерационном процессе метода Ньютона, но при другом начальном приближении. Из двух решений локальному минимуму соответствует только образующая, изображенная сплошной линией (все угловые миноры матрицы вторых производных D'' в соотношениях (3) положительны). Другое численное решение, локальным минимумом не является - угловые миноры матрицы D'' имеют разные знаки, и значение силы сопротивления для него несколько больше. Рис. 2. Образующие оптимальных тел вращения Параметры решений приведены в табл. 1, жирным шрифтом выделены искомые параметры полинома Безье. Таблица 1 Параметры оптимального тела в рамках МЛВ N P0(r,z) P1(r,z) P2(r,z) φ0 φ2 D 3 (0,096; 0) (0,419; 0,356) (1; 2,236) 0,834 1,271 0,1485 3 (0,061; 0) (0,280; 0) (1; 2,236) 0 1,259 0,1494 2 (0; 0) (0,298; 0) (1; 2,236) 0 1,267 0,1507 Обозначим оптимальные значения в первой строчке табл. 1 верхним индексом * (= 0,096, …, = 0,1485). На рис. 3 представлены линии равного уровня целевой функции полученные при фиксированных значениях радиуса притупления = 0,046; 0,052 и =0,096. Углы и принимают значения в интервалах [0; 1] и [1,25; 1,3] соответственно. Изолинии даны с шагом = 0,03 до значения = 1,02 (штриховая линия на рис. 3, а, б) и с шагом =0,2 при >1,02; на рис. 3, в используется шаг 0,1 до значения =1,3 и вдвое больший при >1,3, штриховая линия также соответствует линии раздела. Видно, что при 0,052 функция по переменным и имеет два минимума: краевой и локальный (см. рис. 3, б). При = 0,046 имеется только краевой минимум при = 0 (см. рис. 3, а), при = достигается только локальный минимум (см. рис. 3, в), являющийся также и глобальным. Отметим также, что минимум по переменной может быть локализован лишь при 0,01. а б в Рис. 3. Линии равного уровня значений целевой функции, полученные при фиксированных значениях радиуса притупления При неудовлетворительном начальном приближении итерационный процесс метода Ньютона (3) может сходиться к другому решению (строка 2 табл. 1) или расходиться. В лучшем случае решение системы (2) может быть получено за 5-7 итераций при числе расчетов целевой функции порядка 100. Рассмотрим в связи с этим применение метода циклического покоординатного спуска с переменным шагом варьирования параметров δ. Искомыми параметрами будем считать радиус притупления и углы и характеристического треугольника полинома Безье. Начальное значение шага = 0,05, начальное значение =0; = 1; = 1,35, что не является хорошим начальным приближением метода Ньютона. Последовательные приближения осуществляем, начиная с параметра , значение которого уменьшается от начального 1,35 до 1,25 (рис. 4, а, кривая 1), точки соответствуют вычисленным значениям целевой функции , отложенным по оси ординат. Варьирование параметра (рис. 4, б) приводит к установлению = 0 (см. также рис. 3, а), радиус притупления = 0,05 (рис. 4, в). Результаты второго цикла вычислений с = 0,05 шагом и полученными параметрами показаны на рис. 4, а-в кривыми 2. Отмечается уменьшение относительного значения сопротивления 1 %, но значения параметров далеки от оптимальных значений (приблизительно соответствуют приведенным в строке 2 табл. 1). Рис. 4. Изменение целевой функции в зависимости от параметров образующей на различных этапах покоординатного спуска (расчет в рамках МЛВ) Дальнейшие изменения в решении происходят при уменьшении шага варьирования 0,01. В этом случае в одном из расчетов значение становится большим 0,52 (см. рис. 3, б), и параметр достаточно быстро принимает значение, близкое (см. кривая 3 на рис. 4, б), как и остальные параметры в 3-5 циклах (рис. 4, кривые 3-5, светлые точки). В целом можно отметить, что общее число расчетов в методе спуска с точностью = 0,01 даже при неудовлетворительном начальном приближении оказалось не превышающим таковое в методе Ньютона. Таким образом, для расчетов оптимальных форм на базе уравнений механики сплошных сред в осесимметричной постановке можно рекомендовать аппроксимацию образующей полиномом Безье с тремя параметрами и метод оптимизации на основе циклического покоординатного спуска. 3. Результаты расчетов в осесимметричной постановке Система уравнений модели нелинейной сжимаемой жидкости в эйлеровых переменных в цилиндрической системе координат записывается следующим образом: , , , , , где и - начальная и текущая плотность; d/dt - полная производная по времени; символ после запятой означает дифференцирование по соответствующей переменной, константы А и определяются ударной адиабатой (6). Для определения сопротивления тела решается задача обтекания методом установления. В начальный момент времени давление в грунте равно нулю, скорость частиц грунта равна скорости удара и направлена вдоль оси вращения тела. Внешние границы расчетной области грунта отнесены достаточно далеко, чтобы возмущения, отражающиеся от них, не исказили численное решение в области контакта тела и среды. Неподвижный ударник считается жестким. На головной части ударника, контактирующей с грунтовой средой, принимается условие «непроницаемости» по нормали со «скольжением по касательной». Определенная таким образом сила сопротивления соответствует силе на квазистационарной стадии внедрения с постоянной скоростью, полученной в прямом ударе. Численные расчеты проводятся в рамках методики, основанной на модифицированной схеме Годунова первого порядка точности, реализованной в пакете прикладных программ «Динамика-2» [36]. Применимость известной численной методики и пакета программ к расчету параметров процессов проникания жестких и деформируемых ударников в мягкие грунтовые среды показана ранее [37]. Численные расчеты проводятся при тех же значениях параметров среды, что и в п. 2. Анализ точности используемой численной методики показал практически линейную монотонную сходимость решений (рис. 5) в зависимости от шага сетки , где H - длина образующей тела; n - число узлов разностной сетки. На рис. 5, а треугольниками, ромбами и точками обозначены соответственно рассчитанные численно значения сопротивления конуса и оптимальных тел, полученных в рамках МЛВ () и в осесиметричных численных расчетах (, см. далее). а б Рис. 5. Анализ сходимости результатов численных расчетов проникания тел вращения при изменении шага сетки: а - безразмерные силы сопротивления внедрению; б - относительные погрешности определения сопротивления Значение сопротивления каждого тела i =1, 2, 3 отнесено к сопротивлению конуса , определенного в рамках МЛВ: . Пунктирные линии соответствуют линейным аппроксимациям методом наименьших квадратов. Обозначим прогнозируемые при h = 0 значения сопротивления для каждого тела , i =1, 2, 3. На рис. 5, б (обозначения, как на рис. 5, а) приведены относительные погрешности определения сопротивлений , i = 1, 2, 3. При одинаковом размере ячейки разностной сетки h (выбирались близкие к квадратным ячейки сетки) отмечается большая погрешность для затупленного тела. Далее приводятся результаты определения параметров образующей оптимального тела при проникании в грунт на основе численных осесимметричных расчетов. Начальное приближение применяемого метода циклического покоординатного спуска приведено в строке 1 табл. 1, =0,05, порядок варьирования параметров ---, как в п. 2. На рис. 6 приведены значения целевой функции в зависимости от изменяемых в процессе покоординатного спуска параметров. Видно, что процесс последовательных приближений качественно близок к приведенному ранее на рис. 4. В первом цикле вычислений значение параметра уменьшается до 0 (рис. 6, б), видно, что сопротивление незначительно меняется при 0,2 < < 0,3 (рис. 6, в). Рис. 6. Изменение целевой функции в зависимости от параметров образующей на различных этапах покоординатного спуска (расчеты в осесимметричной постановке) Далее полагаем = 0,2. Результаты второго цикла вычислений с полученными параметрами показаны на рис. 6, а, б кривыми 2, отмечается уменьшение относительного значения сопротивления . Кривые 3 на рис. 6 соответствуют расчетам при = 0,3. Величина осталась практически без изменений, значение сопротивления тела равно минимальному при φ0 = 0,434. Результаты расчетов сведены в табл. 2. Таблица 2 Параметры оптимального тела в осесимметричной постановке N P0(r,z) P1(r,z) P2(r,z) φ0 φ2 D 3 (0,096; 0) (0,419; 0,356) (1; 2,236) 0,834 1,271 0,1299 3 (0,296; 0) (0,688; 0,179) (1; 2,236) 0,434 1,42 0,0797 3 (0,196; 0) (0,66; 0) (1; 2,236) 0 1,42 0,0797 В первой строке табл. 2 приведены параметры оптимального тела, полученного в п. 2 в рамках МЛВ, отличие составляет значение сопротивления D, которое несколько уменьшилось. Вторая и третья строки представляют искомые параметры образующей оптимальных тел, обладающих одинаковым сопротивлением (в рамках погрешности численных расчетов порядка 1 %). Образующие тел с этими параметрами приведены на рис. 2 сплошной и пунктирной кривыми 2 соответственно. Отметим качественно похожее поведение численных решений в постановках на основе МЛВ и механики сплошных сред, но существенное отличие в формах оптимальных тел (см. рис. 2) и их сопротивлении (см. рис. 5, а). Как отмечалось выше, полученное различие в силах сопротивления оптимальных тел для модели грунта - нелинейной жидкости является наибольшим. Различие уменьшается при учете внутреннего и поверхностного трения, а также других особенностей динамического деформирования грунтовых сред, численная реализация которых не представляет труда в рамках известных моделей грунтовых сред [1]. Заключение Представлены результаты численного расчета методом локальных вариаций оптимальных форм осесимметричных ударников при движении в грунтовой среде. Установлено, что методика расчета полей давлений в окрестности головной части тела должна учитывать двумерные эффекты обтекания, однако предварительные исследования целесообразно проводить на базе известных упрощенных моделей, например локального взаимодействия. Показано, что метод циклического покоординатного спуска позволяет с точностью, определяемой погрешностью численных расчетов, вычислять параметры образующей тела вращения, которая может быть представлена полиномом Безье второй степени. Разделы 1, 2 статьи выполнены при частичном финансировании в рамках базовой части государственного задания Минобрнауки РФ (проект № 2014/134 2226), раздел 3 при поддержке РНФ (грант №15-19-10039).

About the authors

V G Bazhenov

Research Institute of Mechanics Lobachevsky Nizhniy Novgorod State University

V L Kotov

Research Institute of Mechanics Lobachevsky Nizhniy Novgorod State University

E Yu Linnik

Research Institute of Mechanics Lobachevsky Nizhniy Novgorod State University

References

Statistics

Views

Abstract - 138

PDF (Russian) - 63