petuhoff Feb 24 at 00:10

Как странные формулы ТАУ заменяют 3D расчеты на СуперЭВМ, и помогают Siemens побеждать

Easy

7 min

5.6K

Matlab*Mathematics*Programming*System Analysis and Design*CAD/CAM*

Tutorial

Этот текст – дополнение ко второй части лекции про особые линейные системы.

Для переходного процесса по распространению тепла в полубесконечной стенке в лекции приводятся вот такие фрмулы:

Переходная функция:

$h(t)=Z^{-1}[H(s)]=Z^{-1}\left[\frac{W(s)}{s}\right]=Z^{-1}\left[\frac{e^{-\sqrt{\tau\cdot s}}}{s}\right]=erfc\left[\frac{1}{2}\sqrt{\frac{\tau}{t}}\right]$ $h(t) =erfc\left[\frac{1}{2}\cdot\sqrt{\frac{\tau}{t}}\right]$

Где появляется функция ошибок Гаусса:

$erfc=\frac{2}{\sqrt{\pi}}\int_y^\infty e^{-u^2}du = 1-erfy \\erfy=\frac{2}{\sqrt{\pi}}\int_0^ye^{-u^2}dy$

Весовая в относительном времени $t/ \tau$ :

$w(t)=\frac{1}{2\cdot\sqrt{\pi}}\cdot\frac{1}{\tau}\cdot\frac{\tau}{t}\cdot\sqrt{\frac{\tau}{t}}\cdot e^{-0.25\cdot\frac{ \tau}{t}}$

Что меня не престаёт удивлять и восхищать в теории автоматического управления, так это появления волшебных чисел типа $\pi$ и и одинаковых функций Гаусса в описании абсолютно разных физических процессов. Так, например, в первом случае в лекции у нас плоская и бесконечная стенка, так откуда в решения появилось число $\pi$ ? Тут нет никакой окружности в принципе! Тут просто плоская стенка и температура! Откуда в решении число $\pi$ ? Причем здесь квадратура круга? Хочется, как Станиславский, закричать: «Не верю!» В таких случаях моя рука тянется к пистолету или АК47! Автоматом в этом случае является система структурного моделирования, которую можно использовать чтобы быстро и надежно проверить гармонию алгеброй.

Возьмём простую стенку из стали и применим теоретические формулы к практическому реальному железу. Рассмотрим три точки в стенке на расстоянии 3 мм, 6 мм и 9 мм от поверхности. Заведем параметры стали в глобальные константы проекта, а также рассчитаем относительное время. В лекции видим, что графики в лекции построены в масштабе времени t/τ.

У меня получился такой глобальный скрипт проекта:

Рисунок 1. Расчет параметров для стальной стенки

Соберем сначала формулу для весовой функции, она проще и в ней нет странных функций Гаусса. И запускаем на расчет в течение 1000 секунд. Поскольку в SimInTech все блоки векторные, одна схема может рассчитать все три и более точек по толщине металла.

У меня получился вот такой результат:

Рисунок 2. Расчет весовой фунции для трех точек

Мы видим, что максимумы функции сдвинуты относительно друг друга, поскольку время на графиках реальное, а в лекции на графике было относительное время.

В лекции у нас использовался график относительного времени t/τ. В скрипте у нас есть расчет этого времени для каждой из функций, поэтому, используя фазовый портрет, мы можем построить весовые функции в относительном времени.

В этом случае мы видим, что все максимумы весовой функции для каждой из точек находятся в точке 1/6, как об этом было предсказано в лекции!

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

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

$h(t) =erfc\left[\frac{1}{2}\cdot\sqrt{\frac{\tau}{t}}\right]$

erfc – эта странная функция Гауса в SimInTech тоже есть. А на схеме у нас уже есть место, где рассчитывается $\sqrt{\tau / t}$ Поставим на схему блок «язык программирования», в котором выполним деление на 2 и вычисление функции Гаусса для все точек, которые мы задали в скрипте.

Рисунок 3. Расчет переходной функции через функцию ошибок Гаусса

График относительной температуры для ступенчатого воздействия:

Рисунок 4. График относительной температуры по толщине пластины

Получается что одна формула позволяет посчитать температурное поле по всей пластине. Это же задача для 3D-расчета и СуперЭВМ!

Но самое удивительное нас ждет дальше.

Мы не просто так рассчитывали в скрипте относительное время. Для каждой точки по толщине оно разное. Давайте выведем это время на схему и построим график зависимости относительной температуры от относительного времени с использованием фазного портрета. У меня получилась примерно такая схема.

Рисунок 5. Построение переходного процесса в относительном времени

Я не верею своим глазам! Тут не только совпали максимумы функций, но все три графика слились в один!

Это реально одна формула, которая позволяет рассчитать температуру в любой точке по толщине пластины. В любой! Это просто магия!

Сохраним эту модель под именем САР.prt – система автоматического регулирования.

Чтобы проверить правильность, выполним сравнение с моделью, где у нас идет послойный расчет температуры в каждом слое стальной пластины.

Создаем новую модель с использование библиотеки теплогидравлики, где у нас есть блоки: тепловое граничное условие, плоская толстая пластина.

Ставим эти два блока на схему.

Добавляем слоев, увеличивая кол-во элементов по толщине до 10, остальные параметры оставляем как есть. Общая толщина стенки у нас 0.1 м, и десять слоев получаются с шагом 0.01 м. Отсюда нам для сравнения нужны точки 0.03, 0.06 и 0.09, значения, которые мы использовали для расчета по теории автоматического управления.

Материал стенки по умолчанию у нас сталь простая.

Рисунок 6. Настройка послойной тепловой модели стальной стенки

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

Рисунок 7. Формирования ступенчатого скачака температуры

Запускаем на расчет и получаем вот такой вот график массива температур по толщине.

Рисунок 8. Температуры по толщине стенки.

Cохраняем проект под именем стенка с температурой поверхности.prt

Чтобы вывести два расчета на один график, соберем пакет в SimInTech. Тогда в одном синхронном времении будет считаться и модель стенки в теплогидравлическом коде, и температура по волшебной формуле переходного звена. Массив температур стенки предадим в базу сигналов проекта, откуда его заберем в модель автоматики. Для этого добавим в базу данных сигналов новую категорию с единственным сигналом Т_st

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

Рисунок 9. Настройка пакет для совместного расчета теплогидравлики с автоматикой

В главном скрипте проекта САР.prt добавляем код, который возьмет из базы данных массив, а из массива заберет три точки для сравнения. Поскольку массив у нас начинается с поверхности, то номера в массиве сдвинуты на 1: 0.03 соответствует 4- й элемент массива, 0.06 – 7-й, 0.09 – 10-й.

Далее выполним пересчет в относительную температуру для корректного сравнения.

Рисунок 10. Выборка нужных точек из массива тепловых слоев

Чтобы схема была компактной, я расчет температуры по волшебной формуле убрал в субструктуру, которая на выход выдает нам весовую функцию. Как известно, весовая функция — это производная от переходной. Поэтому на схеме это интегрирующее звено. Интегрируя весовую функцию w(t), мы получаем переходную функцию h(t). В итоге в проекте САР.prt у меня получилась вот такая схема, как на рисунке 1. В зеленой части – расчет из по формуле, в желтой части – вывод на график температуры, рассчитанной в тепловой многослойной модели.

Рисунок 11. Модель сравнения тепловой модели и формулы ТАУ

Запускаем на расчет и видим замечательное совпадение двух расчетов в течение 100 секунд расчета:

Рисунок 12. Сравнение тепловой модели и расчета ТАУ

Графики на белом фоне – результаты расчета по формуле переходной функции, графики на сером – многослойный тепловой расчет.

Если свести расчеты на один график, мы увидим, что отклонения минимальны:

Рисунок 13. Сравненеи двух моделей в относительных температурах

Но самое удивительное, что при послойном тепловом расчете графики температур в относительном времени ровно также собираются в одну кривую!

Рисунок 14. В температуры в многослойной модели в относительном времени

Обратите внимание, в теории использовалось предположение, что наша стенка полубесконечна. Модель тепловая учитывает толщину стенки, а у нас толщина всего 0.1 м = 10 см. Поэтому со временем должны накапливаться отклонения между реальной и теоретической температурой, но согласитесь, 100 секунд расчета и отличное совпадение — это хороший результат. Увеличим время расчета в гидравлической модели до 1000 сек.

Рисунок 15. Сравнение температур в процессе до 1000 сек.

Видим, что в этом случае где-то после 200 секунд начинаются расхождения, причем последняя точка отклоняется больше и начинает отклоняться раньше (синяя линия). И это не удивительно. В теоретической модели должны быть бесконечные слои, и тепловая волна должна все время уходить, отводя энергию следующим слоям, а поскольку у нас нет следующих слоев, происходит нагрев.

Давайте изменим тепловую модель так, чтобы она больше соответствовала бесконечной стенке. Сделаем толщину в 10 раз больше (1 метр), и количество слоев тоже увеличим в 10 раз. Тогда нам не придется пересчитывать положение точек.

В этом случае совпадение многослойного теплового расчета и расчета по волшебной формуле практически идеальное. Судите сами:

Рисунок 16. Сравнение в абсолютных температурах

На одном графике :

Рисунок 17. Сравнение моделей на одно графике

А на графике в относительном времени все кривые теплового расчета выстроились в идеальную линию.

Рисунок 18. Температура в многослойной модели в относительном времени

Это просто фантастика!

Даже на моем ноутбуке время расчета для 100 точек тепловой сетки заметно больше, чем 10, а расчет по волшебной формуле происходит мгновенно!

В заключение хочу заметить, что данная задача с температурой стенки совсем не теоретическая, а, наоборот, очень даже практическая. Именно такие задачи считают при определении термических циклов для теплового оборудования, напряжении на стенки сопла при запуске реактивных двигателей. Сейчас для этого используют сеточные методы и суперкомпьютеры. А раньше использовали только голову и листок миллиметровой бумаги (похоже предки были умнее нас). Я начинаю понимать почем Фон Браун запустил людей на Луну сразу после второй мировой. Лунная программа: 9 лет от старта до посадки на Луне. А спустя 80 лет Илон Маск за 18 лет смог только на орбиту запустить человека. Об этом у меня есть отдельная статья.

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

Рисунок 19. Сравнение точности расчета двух софтов.

Небольшая видео демонстрация сравнения двух моделей:

Модель из видео можно взять здесь...

Но откуда там все-таки взялось число $\pi$ ?!? Есть идеи?

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

Hubs:

If this publication inspired you and you want to support the author, do not hesitate to click on the button