MATLAB è Simulink íà ðóññêîì

https://hub.exponenta.ru/

Femlab

В.Е.Шмелев "Заметки по использованию системы FEMLAB"

В оглавление книги \ К следующему разделу \ К предыдущему разделу

Решение системы телеграфных уравнений относительно напряжения и тока для однородной двухпроводной линии

Постановка задачи

Пусть имеется однородная двухпроводная линия, к началу которой подключен идеальный источник синусоидальной ЭДС, а к концу - резистивная нагрузка Z2 = 0,2 кОм. Исходные данные:

l = 2 км - длина линии;
r0 = 0.2 кОм/км - сопротивление пары проводников, отнесённое к единице длины линии;
L0 = 1 мГн/км - индуктивность на единицу длины линии;
g0 = 0 - проводимость утечки на единицу длины линии (здесь предполагается, что изоляция идеальная);
C0 = 25 нФ/км - ёмкость на единицу длины линии;
- напряжение в начале линии (циклическая частота рад/мкс, частота f = 500 кГц);
Z2 = 0,2 кОм - сопротивление нагрузки в конце линии (линия работает на активную нагрузку).
В момент времени t = 0 к началу линии подключается идеальный источник ЭДС, мгновенное значение которой равно u1 ( t ). В момент подключения источника во всех сечениях линии токи и напряжения были равны нулю (задача с нулевыми начальными условиями).

Электромагнитные процессы в линии описываются системой телеграфных уравнений относительно напряжения и тока:

(1)

Для удобства моделирования в системе FEMLAB систему уравнений (1) запишем в матричной форме:

(2)

Моделирование в GUI-приложении femlab

Запуск femlab приведёт к развёртыванию окна Навигатора моделей с раскрытой закладкой New. Навигатором моделей выберем прикладной режим 1 D / PDE modes/ Coefficient/ Time dependent , нажмём кнопку More .В строку редактирования Dependent variables впишем: u ii. Это означает, что искомая функция u ( x , t )(напряжение) в (1) и (2) в модели femlab будет обозначаться именем u , а функция i( x, t) (ток) - именем ii . Нажатие кнопки OK приведёт к закрытию Навигатора моделей и переводу GUI femlab в одномерный режим моделирования.

Моделирование начинается с режима Draw Mode . В этом режиме прорисовывается геометрия расчётной области. Для создания геометрического объекта, соответствующего линии, выполним команду меню Draw / Specify Geometry , которая разворачивает диалоговое окно, изображённое на рис. 1. Введённые в диалоговое окно числовые значения означают следующее: началом линии является точка с координатой x = -1 км, а концом - точка x = +1 км. Нажатие кнопки OK приведёт к закрытию окна и созданию солидного одномерного объекта.


Рис. 1. Диалоговое окно создания геометрического объекта.

Когда объект создан, можно перейти в режим задания граничных условий Boundary Mode . В этом режиме выполним команду меню Boundary / Boundary Settings . Эта команда разворачивает диалоговое окно ввода и редактирования граничных условий (ГУ). Окно показано на рис. 2.


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

В диалоговом окне, изображённом на рис. 2, точка 1 - начало линии, точка 2 - конец линии. Если в задаче моделирования количество зависимых переменных больше одной, то диалоговое окно редактирования граничных условий состоит из шести закладок, четыре из которых предназначены для ввода матриц ГУ. В нашем случае с помощью этих матриц можно задать режим питания линии и режим нагрузки.

Сначала опишем режим питания линии. Для этого выделим точку 1, тип ГУ - Дирихле. В закладки q и g заходить не обязательно, т.к. там стоят нули по умолчанию. Значения, которые нужно ввести для точки 1, показаны на рис. 3, 4, 5.


Рис. 3. Матрица коэффициентов условия Дирихле в случае задания идеального источника ЭДС


Рис. 4. Матрица правой части ГУ Дирихле, если в данную точку линии включен источник


Рис. 5. Матрица коэффициентов ГУ Дирихле, если в данную точку (конец) линии включен резистор

Матрицу r для точки 2 зададим равной нулю. Следует обратить внимание также на то, что определитель матрицы h в начале и в конце линии равен нулю.

Теперь можно переходить в режим Subdomain Mode и вводить числовые значения параметров линии. Команда меню Subdomain / Subdomain Settings раскрывает диалоговое окно Subdomain Settings , состоящее из десяти закладок (рис. 6). Поскольку уравнения (1), (2) содержат дифференциальные операторы первого порядка, матрицу c примем равной нулю (рис. 6).


Рис. 6. Диалоговое окно ввода коэффициентов PDE

Остальные матрицы введём в соответствии с уравнением (2):
; ; ;; ;.
В закладке Element по умолчанию выбраны Лагранжевы элементы второго порядка.

С помощью команды меню Mesh / Mesh Mode перейдём в режим построения конечноэлементной сетки. Автоматически расчётная область будет разбита на 15 элементов. Это очень грубая сетка для решения поставленной задачи. Командой меню Mesh / Refine Mesh переопределим конечноэлементную сетку. Теперь расчётная область будет разбита на 30 элементов. Это также слишком грубо. Будем выполнять команду Mesh / Refine Mesh до тех пор, пока количество элементов не достигнет 120.

Командой меню Solve / Parameters раскроем диалоговое окно настройки параметров решателя. Раскроем в этом окне закладку Timestepping . В строку редактирования Output times введём 0:0.1:20 (рис. 7) и нажмём OK . Остальные параметры оставим по умолчанию. Командой меню Solve / Solve Problem запустим краевую задачу на выполнение.


Рис. 7. Диалоговое окно настройки параметров решателя (закладка >Timestepping)

Визуализация решения

После завершения решения GUI femlab переходит в режим визуализации. Для одномерных задач возможно три режима визуализации: показ функции решения в фиксированные моменты времени, анимация, показ решения в линейном или точечном сечении расчётной области. Первые два режима настраиваются командой меню Post / Plot Parameters . Эта команда развёртывает диалоговое окно, вид которого показан на рис. 8. Ниспадающее меню Solution at time позволяет выбрать требуемый момент времени моделирования для показа функции решения. На рис. 9 показано распределение напряжения в линии спустя 9 микросекунд с момента подключения источника к линии. Видно, что в этот момент времени волна напряжения ещё не достигла конца линии.


Рис. 8. Диалоговое окно настройки параметров визуализации решения


Рис. 9. Распределение напряжения в линии в момент времени t = 9 мкс

Если нам надо показать распределение тока в линии, то в закладке Curve диалогового окна Plot Parameters в ниспадающем меню Y expression нужно выбрать ii ( ii ). Распределение тока (в миллиамперах) в тот же момент времени показано на рис. 10.


Рис. 10. Распределение тока в линии в момент времени >t = 9 мкс

Значительный интерес представляет также визуализация решения в линейном сечении расчётной области. Пусть для простоты линейное сечение представляет собой всю расчётную область. Данный вид визуализации может осуществляться в двух режимах: 1) отображение в фигуре MATLAB семейства графиков выбранной функции решения в указанные моменты времени, 2) отображение трёхмерного поверхностного графика, одной координатной осью которого является время. Попробуем проиллюстрировать последний из этих режимов.

По команде меню Post / Cross - Section Plot Parameters разворачивается диалоговое окно настройки параметров режимов визуализации в сечениях расчётной области (рис. 11). В закладке General выделим все позиции списка Solution at time . В закладке Line выделим зону № 1 расчётной области (в данной задаче она всегда является выделенной). Установим в этой же закладке флажок Connect lines to surface и нажмём кнопку OK . В фигуре MATLAB будет нарисовано изображение пространственно-временного распределения интересующей нас функции решения. Распределение напряжения в линии показано на рис. 12.


Рис. 11. Диалоговое окно настройки параметров режимов визуализации в сечениях


Рис. 12. Пространственно-временное распределение напряжения в линии

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

Отражение волн напряжения и тока в режиме холостого хода
и короткого замыкания на конце

Чтобы смоделировать режим холостого хода на конце линии, вернёмся опять в режим Boundary Mode. В точке 2 зададим
h22 = 1, остальные коэффициенты матрицы h зададим равными нулю. Повторим расчёт напряжения и тока в линии. На рис. 13 показано распределение напряжения в линии в момент времени t = 13.5 мкс. На рис. 14 - распределение тока в тот же момент времени.


Рис. 13. Распределение напряжения в линии в момент времени t = 13.5 мкс


Рис. 14. Распределение тока в линии в момент времени t = 13.5 мкс

Пространственно-временное распределение напряжения в линии показано на рис. 15, а тока - на рис. 16. Эти трёхмерные графики можно вывести в разные фигуры MATLAB , если перед выполнением команды построения графика воспользоваться функцией figure . Это бывает полезно при необходимости сопоставления графиков. Следует обратить внимание на то, что в любой момент времени ток в конце линии равен нулю.


Рис. 15. Пространственно-временное распределение напряжения в режиме холостого хода на конце


Рис. 16. Пространственно-временное распределение тока в режиме холостого хода на конце

Чтобы смоделировать режим короткого замыкания на конце линии, вернёмся опять в режим Boundary Mode . В точке 2 зададим h11 = 1, остальные коэффициенты матрицы h зададим равными нулю. Повторим расчёт напряжения и тока в линии. Пространственно-временное распределение напряжения в линии показано на рис. 17, а тока - на рис. 18. Следует обратить внимание на то, что в любой момент времени напряжение в конце линии равно нулю.


Рис. 17. Пространственно-временное распределение напряжения в режиме короткого замыкания на конце


Рис. 18. Пространственно-временное распределение тока в режиме короткого замыкания на конце

Реактивная нагрузка на конце линии

Опять перейдём в режим Boundary Mode и выполним команду меню Boundary/ Boundary Settings. Чтобы смоделировать индуктивную нагрузку в конце линии, нужно в диалоговом окне Boundary Settings выбрать точку 2, обнулить для этой точки матрицы h и r , заполнить первый столбец матрицы q единицами. Остальные элементы матриц q и g должны быть равны нулю. Далее нужно открыть закладку Weak диалогового окна Boundary Settings и вписать в строку редактирования dweak выражения    -u_test*0.3*ii_time   -ii_test*0.3*ii_time   Эти выражения соответствуют индуктивности нагрузки, равной 0.3 мГн.

После ввода граничных условий, соответствующих индуктивной нагрузке, можно запускать вычисления. Для проверки правильности моделирования индуктивной нагрузки выведем в две фигуры MATLAB графики напряжения и тока в конце линии. Развернём диалоговое окно Cross - Section Plot Parameters и раскроем закладку Point . В этой закладке включим радиокнопку Point plot , в рамке Select point via включим радиокнопку Vertices , в списке выберем точку 2. В строку редактирования Point plot впишем имя u . В фигуру № 1 будет выведен график напряжения в конце линии. Аналогичную операцию проведём для тока (имя зависимой переменной ii ), предварительно выполнив оператор MATLAB figure (2) . Графики напряжения и тока в конце линии показаны на рис. 19, 20. По графикам видно, что синусоидальная составляющая тока отстаёт от напряжения на четверть периода, т.е. индуктивная нагрузка смоделирована правильно.


Рис. 19. График напряжения в конце линии


Рис. 20. График тока в конце линии

Пространственно-временное распределение напряжения и тока в линии показано на рис. 21, 22. Как видно, от индуктивной нагрузки волны напряжения и тока отражаются почти так же, как и от разомкнутого или короткозамкнутого конца. Кроме того, индуктивная нагрузка обуславливает наличие непериодической составляющей тока в линии и в нагрузке. Особенно хорошо это видно на рис. 20.


Рис. 21. Пространственно-временное распределение напряжения


Рис. 22. Пространственно-временное распределение тока

Аналогичным способом можно смоделировать ёмкостную нагрузку в конце линии. Однако это практически не удаётся из-за абсолютной неустойчивости получаемой вычислительной модели. Я предлагаю другой способ моделирования нагрузки с сосредоточенными параметрами.

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

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

Переключимся в режим Draw Mode и в конец линии добавим солидный геометрический объект, начинающийся в точке x =1 и заканчивающийся в точке x =1+1/60. Теперь можно переключаться в режим Boundary Mode и вводить граничные условия. В начале линии граничные условия изменять не будем. Концом обобщённой линии теперь станет точка 3, для которой зададим условие короткого замыкания.

Для проверки правильности идеи метода смоделируем резистивную нагрузку исходной однородной линии. Роль резистора будет выполнять зона № 2 расчётной области (добавленный кусок линии длиной 1/60 км). В зоне № 2 зададим следующие значения матриц свойств "линии" (сопротивление нагрузки равно 0,2 кОм):

; ;;; ; .

Для экономии вычислительных ресурсов конечноэлементную сетку желательно сгенерировать так, чтобы в зоне № 2 был построен только один элемент, причём сетка в исходной линии должна остаться такой же, как и в предыдущих расчётах. Для таких случаев (и не только) есть команда меню Mesh / Parameters , позволяющая устанавливать значения параметров режима генерации сетки. Выполним эту команду и развернём диалоговое окно ввода параметров сетки Mesh Parameters . В строку редактирования Max element size , general впишем выражение 1/60. Нажмём кнопку OK и переключимся в режим Mesh Mode . В результате с первого раза будет построена равномерная сетка, состоящая из 121 элемента, один из которых будет принадлежать зоне № 2. Параметры решателя оставим такими же, как и в предыдущих расчётах.

Результат расчёта будет таким же, как и на рис. 12. Совпадение результата свидетельствует о правильности идеи метода пространственного распределения сосредоточенной нагрузки.

Теперь можно попробовать смоделировать индуктивную нагрузку. В зоне № 2 зададим следующие значения матриц свойств "линии" (индуктивность нагрузки равна 0,3 мГн):

;;;;;.

На рис. 23, 24 представлены графики напряжения и тока в конце линии. Сравнение этих графиков с рис. 19, 20 показывает неплохое совпадение результатов, причём методом пространственного распределения все расчёты выполняются более устойчиво и быстро. На рис. 25, 26 представлено пространственно-временное распределение напряжения и тока в моделируемой системе.


Рис. 23. Напряжение в конце линии при индуктивной нагрузке


Рис. 24. Ток в конце линии при индуктивной нагрузке


Рис. 25. Пространственно-временное распределение напряжения при индуктивной нагрузке


Рис. 26. Пространственно-временное распределение тока при индуктивной нагрузке

Теперь можно попробовать смоделировать ёмкостную нагрузку. Перейдём в режим Boundary Mode и в точке 3 зададим условие холостого хода. Перейдём в режим Subdomain Mode и в зоне № 2 зададим следующие значения матриц свойств "линии" (ёмкость нагрузки равна 0,3 нФ):

; ;;; ; .

На рис. 27, 28 представлены графики напряжения и тока в конце линии. Напряжение отстаёт от тока на четверть периода. Это свидетельствует о правильности моделирования ёмкостной нагрузки. Численная устойчивость немножко хуже, чем при моделировании индуктивной нагрузки. На рис. 29, 30 представлено пространственно-временное распределение напряжения и тока в моделируемой системе.


Рис. 27. Напряжение в конце линии при ёмкостной нагрузке


Рис. 28. . Ток в конце линии при ёмкостной нагрузке


Рис. 29. Пространственно-временное распределение напряжения при ёмкостной нагрузке


Рис. 29. Пространственно-временное распределение тока при ёмкостной нагрузке

В оглавление книги \ К следующему разделу \ К предыдущему разделу


Поиск по сайту:


Система Orphus