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

Попробуем провести анализ электрических схем по переменному току с использованием системы символьных вычислений 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²+ω₀²)/ω₀, ω₀ = 2•π•f₀

Запишем переходную характеристику как результат умножения входного сигнала на передаточную функцию, в данном случае выполняется подстановка вместо 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];