Python. Расчёт переходных процессов

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

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

Изменения этих величин происходят скачкообразно только при переходе от одного участка к другому. Это чем-то напоминает работу АЦП (Аналого-Цифрового Преобразователя).

При таком методе расчета неизбежно будут возникать ошибки квантования. Но можно заранее выбрать шаг вычислений настолько мелкий, что эти ошибки никак не скажутся на результатах расчёта. Но самое главное — неберущийся интеграл будет всё-таки вычислен, а переходной процес в нелинейной системе — будет всё-таки просчитан (смоделирован).

Сейчас для этих задач используются разные специализированные программы такие как, например, LabView, MicroCAP и другие. Я с этими  программами не знаком, я не могу даже их перечислить.

Молодёжь их во всю юзает. Это удобно — накидал на экране схему, определил параметры, задал начальные условия, нажал на кнопку и получил результат. Что там и как посчиталось — не важно! А результат — вот он! Быстро. Удобно. Достоверно.

А тот аспект, что люди не знают, как это получается — это не важно. Многие ли нынешние водители могут с первого раза показать на датчик массового расхода воздуха в своём автомобиле. Не удивлюсь, если добрая половина водителей об этом вообще что-то слышала. И тем не менее — ездят все более-менее успешно.

Так что незнание своего инструмента — это не критерий того, что человек не имеет право им пользоваться. Имеет. Но квалификация такого юзвера, разумеется, будет ниже, чем того, кто понимает как работает инструмент.

На днях мне довелось порешать одну задачку. Собственно, проблем особых не было, нужно было просто просчитать переходной процесс в одном изделии. Не буду раскрывать всех тайн: умный и так поймёт, а несведущему и знать не положено.

Имеется кусок обычного провода типа ПВ сечением 4 мм². Длина провода равна 40 м. Причем провод расположен в пространстве (например на полу или на земле) в виде одного витка диаметром 10 м, и примерно 5-метровым «хвостом», которым этот виток подключается к аппаратуре.

Кроме того имеется батарея конденсаторов, общей ёмкостью 5000 мкФ, заряженная до напряжения 300 В. Уже интересно!

И в какой-то момент времени батарея кондёров подключается к этому витку. Что-то должно произойти весьма занимательное…

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

Ec = C * U² /2

Для сравнения. Энергия выстрела из пистолета составляет от 300 до 600 Дж, в зависимости от типа. Разрядить эту батарею конденсаторов — это все равно, что бахнуть из тяжелого пистолета. Ого-го!

Схема устройства предельно простая:

rlc

Здесь кусок провода и другие омические компоненты схемы представлены в виде сопротивления R. Остальное, думаю итак понятно.

Первое, что приходит в голову — «Блин, это же резонансный контур!» Было бы неплохо определить его резонансную частоту.

Нет проблем!

Емкость известна, а индуктивность сейчас найдём по формуле

L = 0.0002 π D (ln(8D/d) — 1.75).

Загружаем Питон и вводим следующие команды:

>>> import math
>>> D = 10000
>>> s = 4
>>> d = math.sqrt(4 * s /math.pi)
>>> L = 0.0002 * math.pi * D * (math.log(8 * D / d) - 1.75)
>>> L
54.8261481920958

Здесь D — диаметр витка (мм), s — сечение провода (мм²), d — диаметр провода.

Индуктивность равна 55 мкГн. Индуктивность «хвоста» не учитываем, так как токи в его проводах текут в противоположные стороны, а физическая площадь между этими проводами ничтожно мала.

Продолжаем. По известной формуле

f = 1 / (2 π √LC)

находим резонансную частоту:

>>> C = 5000E-6
>>> f = 1 / (2 * math.pi * math.sqrt(1E-6 * L * C))
>>> f
303.9773759368137

300 герц — ого! Оказывается резонансная частота контура лежит в диапазоне звуковых волн. Значит мы можем услышать как будут «петь» всякие гаечки-заклёпки.

Теперь найдем сопротивление провода. Известно, что удельное сопротивление меди равно 0.0175 Ом * мм² / м. Тогда сопротивление провода будет равно:

>>> ro = 0.0175
>>> l = 10
>>> R = ro * l / s
>>> R
0.043750000000000004

То есть примерно 0.044 ома.

Но у нас в реальной установке присутствуют переходные разъёмы, другие провода, а так же сопротивление контактов самого выключателя. (Вот тут я дико смухлевал, так как вместо механического выключателя используется мощный MOSFET.)

Я предлагаю взять омическое сопротивление равное 0.2 Ом. Большой ошибки не будет.

На этом предварительные расчёты закончим и перейдём к созданию программы.

Как я уже упоминал, на коротком шаге вычислений будем считать, что напряжение на конденсаторах не будет изменяться, так же не будет изменяться ток через виток провода. С другой стороны, полагаясь на то, что напряжение на конденсаторе неизменно (существенно не изменяется за время dt), мы можем определить на сколько измениться ток в витке:

dI = Uc / L * dt

А так же полагаясь на постоянство тока в витке можем определить на сколько измениться напряжение на конденсаторе:

dUc = I / C * dt

Только тут есть один специфический момент, который обязательно нужно учитывать. Когда мы определяем степень изменения тока в витке, мы должны учитывать падение напряжения на сопротивлении провода. Иначе говоря, правильная формула для расчета должна быть такая:

dI = (Uc — Ur) / L * dt

, где Ur — падение напряжения на сопротивлении провода, возникающее в следствие протекания по нему тока I.

Что касается вычисления изменения напряжения на конденсаторе, то учитывать сопротивление в этом случае не надо — ток во всей цепи один и тот же!

Нам осталось определиться с величиной шага времени (dt) и конечным временем вычислений.

Сейчас я опять немного смухлюю. Я уже имею практический опыт расчёта этой цепи и уже знаю какие значения нужны. А для всех остальных я предлагаю «пойти логическим путём» (с).

Мы знаем, что резонансная частота контура составляет 300 Гц. Это значит, что период колебаний равен примерно 3.3 мс. Давайте возьмём конечное время вычислений равным 5 мс.

Что касается шага вычислений, то я предлагаю взять значение 1 мкс. Это достаточно мелко, для того чтобы расчёты не «поплыли» из-за накопления ошибок квантования. А с другой стороны, мы получим количество шагов вычислений равное 5000 (5 мс / 1 мкс); с практической точки зрения, это вполне приемлемое количество.

Вот готовая программа:

#!/usr/bin/env python3
# coding:utf-8

'''
rlc.py

Программа для расчёт переходных процессов в последовательной RLC-цепи
'''

def lc():
    L = 55E-6       # Индуктивность петли (витка), Гн
    C = 5000E-6     # Ёмкость батареи электролитических конденсаторов, Ф
    R = 0.2         # Омическое сопротивление проволоки, Ом
    U0 = 300        # Начальное напряжение на конденсаторах, В
    I0 = 0          # Начальное значение тока, А
    t = 0           # Начальное время, с
    dt = 1.0e-6     # Шаг времени, с
    Tmax = 5.0E-3   # Конечное время вычисления, с

    # Начальные значения
    u = U0
    i = I0

    while (t < Tmax):
        uR = i * R
        di = (u - uR) / L * dt
        du = i / C * dt
        i += di
        u -= du
        t += dt

        # Печать результатов расчёта (время, напряжение, ток)
        print("{0:f} {1:f} {2:f}".format(t, u, i))

if __name__ == "__main__":
    lc()

В результате её работы мы получим три колонки чисел:

0.000001 300.000000 5.454545
0.000002 299.998909 10.889256
0.000003 299.996731 16.304185
0.000004 299.993470 21.699383
0.000005 299.989131 27.074902
0.000006 299.983716 32.430796
0.000007 299.977229 37.767115
0.000008 299.969676 43.083912
0.000009 299.961059 48.381237
...
0.004997 -0.004809 0.277669
0.004998 -0.004865 0.276571
0.004999 -0.004920 0.275477
0.005000 -0.004975 0.274386

В первой колонке выведено время (в секундах), во второй – напряжение на конденсаторе (в
вольтах), в третьей – ток в петле (в амперах).

Запускать программу лучше командой:

$ rlc.py > rlc.dat

что позволит нам получить данные не на экран терминала, а в файл rlc.dat.

— Ичё щас с этим файлом делать? — зададут вопрос несообразительные студенты.

Да, ничего! Запускаем gnuplot и вводим следующую команду:

gnuplot> plot "rlc.dat" u 1:2 w l t "U(t)", "" u 1:3 w l t "I(t)"

В результате на экране получаем график переходного процесса:

rlc

Как мы можем видеть, переходной процесс вопреки нашим ожиданиям апериодический, а не колебательный. Заканчивается переходной процесс в течение каких-то 3-4 мс. Амплитуда тока превышает значения 1000 ампер. Ого! Как тяжело работать бедному MOSFET-у.

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

Ну и самое главное — результат достигнут простыми средствами.

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

gnuplot> set terminal png
gnuplot> set output "RLC.png"

, это перенаправит вывод gnuplot в файл RLC.png.

Надеюсь, вам было не только интересно почитать, но вы получили также полезные знания.

Реклама

3 responses to “Python. Расчёт переходных процессов

  1. Пользовался Microcap и Proteus для расчетов в windows. Для расчета моделей автоматических систем юзал VisSim. А таким образом, я понимаю, можно свести все к одному ПО? Но необходимы знания Python и прочные знания по математике. Нужно попробовать, спасибо.

    • Я не совсем понял Ваш вопрос.

      Конечно, для моделирования процессов, происходящих в электрических схемах, безусловно, нужны знания. Это даже не обсуждается.

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

      Что касается программирования, ну такие задачи можно решать почти на любом языке высокого уровня — на Питоне, на Бэсике, на Паскале, на Си, на С++… особых (глубоких) знаний языка тоже не требуется.

      Сложнее, наверно, со знаниями основ электротехники (ТОЭ) и электроники. Ну, тут уж ничего не поделаешь. Всё-таки это профильные знания, а не просто расчеты чего-нибудь.

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

      Я, вот, не знаю Microcap и Виндовса у меня тоже нет. Поэтому я пошел таким путём.

  2. Для меня удобство в том, что ненужно вороха программ для выполнения схожих задач. Линукс как раз и подкупает универсальностью, только разобраться немного нужно. gnuplot отличная вещь. Как то снимал показания температур с нагревателя, вручную. После, для создания модели, данные переносил в ексель строил график. Моделировал, сравнивал графики. Это уже 2 немаленьких софтины + ручной труд и расчеты на бумажке. Сейчас бы решил все намного проще. И спасибо, что то не додумался, можно же использовать любой язык программирования.

Добавить комментарий

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход / Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход / Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход / Изменить )

Google+ photo

Для комментария используется ваша учётная запись Google+. Выход / Изменить )

Connecting to %s