Learn › Engineering

Multi-domain simulation: how the simulator computes

The multi-domain simulator puts circuits, motors, heaters, tanks, sensors and control blocks in the same schematic and solves them all at once in time. It has no detailed PDF solution because there is no «hand solution»: it is a numerical integrator. This topic explains what it does inside, so you can read its results, trust them and recognize its limits.

Go to the simulator

One language for several domains

Every physical system it simulates is described with two kinds of variables: an effort or difference variable (voltage, angular speed, temperature, level or pressure) and a flow variable (current, torque, heat flow, flow rate). With them, the elements of each domain look like those of an electric circuit, and their equations have the same form:

Domain	Effort · flow	Resistance	Stores (capacity)	Storage equation
Electrical	V · i	$R$ : $V = R i$	capacitor $C$	$i = C \frac{d V}{d t}$
Mechanical (rotation)	ω · T	friction $b$ : $T = bω$	inertia $J$	$T = J \frac{d ω}{d t}$
Thermal	T · q	$R_{t h}$ : $Δ T = R_{t h} q$	heat capacity $C_{t h}$	$q = C_{t h} \frac{d T}{d t}$
Hydraulic	h · Q	outlet $Q = k h$ (nonlinear)	tank area $A$	$Q = A \frac{d h}{d t}$

That is why an RC circuit, a heater that warms up and a tank that fills respond alike, with an exponential of time constant $τ = R \cdot C$ (or $R_{t h} C_{t h}$ , or $2 A h / k$ ). The components that convert energy join two domains: the DC motor produces a torque $T_{m} = K i$ and at the same time generates a back-EMF $K ω$ :

L \frac{d i}{d t} + R i + K ω = v J \frac{d ω}{d t} + b ω + T_{L} = K i

DC motor: the current produces torque, the speed produces back-EMF, and both equations are solved together.

How it builds the system of equations

The unknowns are the voltages of the nodes of the schematic and the state variables of the components that store energy (a capacitor's voltage, an inductor's current, a motor's current, speed and angle, a heater's temperature, a tank's level…). They all go into a vector $z$ . Each component contributes two things:

The currents that enter or leave through its terminals. At each node their sum must be zero (Kirchhoff's current law): that is nodal analysis.
Its internal equations, written as «residual = 0»: for a capacitor $i - C \overset{v}{˙} = 0$ ; for an inductor $L \overset{}{˙} + R_{s} i - v = 0$ ; for a heater $C_{t h} \dot{T} - [P - (T - T_{a}) / R_{t h}] = 0$ .

Putting everything together gives a system of algebraic and differential equations with the general form

F (z, \overset{z}{˙}, t) = 0

This format accepts any component without special treatment: a diode (Shockley equation), a transistor, an op-amp with saturation or a sensor with its own delay are just more equations in $F$ . Voltage sources have a small internal resistance (0.01 Ω by default), like real ones, and a tiny conductance from every node to ground ( $1 0^{- 12}$ S) is added so that an isolated node does not leave the system without a solution.

Time integration: BDF2 and Newton

Time advances in steps $h$ . In each step the derivatives $\overset{z}{˙}$ are replaced by a formula with the values of the previous steps, and what remains is an algebraic system in $z_{n + 1}$ . The simulator uses the second-order backward differentiation formula (BDF2), which takes the slope at the new point of the parabola through the last three:

\overset{z}{˙}_{n + 1} \approx \frac{3 z _{n + 1} - 4 z _{n} + z _{n - 1}}{2 h}

BDF2: the derivative at tₙ₊₁ is the slope of the parabola through the last three points.

BDF2 is implicit: the unknown appears on both sides, so $F (z_{n + 1}) = 0$ must be solved in every step. It is done with Newton's method: the residual $F$ is evaluated, its Jacobian $J = \partial F / \partial z$ is computed (by differentiating each component numerically), $J Δ z = - F$ is solved and $z$ is corrected until the correction is negligible. Implicit methods are stable even when the system has very fast dynamics (a circuit with microsecond constants inside a thermal process of minutes): the step does not have to be smaller than the fastest constant.

One step of the simulator. The first step (no history) uses backward Euler; if Newton does not converge, the step is split into two halves, up to 8 times.

The step size. By default it is the simulation length divided by 2000, and it is reduced if there are fast sources, sampling or delays. It is also adjusted to the circuit's oscillations: BDF2 damps an oscillation of frequency $ω$ by a factor of about $(ωh)^{4} /4$ per step and delays it by about $(2/9) (ωh)^{3}$ radians, so the simulator chooses $h$ so that, over the whole run, that numerical damping stays under 1% and the phase lag under 0.01 rad. Digital components (clocks, counters, the 555, comparators with hysteresis) change through events at the end of each accepted step.

Linearization: poles, stability and G(s)

At the end of the simulation (or at the chosen instant) the simulator linearizes the system around the state it ended in: it differentiates the equations numerically with respect to the states and the input, and obtains a linear model

\overset{x}{˙} = A x + B u y = C x + D u G (s) = C (s I - A)^{- 1} B + D

The eigenvalues of A are the poles of the whole system; if they all have negative real part, that operating point is stable. The same matrix gives the own poles of each subsystem (the network of capacitors and inductors, a motor, a PID…), the transfer function between the input and output you choose, its step response, the Bode diagram and the root locus with respect to any parameter (it is relinearized for each value).

A linear model is only valid near the point where it was linearized. A diode, a saturated motor, an ON/OFF control or a tank's

k h

outlet change behavior far from that point: that is why the simulator compares, when it can, the response of the linear model with the nonlinear simulation, and warns when the poles say «stable» but the simulation oscillates (a limit cycle).

How to read and check the results

Units and order of magnitude: before looking at decimals, check with a quick calculation (V/R, $K V / (R b + K^{2})$ , $P R_{t h}$ ) that the result makes sense.
Duration: simulate at least 5 time constants of the slowest dynamics to see the final value; the poles card gives you those constants ( $τ = 1/∣ σ ∣$ ).
Step: if you doubt a fine detail (a peak, a ripple), repeat with a smaller step; if the result does not change, it comes from the model and not from the numerical method.
Warnings: the simulator warns about components beyond their limits (a resistor that would burn, a stalled motor, a tank that overflows) and about badly posed schematics (shorted sources or sources in parallel with different values).

Common mistakes

Forgetting the ground: without a reference node the voltages are not defined.
Reading a «final» value before the system settles: look at the whole response or lengthen the simulation.
Using the transfer function far from the operating point of a nonlinear system.
Expecting the accuracy of an exact calculation at the instant of a switching event: right after a jump in the source, BDF2 mixes the pre-jump history for a couple of steps.
Confusing the pole of an integrator (at 0, «marginal») with an instability: what matters is the closed loop.

Worked examples

Example 1Charging and discharging a capacitor

A pulse source from 0 to 5 V at 2 Hz feeds a 10 kΩ resistor in series with a 10 µF capacitor. Find the time constant, the pole and the maximum and minimum capacitor voltages in the periodic regime.

Solution

$τ = RC = 1 0^{4} \cdot 1 0^{- 5} = 0.1$ s; the pole is at $- 1/ τ = - 10$ s⁻¹ and $G (s) = V_{C} / V = \frac{10}{s + 10}$ .
Each half period lasts 0.25 s = 2.5τ. In the periodic regime, the charge goes from $V_{min}$ to $V_{ma x}$ and the discharge comes back: $V_{ma x} = 5 + (V_{min} - 5) e^{- 2.5}$ and $V_{min} = V_{ma x} e^{- 2.5}$ .
Solving: $V_{ma x} = \frac{5 ( 1 - e ^{- 2.5} )}{1 - e ^{- 5}} = 4.621$ V and $V_{min} = 4.621 e^{- 2.5} = 0.379$ V.
Step response of the linear model: 10–90% rise in $τ ln 9 = 220$ ms and 2% settling in $τ ln 50 = 391$ ms.

τ = 0.1

s, pole

- 10

; in the periodic regime

V_{C}

swings between 0.379 and 4.621 V.

Solve it in Simuveliq

In Engineering → Multi-domain simulation, in Examples, «Load example» → Capacitor charge and discharge (RC). The Oscilloscope shows the source wave and the capacitor's (peaks of 4.62 V and minima of 0.38 V with the 1 ms step).
Stability: whole system and subsystems gives the pole −10 with τ = 100 ms and «Follows the periodic source»; Control analysis gives G(s) = 10/(s + 10) and, for the linear model, rise 220 ms, settling 392 ms and final value 1.

Open in the simulator →

Example 2Underdamped RLC circuit

A 1 V pulse feeds R = 10 Ω, L = 10 mH and C = 100 µF in series (the source has 0.01 Ω inside). Find the poles, the damping, the overshoot of the capacitor voltage and the frequency of the oscillation.

Solution

Loop: $L \overset{q}{¨} + R \overset{q}{˙} + q / C = v$ , so $\frac{V _{C}}{V} = \frac{1/ L C}{s ^{2} + ( R / L ) s + 1/ L C}$ .
$ω_{n} = 1/ L C = 1000$ rad/s. With $R = 10.01$ Ω (10 + 0.01 from the source): $σ = R /2 L = 500.5$ s⁻¹, $ζ = σ / ω_{n} = 0.5005$ and $ω_{d} = ω_{n}^{2} - σ^{2} = 865.7$ rad/s. Poles: $- 500.5 \pm 865.7 j$ .
$M_{p} = e^{- π ζ / 1 - ζ^{2}} = 16.3%$ : the capacitor reaches 1.163 V; $t_{p} = π / ω_{d} = 3.63$ ms. The oscillation has $f_{d} = ω_{d} /2 π = 137.8$ Hz.

Poles

- 500.5 \pm 865.7 j

(

ζ = 0.5

ω_{n} = 1000

rad/s);

V_{C}

peaks at 1.163 V 3.6 ms after the edge.

Solve it in Simuveliq

«Load example» → Underdamped RLC circuit (complex poles). The dominant pole is −500.5 ± 866j with ζ = 0.5 and ωn = 1000 rad/s; on the oscilloscope the capacitor voltage reaches 1.163 V after each rising edge.
Change the source's internal resistance to 0 or R to 20 Ω and watch the poles move (with R = 20 Ω, ζ ≈ 1: there is no overshoot anymore).

Open in the simulator →

Example 3DC motor at 12 V

A motor with

R = 2

Ω,

L = 1

mH,

K = 0.05

V·s/rad,

J = 1 0^{- 4}

kg·m² and

b = 1 0^{- 4}

N·m·s is connected to a 12 V battery (0.01 Ω inside). Find the final speed, the final and starting current, and the poles.

Solution

In steady state the derivatives are zero: $R i + K ω = V$ and $b ω = K i$ . With $R = 2.01$ Ω: $ω = \frac{K V}{R b + K ^{2}} = \frac{0.05 \cdot 12}{2.01 \cdot 1 0 ^{- 4} + 0.0025} = 222.1$ rad/s $= 2121$ rpm, and $i = bω / K = 0.444$ A.
At start $ω = 0$ and there is no back-EMF: the current rises toward $V / R = 5.97$ A, held back by the inductance (the simulated peak is 5.83 A). That is why motors are protected at start-up.
Poles: $(L s + R) (J s + b) + K^{2} = 0$ → $s^{2} + 2011 s + 27010 = 0$ → $s = - 13.52$ (mechanical, τ = 74 ms) and $s = - 1997$ (electrical, 0.5 ms). The slow one dominates: the speed reaches 63% in about 74 ms.
Voltage-to-speed gain: $\frac{2121}{12} = 176.8$ rpm/V.

ω_{\infty} = 2121

rpm,

i_{\infty} = 0.444

A, start-up ≈ 5.8–6 A; poles −13.52 and −1997.

Solve it in Simuveliq

The button opens a motor M1 with a 12 V battery. The oscilloscope shows the speed, which reaches about 2120 rpm at 0.5 s, and the current, with its 5.83 A start-up peak dropping to 0.45 A.
Stability: poles −13.52 and −1997 in the «DC motor» subsystem (the angle θ only accumulates and does not count). Control analysis from V1 to the speed: G(s) = 4 775 000/(s² + 2011s + 27 010), DC gain 176.8 and final value of the linear model 176.8.

Open in the simulator →

Example 4Heater: a first-order thermal system

A 12 Ω heater with

C_{t h} = 20

J/K and

R_{t h} = 5

K/W, in a 25 °C ambient, is connected to 12 V. Find the power, the final temperature, the time constant and the temperature at 100 s.

Solution

Voltage across the heater: $12 \cdot 12/12.01 = 11.990$ V, so $P = V^{2} / R = 11.98$ W.
Energy balance: $C_{t h} \dot{T} = P - \frac{T - T _{a}}{R _{t h}}$ . It is a first-order system with $τ = R_{t h} C_{t h} = 100$ s and final $T_{\infty} = T_{a} + P R_{t h} = 25 + 59.9 = 84.9$ °C.
$T (100 s) = 25 + 59.9 (1 - e^{- 1}) = 62.9$ °C.
The power depends on $V^{2}$ (nonlinear). Linearizing at 12 V: $\frac{d T _{\infty}}{d V} = R_{t h} \frac{2 V _{h}}{R} \cdot \frac{12}{12.01} = 9.98$ °C/V: one more volt would raise it about 10 °C.

P = 11.98

τ = 100

T_{\infty} = 84.9

°C,

T (100 s) = 62.9

°C; linear gain 9.98 °C/V.

Solve it in Simuveliq

The button opens heater H1 with its battery; the run lasts 600 s (6τ) and ends at 84.75 °C, about to reach 84.9. The oscilloscope gives the power, 11.98 W.
Stability: pole −0.01 with τ = 100 s, «Thermal time constant τ = C·Rth». Control analysis from V1 to the temperature: G(s) = 0.09983/(s + 0.01), DC gain 9.983.

Open in the simulator →

Example 5Tank: linearizing a nonlinear outlet

A tank of area

A = 0.01

m² drains through an outlet

Q = k h

with

k = 4 \times 1 0^{- 4}

m^2.5/s, and a PID with a pressure sensor holds the level at 0.5 m. What time constant does the tank alone have around that level?

Solution

Volume balance: $A \dot{h} = Q_{in} - k h$ . It is nonlinear because of the root.
Near $h_{0}$ the curve is replaced by its tangent: $k h \approx k h_{0} + \frac{k}{2 h _{0}} (h - h_{0})$ . With $δ h = h - h_{0}$ : $A δ \dot{h} = δ Q_{in} - \frac{k}{2 h _{0}} δ h$ .
Pole: $s = - \frac{k}{2 A h _{0}}$ . At the level the simulation reaches, $h_{0} = 0.499$ m: $s = - \frac{4 \times 1 0 ^{- 4}}{2 \cdot 0.01 \cdot 0.7064} = - 0.0283$ s⁻¹, $τ = 35.3$ s. At another level the pole changes: the tank is slower the fuller it is.
That is why the linear model and the simulation do not fully agree: the reference step starts with the tank far from 0.5 m (overshoot 15.6% in the linear model and 5.31% in the simulation).

Around 0.5 m the tank has a pole at −0.0283 s⁻¹ (

τ = 35.3

s), which depends on the level.

Solve it in Simuveliq

«Load example» → Tank level control with a pressure sensor. In the subsystem table, «TK1 Tank — Level (linearized k√h) and valve» has the poles −0.02831 (level) and −2 (the valve, with a 0.5 s response).
In Control analysis, the table compares the nonlinear simulation with the linear model: overshoot 5.31% versus 15.6%.

Open in the simulator →

Multi-domain simulation: how the simulator computes

One language for several domains

How it builds the system of equations

Time integration: BDF2 and Newton

Linearization: poles, stability and G(s)

How to read and check the results

Common mistakes

Worked examples

Example 1Charging and discharging a capacitor

Solution

Solve it in Simuveliq

Example 2Underdamped RLC circuit

Solution

Solve it in Simuveliq

Example 3DC motor at 12 V

Solution

Solve it in Simuveliq

Example 4Heater: a first-order thermal system

Solution

Solve it in Simuveliq

Example 5Tank: linearizing a nonlinear outlet

Solution

Solve it in Simuveliq

Related topics