Динамические задачи строительной механики с отрицательным временем

Полный текст

1. Введение Динамические задачи строительной механики предполагают моделирование поведения системы при нарастающем времени. В данном исследовании мы рассматриваем обратную в определенном смысле задачу, когда уравнения движения используются для восстановления предыдущих по времени состояний системы. Отметим, что подобные задачи известны в механике, например, в небесной механике развиты методы предсказаний не только будущих солнечных затмений, но и определение времени ранее наблюдавшихся затмений, что позволяет уточнять хронологию исторических событий. Как пример задачи с отрицательным временем для упругой системы приведем аналитическое решение задачи о затухающих колебаниях системы с одной степенью свободы (рис. 1). Известное частное решение которой имеет вид y = Be-ntsinω .t (1). а б Рис. 1. Затухающие колебания системы с одной степенью свободы: а - график аналитического решения; б - система с одной степенью свободы И с т о ч н и к : выполнено В.Б. Зылевым, А.В. Штейном Figure 1. Damped vibration of a single degree-of-freedom system: a - analytical solution; б - single degree-of-freedom system S o u r c e : compiled by V.B. Zylev, A.V. Shtein В данном случае для восстановления предыдущих состояний системы (светлая часть графика на рис. 1, а) достаточно подставить в решение (1) отрицательное значение времени. В настоящее время широко используются численные методы интегрирования уравнений движения. Наибольшее количество публикаций относится к неявным вычислительным схемам (метод Ньюмарка, θ - метод Вилсона и др.) [1-4]. Кроме неявных вычислительных схем используются и явные вычислительные схемы (метод центральных разностей, метод Рунге - Кутта, метод Адамса [5]). Применение любого численного метода неизбежно связано с погрешностями, определение и оценка которых относятся к основным вопросам применимости того или иного вычислительного алгоритма. Методы контроля весьма разнообразны: счет с уменьшенным шагом; сравнение с известными аналитическими решениями, сравнение с экспериментальными данными и т.д. Мы хотим обратить внимание на ту дополнительную возможность, которую здесь предоставляет решение обратной динамической задачи, или, лучше сказать, решение задачи с отрицательным временем. Следует заметить, что термин «обратные задачи» в задачах акустики имеет совершенно другое значение [6-10]. Обратные задачи акустики заключаются в определении упругих и массовых характеристик массива грунта по параметрам распространяющихся в нем волнам деформаций. Штейн А.В., Зылев В.Б. Строительная механика инженерных конструкций и сооружений. 2024. Т.20. № 3. . 276-288 При численном решении динамической задачи погрешность при прямом и обратном течении времени получается различной. Это открывает возможность контроля точности решения. По сравнению с обычным временем отрицательное время в задачах контроля точности имеет то преимущество, что заранее известно состояние системы, к которому она должна прийти (точнее говоря, вернуться). Например, это может быть начальное ненапряженное состояние системы. Практическое приложение решений задач с отрицательным временем мы видим в контроле точности численных методов. Рассмотрим далее особенности выполнения такого контроля. 2. Рассматриваемый численный метод интегрирования Возможность построения решения при обратном течении времени в принципе не связана с рассмотрением какого-либо конкретного численного метода интегрирования уравнений движения, соответственно, любой вычислительный алгоритм может рассматриваться как тестируемый подобным способом. Приводимые ниже примеры решаются нами с использованием явной экстраполяционной схемы шагового интегрирования, которая достаточно подробно изложена в [11-13]. Все рассмотренные в данной статье примеры относятся к плоской задаче теории упругости с возможностью рассмотрения больших перемещений. Использованы треугольные конечные элементы с линейной аппроксимацией перемещений внутри элемента. Сетка разбиения берется с постоянным шагом. Массы предполагаются сосредоточенными в узлах. Условно устойчивая вычислительная схема использует весьма малый шаг по времени, что дает возможность получать волновые решения. Приведем крайне сокращенное описание численного метода, что необходимо для понимания особенностей решения при прямом и обратном течении времени. На рис. 2 показаны графики изменения ускорения (составляющей ускорения по одной из координатных осей) некоторого узла системы. Рассмотрим сначала обычное течение времени. На рис. 2 - это шаг интегрирования по времени. W1 , W2 , W3 - это значения ускорений на уже пройденных участках интегрирования. W0 - ускорение, которое только что вычислено по текущим координатам узлов и значениям скоростей. а б Рис. 2. График изменения составляющей ускорения: а - прямое течение времени; б - обратное течение времени И с т о ч н и к : выполнено В.Б. Зылевым, А.В. Штейном Figure 2. Change of acceleration component with time: a - forward flow of time; б - backward flow of time S o u r c e : compiled by V.B. Zylev, A.V. Shtein ДИНАМИКА КОНСТРУКЦИЙ И СООРУЖЕНИЙ . 2024;20(3):276-288 Через четыре известных значения ускорений проводится кубическая парабола, что обеспечивает получение приближенного выражения для ускорения на предстоящем шаге интегрирования. Используемые при расчете формулы для приращения скорости и смещения за шаг по времени имеют вид ΔV 55W0 59W1 37W2 9W3 ; 24 ΔS V0Δ Δ2 323W0 264W1 159W2 38W3 . (2) 360 Один шаг по времени при обратном течении времени делается аналогично прямому. Значения ускорений в опорном отрезке будут W0/ , W1/ , W W2/ , 3/ . Конечно, эти значения получены при обратном течении времени (см. рис. 2, б). Таким образом, шаг вперед и шаг назад по времени будут иметь каждый свою погрешность, что определяет возможность использовать отрицательное время для контроля точности. 3. Пример численного решения задачи о пластине (балке) при внезапном приложении нагрузки На рис. 3 показана расчетная схема консольной балки (пластины). Размеры балки - высота 2 м, длина 20 м. Балка загружается 11 сосредоточенными силами 200 кН, направленными вниз. Точки приложения сил по верхней границе балки отмечены на рис. 3, а. Глухая заделка закрепляет пластину в правом сечении балки. Конечные элементы - равносторонние прямоугольные треугольники в 20 слоев по вертикали и 200 слоев по горизонтали. Таким образом, имеем 8000 конечных элементов. Материал имеет модуль упругости 200 000 МПа, коэффициент Пуассона 0,25. Удельная плотность 7850 кг/м3. Толщина расчетного слоя 0,01 м. Представлена плоская деформация. Рассматриваются большие перемещения, но соотношения между деформациями и напряжениями подчиняются закону Гука, таким образом, задача условно решается с учетом лишь геометрической нелинейности, независимо от уровня напряжений. Рис. 3. Внезапное приложение нагрузки к консольной балке: а - расчетная схема с разбивкой на конечные элементы, статический прогиб; б - полученные в расчете максимальные отклонения вниз; в - полученные в расчете максимальные отклонения вверх И с т о ч н и к : выполнено В.Б. Зылевым, А.В. Штейном Figure 3. Instantaneous load applied to a cantilever beam: a - finite element model; б - maximum upwards deflection; в - maximum downwards deflection S o u r c e : compiled by V.B. Zylev, A.V. Shtein Штейн А.В., Зылев В.Б. Строительная механика инженерных конструкций и сооружений. 2024. Т.20. № 3. . 276-288 В расчетной модели учитываются силы внешнего сопротивления. В каждый узел прикладывается противоположно его скорости сила сопротивления -kmy , где m - масса узла. Значение k принято равным 10 с-1. На рис. 4 показан график изменения во времени перемещения отмеченной на рис. 3, а точки торцевого поперечного сечения стержня по вертикали. Красными стрелками (направленными направо) на рис. 4 отмечено обычное течение времени, голубой стрелкой (направленной налево) отмечено обратное течение времени. Рис. 4. График изменения перемещения торцевой точки стержня по вертикали И с т о ч н и к : выполнено В.Б. Зылевым, А.В. Штейном Figure 4. Vertical displacement of the end of the bar with time S o u r c e : compiled by V.B. Zylev, A.V. Shtein Итак, как следует из графика на рис. 4, в течение одной секунды с момента приложения нагрузки почти достигается положение равновесия балки, которое показано на рис. 3, а, где перемещения показаны в натуральном масштабе. После момента времени 1 с обычный расчет заканчивается и начинается расчет при обратном течении времени (второй участок графика на рис. 4). В нашей компьютерной программе график при обратном течении времени строится с нарастанием времени направо. При таком способе построения графиков не происходит наложения графиков, что обеспечивает их удобное сравнение и взаимный контроль. В идеале график при обратном течении времени будет зеркальным отображением графика при прямом течении времени, что и наблюдается на рис. 4. На втором участке (рис. 4) можно наблюдать возвращение системы в исходное недеформированное состояние (график коснулся оси времени). Если в этот момент снять нагрузку, исходное равновесие восстановится. В данном примере продемонстрирован переход за начальный момент времени при сохранении нагрузки. Поскольку при прямом течении времени амплитуды убывали (затухающие колебания), при обратном течении времени они будут неограниченно возрастать, что можно видеть на рис. 4. В данном случае на возрастание амплитуды отведено 0,6 с, при этом амплитуды стали нереально большими, система перешла в состояние весьма больших прогибов (см. рис. 3, б), которые тем не менее отображаются расчетной моделью с точки зрения геометрических соотношений. Далее (см. рис. 4) опять включается расчет с обычным течением времени, при этом система вновь проходит исходное равновесное состояние и далее повторяет движение, соответствующее обычному динамическому расчету. Насколько можно видеть из графиков на рис. 4, описанная картина движения правильно отображается расчетом. ДИНАМИКА КОНСТРУКЦИЙ И СООРУЖЕНИЙ . 2024;20(3):276-288 Для более точной оценки погрешности численного решения рассмотрим данные, приведенные в табл. 1, где даны локальные минимумы перемещений для точек 1, 2, 3 (см. рис. 4). Все эти точки соответствуют тождественным состояниям системы. Таблица 1 Локальные минимумы на графике изменения перемещения торцевой точки балки для идентичных состояний 1, 2, 3 (рис. 4) Номер локального минимума на рис. 4 1 2 3 Значения перемещения по вертикали, метры -3,833306 -3,833305 -3,833432 И с т о ч н и к : выполнено выполнено В.Б. Зылевым, А.В. Штейном Table 1 Local minima on the bar end displacement curve for identical states 1, 2, 3 (Figure 4) Local minima ID in Figure 4 1 2 3 Vertical displacement, m -3.833306 -3.833305 -3.833432 S o u r c e : compiled by V.B. Zylev, A.V. Shtein а б в Рис. 5. Отрезок времени 0,25 с: а - вертикальное перемещение узла (фрагмент рис. 4); б - составляющая скорости узла по вертикали; в - составляющая ускорения по вертикали И с т о ч н и к : выполнено В.Б. Зылевым, А.В. Штейном Figure 5. Time interval of 0.25 s: a - vertical nodal displacement (Figure 4 fragment); б - vertical component of nodal velocity; в - vertical component of acceleration S o u r c e : compiled by V.B. Zylev, A.V. Shtein Как видно из табл. 1, отличия в прогибах имеются лишь в пятом знаке, что явилось неожиданностью даже для авторов статьи. Решение велось с шагом по времени = 0,00000256 с, таким образом, до третьего максимума компьютер выполнил около 1 миллиона шагов. Отметим, что шаг в нашей компьютерной программе назначается автоматически на основании фор- мулы =Тпарц , где Тпарц - наименьший парци20 альный период системы [12]. Наша явная вычислительная схема является условно устойчивой. Попытка увеличения этого шага в 2,5 раза при- водит к неустойчивому вычислительному про- цессу. Таким образом, данные, приведенные в табл. 1, соответствуют стандартному значению шага по времени. Сравнение напряжений при прямом и обратном счете приводит к большим погрешностям, которые содержатся уже в четвертой значащей цифре. Можно констатировать, что выполненный тест дал вполне удовлетворительный результат для напряжений и перемещений. Теперь рассмотрим скорость по вертикали того же узла и соответствующую составляющую ускорения. На рис. 5 показан временной интервал 0,25 с от момента внезапного приложения наг- рузки и обратный счет тоже 0,25 с. Графики на рис. 5, а это перемещение, скорость и ускорение торцевого узла по вертикали. Таким образом, здесь имеется график функции, ее первой и вто- Штейн А.В., Зылев В.Б. Строительная механика инженерных конструкций и сооружений. 2024. Т.20. № 3. . 276-288 рой производной, что на первый взгляд может вызвать сомнения. Однако рассмотрение графика перемещения в сильно увеличенном масштабе эти сомнения аннулирует, так как показывает, что эта кривая отнюдь не является гладкой. Мелкие неровности на кривой перемещений отображают высокочастотные колебания узла, без моделирования которых решить задачу не получится. Подробнее этот вопрос рассмотрен в [14]. В табл. 2 даны полученные на компьютере значения экстремумов на графиках скорости и ускорения. Таблица 2 Минимальное значение скорости и максимальное значение ускорения на отрезке времени 0,25 с Шаг численного интегрирования t c, Минимальное значение скорости на исследуемом отрезке времени м/с Максимальное значение ускорения на исследуемом отрезке времени, м/с2 Прямой ход Обратный ход Прямой ход Обратный ход =t 0,00000256 -53,26545 -53,21240 43244,92 29328,13 t = 0,00000128 2 -53,28352 -53,26672 43609,32 42302,79 t = 0,00000016 16 -53,29225 -53,29226 43625,28 43625,32 И с т о ч н и к : выполнено В.Б. Зылевым, А.В. Штейном Table 2 Minimum value of velocity and maximum value of acceleration over the time interval of 0.25 s. Numerical integration step ∆

Об авторах

Алексей Владимирович Штейн

Российский университет транспорта

Email: avsh7714@yandex.ru
ORCID iD: 0009-0003-2232-5121
SPIN-код: 3150-4438

доцент кафедры строительной механики

Москва, Россия

Владимир Борисович Зылев

Российский университет транспорта

Автор, ответственный за переписку.
Email: zylevvb@ya.ru
ORCID iD: 0000-0001-5160-0389

доктор технических наук, заведующий кафедрой строительной механики

Москва, Россия

Список литературы