Символический анализ цепей переменного тока

В установившемся режиме

Попробуем провести анализ электрических схем по переменному току с использованием системы символьных вычислений Maxima. Представим себе произвольную схему, после чего:

• добавляем источник синусоидальной ЭДС в общем виде A•sin(2•π•f₀•t+ψ₀)

• имеем в виду что 2•π•f₀•t+ψ₀ — это фаза, а ψ₀ — начальная фаза

• что f₀ — циклическая частота [Гц], ω₀ — круговая [рад•с⁻¹],

• выбираем необходимые переменные состояния, обозначаем токи и напряжения

• рисуем контуры так, чтобы на каждом элементе подключённого к двум ветвям контуры проходили хотя бы пару раз, обозначаем узлы с охватом всех ветвей

Символический анализ

Используем систему компьютерной алгебры Maxima. Для этого составляем первую систему уравнений. По умолчанию при составлении системы уравнений предполагается, что каждое значение в списке равно нулю. Знак равенства можно не использовать. a=b записать как
a-b=0 или просто a-b. Таким образом, можно обойтись без знака равенства, полагая результатом каждого выражения в книге ноль (может быть, поэтому, на самом деле похожи
— и =).

eq:[ /*Правило Кирхгофа для напряжений контуров*/ -E1 + iR1*R1+iL*p*L+iC/(p*C), /*Уравнение по первому контуру I*/ +iL*p*L-iR2*R2, /*Уравнение по второму контуру II*/ -iC/(p*C)+iR3*R3, /*Уравнение по третьему контуру III*/ /*Правило Кирхгофа для токов узлов*/      iR1-iR2-iL, /*Уравнение по первому узлу (1)*/     iL-iC-iR3, /*Уравнение по первому узлу (2)*/    /*Вычисляемая переменная состояния - выход*/    Un=iR3*R3   ];

далее производится решение и выделение необходимой переменной состояния Un. Небольшая заметка

• Если в Maxima имеется выражение V1:[a=1,b=x,c=…], то для того чтобы это подставить в некоторое выражение V2:a+b можно использовать конструкцию с оператором «запятая» V3:V2,V1 которая эквивалентна subst(V1,V2). В основном subst используется внутри функций, циклов и условий

/*решаем систему уравнений eq относительно переменных  состояния с выбором Un в качестве ответа*/ sln:Un,solve(eq,[iR1,iR2,iR3,iC,iL,Un]);

Получается передаточная функция второго порядка. Следует уточнить, что в данном случае используется термин «передаточная функция» (как отношение выхода ко входу), несмотря на то что в выражении содержится входной сигнал E1. Для упрощения можно полагать этот сигнал как произвольный множитель, входящий в соотношение, в частности, равный единице или нормированный к номинальной величине. Переходная характеристика же содержит конкретный вид сигнала E1, в данном случае синусоидальную величину. Передаточная функция используется для частотного анализа, переходная характеристика — для последующего анализа во временной форме и для нахождения значения в установившемся режиме.

Запишем входные сигналы во времени для варианта постоянного (единичная ступенька) и переменного тока (синус, начинаемый с нуля), а также найдём их изображения по Лапласу:

/*Для постоянного тока - постоянный уровень U*/ V1:U; /*Для переменного тока - синусоидальная величина с амплитудой A*/ V2:A*sin(2*%pi*f_0*t+ψ_0); /*Преобразование Лапласа от постоянной*/ V3:laplace(V1,t,p); /*Преобразование Лапласа от синуса*/ V4:laplace(V2,t,p);

В результате будет иметься следующее:

Для нахождения значения в установившемся режиме необходимо использовать следующие пределы с «магическим» обведённым коэффициентом для переменного тока, который работает как показано далее.

Код для получения формулы на рисунке выше

/*Магический коэффициент*/ V5:(p^2+4*%pi^2*f_0^2)/(2*%pi*f_0); /*Произведение при заданной форме входного сигнала*/ V6:V5*V4;

Итак, необходимым для анализа коэффициентом является выражение, которое умножается на произведение заданной формы сигнала E1(p) в операторной форме на передаточную функцию системы W(p). Также коэффициент можно записать как (p²+ω₀²)/p

Запишем переходную характеристику как результат умножения входного сигнала на передаточную функцию, в данном случае выполняется подстановка вместо E1:

/*Переходная характеристика*/ V7:sln,E1=V4;

Результатом будет выражение V7

Далее переходная характеристика умножается на «Магический коэффициент»:

/*Переходная характеристика с коэффициентом*/ V8:factor(V7*V5);

Получается следующее выражение V8, равное

Как видно, «магический коэффициент» предназначен для устранения сингулярности в знаменателе, образованной преобразованием Лапласа от синусоидального входного сигнала. Так как если в p² + (2•π•f₀)² подставить p = j•2•π•f₀, то получится (j•2•π•f₀)² +(2•π•f₀)² = -(2•π•f₀)² +(2•π•f₀)² = 0 деление на ноль. Далее можно подставить значение p для перехода в комплексно-частотную область для заданной точки.

/*Значение оператора p в заданной точке и переход в КЧХ*/ V9:p=%i*2*%pi*f_0; /*Значение выходного сигнла в заданной комплексно-частотной точке*/ V10:V8,V9;

Получается значение в установившемся режиме как комплексно-частотная точка

Итак, общая схема для нахождения значения в установившемся режиме представлена на рисунке

Для нахождения амплитуды и фазы выходного сигнала достаточно взять модуль и аргумент. Перед этим, также задав предположение о положительном значении величин, чтобы система символических вычислений могла оптимизировать корни. √x² = |x|, но √x² = x если x≥0

/*Предположение о неотрицательности параметров*/ assume(A>0,R1>0,R2>0,R3>0,f_0>0,C>0,L>0); V10A:cabs(V10); /*Амплитуда выходного сигнала в заданной частотной точке*/ V10P:carg(V10); /*Фаза выходного сигнала в заданной частотной точке*/

В итоге получается следующий ответ

Подставим численные параметры в данные выражения и рассчитаем получившуюся амплитуду и фазу выходного синусоидального сигнала, например, на частоте 2 кГц f₀=2e3 и фазой входного сигнала -40 эл. град. Кратко входной сигнал можно записать как 1∠-40°

Для перевода в радианы умножим на π/180, чтобы выразить %pi значением используется команда float. Следует отметить что %pi в Maxima выглядит как π, но при этом имеет численный атрибут.

/*Задание численных параметров*/ Num:[R1=1,R2=100,R3=30,L=1e-3, C=10e-6,ψ_0=float(-40/180*%pi),A=1,f_0=2e3]; 

/*Получение значения амплитуды*/ V11A:float(V10A),Num; /*Получение значения фазы*/ V11P:float(V10P),Num; /*из радиан в градусы*/ float((V11P)*180/%pi);

В итоге получаются следующие значения: амплитуда (выражение V11A) 1.2736 В, величина в радианах составляет (выражение V11P) -3.06591, что соответствует в градусах -175.66°. Таким образом, выходное напряжение можно записать как 1.27∠-176°, что и является искомым ответом. Данный метод можно использовать для любой схемы или произвольной передаточной функции.

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

Построим теперь те же самые характеристики в Maxima.

/*Подставить численные значения в переходную характеристику*/ V12:V7,Num; ratprint:false;/*чтобы избежать предупреждение о преобразовании в рациональные числа*/ /*Найти обратное преобразование Лапласа от переходной характеристики*/ V13:float(ilt(V12,p,t)); /*полюса передаточной функции - корни знаменателя*/ rectform(float(solve(denom(V12),p)));

Выводом будут следующие выражения. Обратить внимание что 1/2.718…^t это exp(-t)

Важно отметить, что если если p[i] — произвольный корень знаменателя, то

• 1/Re(p[i]) — постоянная времени составляющей переходного процесса (секунда) или величина обратная действительной части корня

• Im(p[i]) — круговая частота переходного процесса (радиан•с⁻¹) или величина мнимой части корня, а Im(p[i])/(2π) — циклическая частота составляющей переходного процесса (Гц)

/*Подстановка численных значений во входной сигнал*/ V14:V2,Num; /*Подстановка полученных значений амплитуды и начальной фазы в выходной сигнал*/ V15:V11A*sin(2*%pi*f_0*t+V11P),Num;

В итоге имеем два сигнала — входной 1∠-0.698 (радиан) и выходной 1.274∠-3.066, с подставленными и найденными ранее численными параметрами
1•sin(4000.0•π•t-0.6981317007977318) — вход, 1.2736•sin(4000.0•π•t-3.06591) — выход в установившемся режиме

/*построить график*/ wxplot2d([V13,V14,V15],[t,0,0.003],[legend,false]),wxplot_size=[480,320];

Далее построим частотные характеристики. Для этого необходимо подставить p = j•2•π•f в передаточную функцию для перехода к комплексно-частотной характеристики, содержащей мнимую и действительную части от варьируемой частоты.

/*Переход в комплексно-частотную область для любой частотной точки f*/ V16:p=%i*2*%pi*f; /*Подстановка в подстановку: задаём ЭДС формально равное единице  для получения передаточной функции*/; V17:sln,[E1=A,V16]; /*Подстановка численных параметров*/ V18:V17,Num; /*Модуль */ V18A:float(cabs(V18)); /*Аргумент*/ V18P:float((carg(V18))/%pi*180); /*Построение графика АЧХ*/ wxplot2d(V18A,[f,1,100e3],[legend,false],[logx,true],[logy,true]),wxplot_size=[480,320]; /*Построение графика ФЧХ*/ wxplot2d(V18P,[f,1,100e3],[legend,false],[logx,true]),wxplot_size=[480,320];

В результате получаются следующие графики, чуть отличающиеся от моделирования в SPICE. Обратить внимание на количество точек «перегиба» и производные изменения. По всей видимости, численная модель имеет некоторые особенности, что требует дальнейших исследований.

Напоследок чтобы проверить схему можно рассмотреть установившийся режим для передаточной функции на постоянном токе, подставляя p=0 в выражение sln командой factor(sln),p=0 можно получить E1•R3/(R3+R1) — правило делителя, при этом индуктивность становится закороткой а ёмкость уходит в разрыв.
Общий скрипт

kill(all); eq:[ /*Правило Кирхгофа для напряжений контуров*/ -E1 + iR1*R1+iL*p*L+iC/(p*C), /*Уравнение по первому контуру I*/ +iL*p*L-iR2*R2, /*Уравнение по второму контуру II*/ -iC/(p*C)+iR3*R3, /*Уравнение по третьему контуру III*/ /*Правило Кирхгофа для токов узлов*/      iR1-iR2-iL, /*Уравнение по первому узлу (1)*/     iL-iC-iR3, /*Уравнение по первому узлу (2)*/     Un=iR3*R3 /*Вычисляемая переменная состояния - выход*/    ]; /*решаем систему уравнений eq относительно выбранных переменных состояния*/ sln:Un,solve(eq,[iR1,iR2,iR3,iC,iL,Un]); /*Для постоянного тока - постоянный уровень U*/ V1:U; /*Для переменного тока - синусоидальная величина с амплитудой A*/ V2:A*sin(2*%pi*f_0*t+ψ_0); /*Преобразование Лапласа от постоянной*/; V3:laplace(V1,t,p); /*Преобразование Лапласа от синуса*/ V4:laplace(V2,t,p); /*Магический коэффициент*/ V5:(p^2+4*%pi^2*f_0^2)/(2*%pi*f_0); /*Произведение при заданной форме входного сигнала*/ V6:V5*V4; /*Переходная характеристика*/ V7:sln,E1=V4; /*Переходная характеристика с коэффициентом*/ V8:factor(V7*V5); /*Значение оператора p в заданной точке и переход в КЧХ*/ V9:p=%i*2*%pi*f_0; /*Значение выходного сигнла в заданной комплексно-частотной точке*/ V10:V8,V9; /*Предположение о неотрицательности параметров*/ assume(A>0,R1>0,R2>0,R3>0,f_0>0,C>0,L>0); V10A:trigsimp(cabs(V10)); /*Амплитуда выходного сигнала в заданной частотной точке*/ V10P:carg(V10); /*Фаза выходного сигнала в заданной частотной точке*/ /*Задание численных параметров*/ Num:[R1=1,R2=100,R3=30,L=1e-3, C=10e-6,ψ_0=float(-40/180*%pi),A=1,f_0=2e3]; /*Получение значения амплитуды*/ V11A:float(V10A),Num; /*Получение значения фазы*/ V11P:float(V10P),Num; /*из радиан в градусы*/ float((V11P)*180/%pi); /*Подставить численные значения в переходную характеристику*/ V12:V7,Num; /*Найти обратное преобразование Лапласа от переходной характеристики*/ ratprint:false;/*чтобы избежать предупреждение о преобразовании в рациональные числа*/ V13:float(ilt(V12,p,t)); /*полюса передаточной функции - корни знаменателя*/ rectform(float(solve(denom(V12),p))); /*Подстановка численных значений во входной сигнал*/ V14:V2,Num; /*Подстановка полученных значений амплитуды и начальной фазы в выходной сигнал*/ V15:V11A*sin(2*%pi*f_0*t+V11P),Num; wxplot2d([V13,V14,V15],[t,0,0.003],[legend,false]),wxplot_size=[480,320]; /*Переход в комплексно-частотную область для любой частотной точки f*/ V16:p=%i*2*%pi*f; /*Подстановка в подстановку: задаём ЭДС формально равное единице  для получения передаточной функции*/; V17:sln,[E1=A,V16]; /*Подстановка численных параметров*/ V18:V17,Num; /*Модуль */ V18A:float(cabs(V18)); /*Аргумент*/ V18P:float((carg(V18))/%pi*180); /*Построение графика АЧХ*/ wxplot2d(V18A,[f,1,100e3],[legend,false],[logx,true],[logy,true]),wxplot_size=[480,320]; /*Построение графика ФЧХ*/ wxplot2d(V18P,[f,1,100e3],[legend,false],[logx,true]),wxplot_size=[480,320];