Analog simulators, like LTSpice, are valuable tool when doing mostly analog design.
Sometimes the requirement might be opposite, most of the design is digital, but some analog behavior is needed in the simulation test bench.
This circuit, just an example, has known voltage feed (V(net_a)from voltage source (V1), and our interest is to know the voltage output (V(net_b)). From the simulation results, we can see that the output waveform is definitely not trivial (voltage of net_b and current of C1 showed).
One way to do this kind of simulation, is to do it in small time steps, and assume that all voltages and currents behave linearly between the steps.
The problem though, is that huge number of steps are needed, and it consumes a lot of CPU time to get the results.
This tutorial shows analytical way of solving the same problem.
To properly analyze the circuit, we must fix our convention of voltage polarities across elements (voltage source, resistors, inductor and capacitor) and current directions through elements. They can be arbitrary, as long as current arrows point from + to -.
Now we need equations describing this circuit (remember Kirchhoff’s law?).
First the voltage loops (three of them).
The voltages are divided to time dependent (V_element(t)) and fixed (V_V1). The fixed voltage is not fixed for the whole simulation, but piecewise fixed. I.e. It can change, but stays fixed (horizontal in waveform view) between change points.
First voltage loop:
In equation (eq_loop_1) this voltage loop is:
V_V1-V_L1(t)-V_R2(t)=0
Basically this says, starting from ground (green bar), and going (any route) back to starting point (ground), and adding all voltage differences en route, the sum of voltages must be zero.
V_V1 stands for (piecewise fixed) voltage of voltage source, V_L1(t) stands for time dependent voltage of L1, etc.
When the green arrow goes from minus to plus (yellow voltage markers), the voltage is added in the equation, and when the green arrow goes from plus to minus, the voltage is subtracted in the equation.
Second voltage loop:
The equation (eq_loop_2) for this loop is:
V_L1(t)-V_C1(t)=0
Third voltage loop:
The equation (eq_loop_3) for this loop is:
V_C1(t)-V_R1(t)=0
Now we create current equations. For each net (continuous wire), the sum of all currents coming to the net from all directions must be zero.
When the white arrow is same direction with orange arrow, the current is added in the equation.
When the white arrow is opposite direction with orange arrow, the current is subtracted in the equation.
First net:
I_V1(t) (subtractive), I_L1 (subtractive), I_C1 (subtractive), I_R1 (subtractive)
The equation (eq_net_a) for this net is:
-I_V1(t)-I_L1(t)-I_C1(t)-I_R1(t)=0
Second net:
The equation (eq_net_b) for this net is:
-I_R2(t)+I_L1(t)+I_C1(t)+I_R1(t)=0
Now we need to define the time dependent equations for each element in the circuit:
Resistors are easy. The Ohm’s law states that V=R*I.
Resistor R1 equation (eq_R1) is:
V_R1(t)=R_R1*I_R1(t)
The R_R1 is the resistance of R1. It is not time dependent, but can be changed during the simulation same way the V1 voltage can be changed.
Resistor R2 equation (eq_R2) is:
V_R2(t)=R_R2*I_R2(t)
For the inductors, the equation is differential. It can be expressed either in integral form, or in derivative form. We use the latter.
Inductor L1 equation (eq_L1) is:
Capacitors are also expressed with differential equations.
Capacitor C1 equation (eq_C1) is:
Now we have 9 equations and 9 time dependent variables in the equations.
Solving the group of equations is done by utilizing Maxima.
Start Maxima (e.g. wxMaxima) in your computer (unfortunately online Maxima, e.g. this, cannot be used at this phase) and paste the following code to it (don’t panic, it probably gives an error at this point). You need to press CTRL-Enter to instruct Maxima to execute the code you entered.
eq_loop_1: V_V1-V_L1(t)-V_R2(t)=0;
eq_loop_2: V_L1(t)-V_C1(t)=0;
eq_loop_3: V_C1(t)-V_R1(t)=0;
eq_net_a: -I_V1(t)-I_L1(t)-I_C1(t)-I_R1(t)=0;
eq_net_b: -I_R2(t)+I_L1(t)+I_C1(t)+I_R1(t)=0;
eq_R1: V_R1(t)=R_R1*I_R1(t);
eq_R2: V_R2(t)=R_R2*I_R2(t);
eq_C1: I_C1(t)=C_C1*'diff(V_C1(t), t);
eq_L1: V_L1(t)=L_L1*'diff(I_L1(t), t);
assume(R_R1>0);
assume(R_R2>0);
assume(C_C1>0);
assume(L_L1>0);
eq_set: [eq_loop_1, eq_loop_2, eq_loop_3, eq_net_a, eq_net_b, eq_R1, eq_R2, eq_C1, eq_L1];
x_set: [I_V1(t), V_R1(t), I_R1(t), V_C1(t), I_C1(t), V_L1(t), I_L1(t), V_R2(t), I_R2(T)];
desolve(eq_set, x_set);The lines starting with eq_ define the 9 equations we formulated earlier. The 9 unknown time dependent variables (ending with (t)) are embedded in the equations.
The lines with assume() are telling Maxima the allowed ranges of different (non time dependent) variables. We know that the values of elements (R, L, C) need to be non-negative (to be physically consistent). In this case we stick to positive values. This might not be necessary, but it gives Maxima some insight which solutions (if many of them are found) are actually valid.
The line starting with eq_set: is making a collection of all equations together.
The line starting with x_set: is making a collection of all unknowns together.
desolve() is Maxima’s procedure to solve group of differential equations using Laplace transform.
But for some reason (unknown to me) this procedure doesn’t work for this set of equations.
So we need to go closer to basics, let’s solve the S-domain equations (equation gone through Laplace transform is said to be in S-domain), and try getting inverse Laplace transform to get back to time domain (equations where t is describing time is said to be in time-domain).
eq_loop_1: V_V1-V_L1(t)-V_R2(t)=0;
eq_loop_2: V_L1(t)-V_C1(t)=0;
eq_loop_3: V_C1(t)-V_R1(t)=0;
eq_net_a: -I_V1(t)-I_L1(t)-I_C1(t)-I_R1(t)=0;
eq_net_b: -I_R2(t)+I_L1(t)+I_C1(t)+I_R1(t)=0;
eq_R1: V_R1(t)=R_R1*I_R1(t);
eq_R2: V_R2(t)=R_R2*I_R2(t);
eq_C1: I_C1(t)=C_C1*'diff(V_C1(t), t);
eq_L1: V_L1(t)=L_L1*'diff(I_L1(t), t);
assume(R_R1>0);
assume(R_R2>0);
assume(C_C1>0);
assume(L_L1>0);
L_loop_1: laplace(eq_loop_1, t, s);
L_loop_2: laplace(eq_loop_2, t, s);
L_loop_3: laplace(eq_loop_3, t, s);
L_net_a: laplace(eq_net_a, t, s);
L_net_b: laplace(eq_net_b, t, s);
L_R1: laplace(eq_R1, t, s);
L_R2: laplace(eq_R2, t, s);
L_C1: laplace(eq_C1, t, s);
L_L1: laplace(eq_L1, t, s);
eq_L_set: [L_loop_1, L_loop_2, L_loop_3, L_net_a, L_net_b, L_R1, L_R2, L_C1, L_L1];
x_L_set: [laplace(I_V1(t), t, s), laplace(V_R1(t), t, s), laplace(I_R1(t), t, s), laplace(V_C1(t), t, s), laplace(I_C1(t), t, s), laplace(V_L1(t), t, s), laplace(I_L1(t), t, s), laplace(V_R2(t), t, s), laplace(I_R2(t), t, s)];
L_solved: solve(eq_L_set, x_L_set);
solved_I_V1: ilt(L_solved[1][1],s,t);
solved_V_R1: ilt(L_solved[1][2],s,t);
solved_I_R1: ilt(L_solved[1][3],s,t);
solved_V_C1: ilt(L_solved[1][4],s,t);
solved_I_C1: ilt(L_solved[1][5],s,t);
solved_V_L1: ilt(L_solved[1][6],s,t);
solved_I_L1: ilt(L_solved[1][7],s,t);
solved_V_R2: ilt(L_solved[1][8],s,t);
solved_I_R2: ilt(L_solved[1][9],s,t);During processing the ilt()s Maxima asks 9 times:Is 4C_C1R_R1^2R_R2^2-L_L1R_R2^2-2L_L1R_R1R_R2-L_L1R_R1^2 positive, negative or zero?
Enter positive (and press CTRL-Enter) for each question.
L_loop_1 to L_L1 are equations in S-domain (laplace(f(t),t,s) converts function f(t) from time (t) domain to S-domain).
eq_L_set is the collection of all equations in S-domain.
x_L_set is the collection of all unknowns in S-domain. Note that as we don’t know (unless we peek the midway results) the S-domain representation of each equation alone, we are using the Laplace transformed equations (whatever they might be) as unknowns.
L_solved is a result got by solving all equations in S-domain. Note that S-domain equations are polynomial, so Maxima’s basic solve() procedure attacks them with ease.
Voilà, we got some results:
%e stands for e (2.718…)
C_C1 is printed as CC1.
The formula is long, but readable. Two special points to note:
IL1(0) and VC1(0). These are the initial values for the current through L1 and voltage across C1 when time was 0 seconds.
Other way to state this: If you know the voltage across capacitor C1 at time 0 seconds, and current through inductor L1 at time 0 seconds, you can accurately predict the behavior of current through V1 (IV1(t)) at any time (either future or past).
And the same goes with 8 other results.
You can also see the formula of negative/zero/positive question as part of argument to sin(). The argument of sin consists of fixed (non time dependent) values, and t (time dependent by definition).
The fixed part works as a frequency, which is depending on component values only. With certain values of components (typically big inductor compared to capacitor, I haven’t checked for this model) the argument to square root goes negative, meaning the circuit doesn’t show any oscillation effect with those component values.
If you enter negative to the question, sin() is changed to it’s hyperbolic version, sinh(), and the argument of square root is negated to positive.
Let’s cleanup the Maxima code a bit more, by adding lineassume(4C_C1R_R1^2R_R2^2-L_L1R_R2^2-2L_L1R_R1R_R2-L_L1R_R1^2 > 0);
after other assume()s:
eq_loop_1: V_V1-V_L1(t)-V_R2(t)=0;
eq_loop_2: V_L1(t)-V_C1(t)=0;
eq_loop_3: V_C1(t)-V_R1(t)=0;
eq_net_a: -I_V1(t)-I_L1(t)-I_C1(t)-I_R1(t)=0;
eq_net_b: -I_R2(t)+I_L1(t)+I_C1(t)+I_R1(t)=0;
eq_R1: V_R1(t)=R_R1*I_R1(t);
eq_R2: V_R2(t)=R_R2*I_R2(t);
eq_C1: I_C1(t)=C_C1*'diff(V_C1(t), t);
eq_L1: V_L1(t)=L_L1*'diff(I_L1(t), t);
assume(R_R1>0);
assume(R_R2>0);
assume(C_C1>0);
assume(L_L1>0);
assume(4*C_C1*R_R1^2*R_R2^2-L_L1*R_R2^2-2*L_L1*R_R1*R_R2-L_L1*R_R1^2 > 0);
L_loop_1: laplace(eq_loop_1, t, s);
L_loop_2: laplace(eq_loop_2, t, s);
L_loop_3: laplace(eq_loop_3, t, s);
L_net_a: laplace(eq_net_a, t, s);
L_net_b: laplace(eq_net_b, t, s);
L_R1: laplace(eq_R1, t, s);
L_R2: laplace(eq_R2, t, s);
L_C1: laplace(eq_C1, t, s);
L_L1: laplace(eq_L1, t, s);
eq_L_set: [L_loop_1, L_loop_2, L_loop_3, L_net_a, L_net_b, L_R1, L_R2, L_C1, L_L1];
x_L_set: [laplace(I_V1(t), t, s), laplace(V_R1(t), t, s), laplace(I_R1(t), t, s), laplace(V_C1(t), t, s),
laplace(I_C1(t), t, s), laplace(V_L1(t), t, s), laplace(I_L1(t), t, s), laplace(V_R2(t), t, s), laplace(I_R2(t), t, s)];
L_solved: solve(eq_L_set, x_L_set);
solved_I_V1: ilt(L_solved[1][1],s,t);
solved_V_R1: ilt(L_solved[1][2],s,t);
solved_I_R1: ilt(L_solved[1][3],s,t);
solved_V_C1: ilt(L_solved[1][4],s,t);
solved_I_C1: ilt(L_solved[1][5],s,t);
solved_V_L1: ilt(L_solved[1][6],s,t);
solved_I_L1: ilt(L_solved[1][7],s,t);
solved_V_R2: ilt(L_solved[1][8],s,t);
solved_I_R2: ilt(L_solved[1][9],s,t);The effect of added assume() is that we don’t need to tell whether the long formula is negative, zero, or positive. assume() does it for us in advance.
To get nicer looking results, you can also enter the code to online Maxima.
What these results actually mean?
These 9 formulas tell us how the 9 unknown (time dependent) variables vary over time (t seconds), based on:
– Inductor current at time zero seconds.
– Capacitor voltage at time zero seconds.
– Voltage of V1 (constant from zero seconds to our time of interest, t)
– Resistance of R1 (constant from zero seconds to our time of interest, t)
– Resistance of R2 (constant from zero seconds to our time of interest, t).
What if V1, R1 and R2 are changed from time to time?
In the LTSpice simulation above, the voltage of V1 varied the following way:
– At time 0, the voltage is 0 volts.
– At time 1 ms, the voltage is abruptly changed to 10 volts (in the simulation the transition period was 1 ns).
– At time 101 ms, the voltage is abruptly changed to 0 volts.
And our initial conditions were
– Capacitor’s voltage at time 0 is 0 volts.
– Inductor’s current at time 0 is 0 amperes.
To proceed, we do our calculation in periods:
First period. From 0 seconds to 1 ms:
The time dependent values can be calculated using the formulas, using initial values
IL1(0)=0 and VC1(0)=0.
The formulas are used to calculate IL1 and VC1 at the end of the period (at 1 ms).
Second period. From 1 ms to 101 ms:
The initial values used are the ones calculated at the end of previous period.
The time entered to the formulas is t-1 ms (1 ms is the star of this period). This way the formulas treat the start of this period as 0 seconds.
The formulas are used to calculate IL1 and VC1 at the end of the period (at 101 ms).
Third period. From 101 ms to positive infinity:
The initial values used are the ones calculated at the end of previous period.
The time entered to the formulas is t-101 ms (101 ms is the star of this period). This way the formulas treat the start of this period as 0 seconds.
The same change of other fixed values (RR1, RR2, CC1, LL1) can be done by creating a period boundary each time any of the fixed values changes.
VHDL of the same. I split this to two files.
circuit_math_pkg.vhdl:
The long formulas were copy-pasted from wxMaxima’s output, and search-replaced to fit the need.
library ieee;
use ieee.std_logic_1164.all;
use ieee.numeric_std.all;
use ieee.math_real.all;
package circuit_math_pkg is
type fixed_values_t is record
-- fixed values:
V_V1: real;
R_R1: real;
R_R2: real;
C_C1: real;
L_L1: real;
end record;
type circuit_state_t is record
timestamp: time;
fixed: fixed_values_t;
-- initial values:
init_V_C1: real;
init_I_L1: real;
-- time dependent values:
I_V1: real;
V_R1: real;
I_R1: real;
V_R2: real;
I_R2: real;
V_C1: real;
I_C1: real;
V_L1: real;
I_L1: real;
end record;
function circuit_evaluated(state: circuit_state_t; projected_time: time) return circuit_state_t;
end;
package body circuit_math_pkg is
function circuit_evaluated(state: circuit_state_t; projected_time: time) return circuit_state_t is
variable rval: circuit_state_t; -- Return value
variable dt_t: time; -- Time from the start of the period.
variable dt_r: real; -- Time from the start of the period (as real number).
variable t: real; -- used in Maxima generated formulas.
begin
dt_t := projected_time - state.timestamp;
dt_r := real(dt_t/1 ps) * 1.0e-12;
t := dt_r;
rval := state;
rval.timestamp := projected_time;
-- Using Notepad++ Replace with Regular expression to change from Maxima format to VHDL format:
-- Find what: [A-Z]_[A-Z0-9]*
-- Replace with: state.$0
--
-- plus a lot of manual editing...
rval.I_V1 :=
(exp(-((state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*t)/(2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2))*(((2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*(state.fixed.L_L1*state.fixed.R_R1*state.fixed.V_V1-state.init_I_L1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2)-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*(state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1))*sin((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2)))/(sqrt(state.fixed.L_L1)*sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0))+state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*cos((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2))))/(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2)-state.fixed.V_V1/state.fixed.R_R2;
rval.V_R1 :=
(exp(-((state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*t)/(2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2))*(((2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*(state.fixed.L_L1*state.fixed.R_R1*state.fixed.V_V1-state.init_I_L1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2)-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*(state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1))*sin((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2)))/(sqrt(state.fixed.L_L1)*sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0))+state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*cos((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2))))/(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2);
rval.I_R1 :=
(exp(-((state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*t)/(2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2))*(((2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*(state.fixed.L_L1*state.fixed.V_V1-state.init_I_L1*state.fixed.L_L1*state.fixed.R_R2)-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R2*(state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1))*sin((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2)))/(sqrt(state.fixed.L_L1)*sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0))+state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R2*cos((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2))))/(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2);
rval.V_C1 :=
(exp(-((state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*t)/(2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2))*(((2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*(state.fixed.L_L1*state.fixed.R_R1*state.fixed.V_V1-state.init_I_L1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2)-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*(state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1))*sin((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2)))/(sqrt(state.fixed.L_L1)*sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0))+state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*cos((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2))))/(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2);
rval.I_C1 :=
(exp(-((state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*t)/(2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2))*(((-(state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.V_V1+(-state.init_I_L1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1)*state.fixed.R_R2-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1)-2.0*state.init_V_C1*state.fixed.C_C1**2.0*state.fixed.L_L1*state.fixed.R_R1**2.0*state.fixed.R_R2**2.0)*sin((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2)))/(sqrt(state.fixed.L_L1)*sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0))+(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.V_V1+(-state.init_I_L1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1)*state.fixed.R_R2-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1)*cos((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2))))/(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2);
rval.V_L1 :=
(exp(-((state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*t)/(2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2))*(((2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*(state.fixed.L_L1*state.fixed.R_R1*state.fixed.V_V1-state.init_I_L1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2)-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*(state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1))*sin((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2)))/(sqrt(state.fixed.L_L1)*sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0))+state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*cos((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2))))/(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2);
rval.I_L1 :=
(exp(-((state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*t)/(2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2))*(((-2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*((state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*state.fixed.V_V1+(-state.init_V_C1*state.fixed.C_C1*state.fixed.R_R1-state.init_I_L1*state.fixed.L_L1)*state.fixed.R_R2**2.0-state.init_I_L1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2)-(state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*(state.init_I_L1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.V_V1))*sin((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2)))/(sqrt(state.fixed.L_L1)*sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0))+(state.init_I_L1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.V_V1)*cos((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2))))/(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2)+state.fixed.V_V1/state.fixed.R_R2;
rval.V_R2 :=
(exp(-((state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*t)/(2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2))*(((2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*(state.init_I_L1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1*state.fixed.V_V1)+state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*(state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1))*sin((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2)))/(sqrt(state.fixed.L_L1)*sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0))-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2*cos((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2))))/(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2)+state.fixed.V_V1;
rval.I_R2 :=
(exp(-((state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)*t)/(2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2))*(((state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*(state.fixed.L_L1*state.fixed.R_R2+state.fixed.L_L1*state.fixed.R_R1)-2.0*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*(state.fixed.L_L1*state.fixed.R_R1*state.fixed.V_V1-state.init_I_L1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2))*sin((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2)))/(sqrt(state.fixed.L_L1)*sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0))-state.init_V_C1*state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*cos((sqrt((4.0*state.fixed.C_C1*state.fixed.R_R1**2.0-state.fixed.L_L1)*state.fixed.R_R2**2.0-2.0*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2-state.fixed.L_L1*state.fixed.R_R1**2.0)*t)/(2.0*state.fixed.C_C1*sqrt(state.fixed.L_L1)*state.fixed.R_R1*state.fixed.R_R2))))/(state.fixed.C_C1*state.fixed.L_L1*state.fixed.R_R1*state.fixed.R_R2)+state.fixed.V_V1/state.fixed.R_R2;
-- Save the new V_C1 and I_L1 as new initial values.
rval.init_V_C1 := rval.V_C1;
rval.init_I_L1 := rval.I_L1;
return rval;
end;
end;
circuit_simulator.vhdl:
library ieee;
use ieee.std_logic_1164.all;
use ieee.numeric_std.all;
use ieee.math_real.all;
use work.circuit_math_pkg.all;
entity circuit_simulator is
end;
architecture tb of circuit_simulator is
function fixed_init_val return fixed_values_t is
variable r: fixed_values_t;
begin
r := (others => 0.0);
r.V_V1 := 0.0; -- Here you can enter other value.
r.R_R1 := 270.0; -- Here you can enter other value.
r.R_R2 := 100.0; -- Here you can enter other value.
r.C_C1 := 33.0e-6; -- Here you can enter other value.
r.L_L1 := 10.0e-3; -- Here you can enter other value.
return r;
end;
function circuit_init_val return circuit_state_t is
variable r: circuit_state_t;
begin
r := (timestamp => 0 ns, fixed => fixed_init_val, others => 0.0);
r.init_V_C1 := 0.0; -- Here you can enter other initial value.
r.init_I_L1 := 0.0; -- Here you can enter other initial value.
return r;
end;
signal circuit_state: circuit_state_t := circuit_init_val;
signal fixed_state: fixed_values_t := fixed_init_val;
signal adc_clk: std_logic;
signal adc_analog_input_value: real := 0.0;
begin
circuit_evaluator_pr: process(fixed_state)
variable new_circuit_state: circuit_state_t;
begin
new_circuit_state := circuit_evaluated(circuit_state, now);
new_circuit_state.fixed := fixed_state;
circuit_state <= new_circuit_state;
end process;
input_waveform_pr: process
begin
fixed_state <= fixed_init_val;
wait for 1 ms;
fixed_state.V_V1 <= 10.0;
wait for 100 ms;
fixed_state.V_V1 <= 0.0;
wait for 99 ms;
fixed_state.V_V1 <= 0.0;
wait;
end process;
adc_clk_pr: process
begin
adc_clk <= '0';
wait for 15 us;
adc_clk <= '1';
wait for 15 us;
end process;
adc_evaluate_pr: process(adc_clk)
variable eval_state: circuit_state_t;
begin
if rising_edge(adc_clk) then
eval_state := circuit_evaluated(circuit_state, now);
adc_analog_input_value <= eval_state.V_R2;
end if;
end process;
end;
Concetrate first on to circuit_simulator.vhdl file
The process:
circuit_evaluator_pr: process(fixed_state)
variable new_circuit_state: circuit_state_t;
begin
new_circuit_state := circuit_evaluated(circuit_state, now);
new_circuit_state.fixed := fixed_state;
circuit_state <= new_circuit_state;
end process;does most of the trick.
Every time signal fixed_state changes (all fixed values, like input voltage are stored there), the new circuit_state is evaluated (calculated by the long formulas), and stored as a state from this point onwards.
The imagined ADC converter needs the output voltage (net_b, which is equal to V_R2(t)), but only when adc_clk has rising edge. The process:
adc_evaluate_pr: process(adc_clk)
variable eval_state: circuit_state_t;
begin
if rising_edge(adc_clk) then
eval_state := circuit_evaluated(circuit_state, now);
adc_analog_input_value <= eval_state.V_R2;
end if;
end process;Calculates the adc_analog_input_value only when needed (rising edge of adc_clk).
You can use any decent VHDL simulator to simulate this model. GHDL (assuming Linux here) is used by entering:
ghdl -a --std=08 circuit_math_pkg.vhdl
ghdl -a --std=08 circuit_simulator.vhdl
ghdl -e --std=08 circuit_simulator
ghdl -r --std=08 circuit_simulator --stop-time=200ms --wave=simulation.ghw
Then the result can be visualized bygtkwave simulation.ghw &
And after setting the gtkwave’s output signals, the output looks like:
This is very close to the waveform given out by LTSpice.
If the ADC clock period is increased to much bigger steps, the analog value is evaluated only on those times, without loosing any accuracy on the actual value of the analog signal.
Please note that this is the about the most complex circuit which have closed form analytical results. If you try to add another capacitor or inductor to the circuit, the inverse Laplace transform is probably not available (it does exist, but not in the language of math we know of).
2 Comments
Philippe · May 10, 2022 at 2:00 pm
Interesting article.
Only V_C1, I_L1 (to update init values) and V_R2 need to be calculated for adc_analog_input_value
Regards,
Philippe
Philippe · May 11, 2022 at 10:16 am
Hello again,
I ‘m surprised that you solve the differential equations (using Maxima) and after, for simulation, you use delta time ( projected_time – state.timestamp ) and change init variables at each step.
In my opinion, this method is complex.
Instead I successed to use original equations to solve the circuit:
equations are (extra variables need to be defined before …):
dt_t := projected_time – state.timestamp;
dt := real(dt_t/1 ps) * 1.0e-12;
rval := state;
rval.timestamp := projected_time;
il := state.V_C1/state.fixed.L_L1 * dt + state.I_L1; — compute new current in L1
ir := state.V_C1/state.fixed.R_R1;
i := (state.fixed.V_V1 – state.V_C1)/state.fixed.R_R2;
ic := i – il – ir;
rval.V_C1 := ic/state.fixed.C_C1 * dt + state.V_C1; — compute new voltage on C1
rval.I_L1 := il;
rval.V_R2 := i * state.fixed.R_R2;
for circuit_simulator.vhdl file, the only modification is to add adc_clk in the sensitivity list of the process “circuit_evaluator_pr”. So each “delta time” lasts 15 ns /2.