Skip to content

Power systems

Power systems: generation, transmission, and recovery of electric energy.

faraday_disk(I=0.0, U=0.0, omega=0.0, *, B=1.5, M=4.8, N=250.0, R1=2.3, R2=109.0, L=0.2, C=0.75)

Return a battery charger driven by a Faraday disk generator whose own current strengthens its field.

A motor turns the conducting disk of radii r1 and r2 with a constant torque M in the field B of a magnet plus the field μ N I / l of a coil that the generated current flows through. With a = B (r2^2 - r1^2) / 2 and b = μ N (r2^2 - r1^2) / (2 l), the current, the polarization voltage of the battery, and the angular speed follow L I' = -R1 I - U + (a + b I) ω - E(I), C U' = I - U / R2, and J ω' = M - alpha sign(ω) - β r2 ω - (a + b I) I, where the battery voltage is E(I) = a0 + a1 I + a3 I^3 and the last term is the braking torque of the Lorentz force.

The generator charges the battery only when its voltage a ω exceeds the open-circuit voltage a0 = 8 of the battery, at ω above 35 for the default field. Below that speed the steady current is a small negative current, the battery discharges through the disk and helps the motor, and the steady speed grows linearly with the torque, as (M - alpha) / (β r2) up to the small Lorentz torque: with the defaults the current settles at -0.040 and the speed at 15.7. The coil acts as a negative resistance b ω against R1 + a1, so above ω = (R1 + a1) / b, about 32 or a torque near 9.5 for the defaults, the steady state loses its stability and the current oscillates around zero, swinging between charging and discharging, while the mean braking torque b <I²> bends the speed below the linear trend. A stronger magnet charges at the default torque: B = 4 gives a steady current of 0.013.

Parameters:

Name Type Description Default
I float

The initial charging current.

 0.0
U float

The initial polarization voltage of the battery.

 0.0
omega float

The initial angular speed of the disk.

 0.0
B float

The field of the permanent magnet.

 1.5
M float

The torque of the drive motor.

 4.8
N float

The number of turns of the feedback coil.

 250.0
R1 float

The series resistance of the coil, the disk, and the battery.

 2.3
R2 float

The polarization resistance of the battery.

 109.0
L float

The inductance of the coil.

 0.2
C float

The polarization capacitance of the battery.

 0.75

Returns:

Name Type Description
system System

The system, with the current charger.I, the polarization voltage charger.U, and the angular speed charger.omega.

Examples:

>>> from skcomplex.systems.engineering.power_systems import faraday_disk
>>> faraday_disk(M=6.0).elements[0].M
6.0
Source code in src/skcomplex/systems/engineering/power_systems/_electromechanics.py 
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
@register_system(
    field='engineering/power_systems',
    tags=('energy', 'electromechanics', 'batteries', 'ode'),
    references=('Trzaska 2023',),
    defaults={'until': 50.0, 'dt': 0.1},
)
def faraday_disk(
    I: float = 0.0,
    U: float = 0.0,
    omega: float = 0.0,
    *,
    B: float = 1.5,
    M: float = 4.8,
    N: float = 250.0,
    R1: float = 2.3,
    R2: float = 109.0,
    L: float = 0.2,
    C: float = 0.75,
) -> System:
    """Return a battery charger driven by a Faraday disk generator whose own current strengthens its field.

    A motor turns the conducting disk of radii `r1` and `r2` with a constant torque `M` in the field `B` of a magnet
    plus the field `μ N I / l` of a coil that the generated current flows through. With `a = B (r2^2 - r1^2) / 2` and
    `b = μ N (r2^2 - r1^2) / (2 l)`, the current, the polarization voltage of the battery, and the angular speed follow
    `L I' = -R1 I - U + (a + b I) ω - E(I)`, `C U' = I - U / R2`, and
    `J ω' = M - alpha sign(ω) - β r2 ω - (a + b I) I`, where the battery voltage is `E(I) = a0 + a1 I + a3 I^3` and the
    last term is the braking torque of the Lorentz force.

    The generator charges the battery only when its voltage `a ω` exceeds the open-circuit voltage `a0 = 8` of the
    battery, at `ω` above 35 for the default field. Below that speed the steady current is a small negative current,
    the battery discharges through the disk and helps the motor, and the steady speed grows linearly with the torque,
    as `(M - alpha) / (β r2)` up to the small Lorentz torque: with the defaults the current settles at -0.040 and the
    speed at 15.7. The coil acts as a negative resistance `b ω` against `R1 + a1`, so above `ω = (R1 + a1) / b`, about
    32 or a torque near 9.5 for the defaults, the steady state loses its stability and the current oscillates around
    zero, swinging between charging and discharging, while the mean braking torque `b <I²>` bends the speed below the
    linear trend. A stronger magnet charges at the default torque: `B = 4` gives a steady current of 0.013.

    Args:
        I:
            The initial charging current.

        U:
            The initial polarization voltage of the battery.

        omega:
            The initial angular speed of the disk.

        B:
            The field of the permanent magnet.

        M:
            The torque of the drive motor.

        N:
            The number of turns of the feedback coil.

        R1:
            The series resistance of the coil, the disk, and the battery.

        R2:
            The polarization resistance of the battery.

        L:
            The inductance of the coil.

        C:
            The polarization capacitance of the battery.

    Returns:
        system:
            The system, with the current `charger.I`, the polarization voltage `charger.U`, and the angular speed
            `charger.omega`.

    Examples:
        >>> from skcomplex.systems.engineering.power_systems import faraday_disk
        >>> faraday_disk(M=6.0).elements[0].M
        6.0
    """
    charger = FaradayDiskElement('charger', I=I, U=U, omega=omega, B=B, M=M, N=N, R1=R1, R2=R2, L=L, C=C)
    return System(elements=[charger])

power_grid(n_nodes=8, *, K=2.75, power=1.0, gamma=0.1, added_line=None, doubled_line=None)

Return a ring power grid in the swing-equation system, where more transmission capacity can break synchrony.

The phase of node i relative to the grid frequency follows phi_i'' = P_i - gamma phi_i' + Σ_j K_ij sin(phi_j - phi_i). The first half of the ring are generators with P = power and the second half consumers with P = -power, joined by lines of capacity K. The default grid synchronises, but adding the line (0, 2) or doubling the capacity of the line (0, 1) makes the flows overload a line, and synchrony is lost, which is Braess's paradox.

Parameters:

Name Type Description Default
n_nodes int

The number of nodes of the ring.

 8
K float

The capacity of every line.

 2.75
power float

The power of every generator, and its negative the power of every consumer.

 1.0
gamma float

The damping.

 0.1
added_line tuple[int, int] | None

A pair of nodes to join by an extra line of capacity K, or None.

 None
doubled_line tuple[int, int] | None

A pair of neighbouring nodes whose line has capacity 2 K, or None.

 None

Returns:

Name Type Description
system System

The system, with the phases node_<i>.phi and frequency deviations node_<i>.phi_d1.

Examples:

>>> from skcomplex.systems.engineering.power_systems import power_grid
>>> len(power_grid(added_line=(0, 2)).interactions)
9
Source code in src/skcomplex/systems/engineering/power_systems/_grids.py 
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
@register_system(
    field='engineering/power_systems',
    tags=('power grid', 'synchronization', 'network', 'braess paradox', 'ode'),
    references=('Witthaut and Timme 2012', 'Datseris and Parlitz 2022'),
    defaults={'until': 300.0, 'dt': 1.0},
)
def power_grid(
    n_nodes: int = 8,
    *,
    K: float = 2.75,
    power: float = 1.0,
    gamma: float = 0.1,
    added_line: tuple[int, int] | None = None,
    doubled_line: tuple[int, int] | None = None,
) -> System:
    """Return a ring power grid in the swing-equation system, where more transmission capacity can break synchrony.

    The phase of node `i` relative to the grid frequency follows
    `phi_i'' = P_i - gamma phi_i' + Σ_j K_ij sin(phi_j - phi_i)`.
    The first half of the ring are generators with `P = power` and the second half consumers with `P = -power`, joined
    by lines of capacity `K`. The default grid synchronises, but adding the line `(0, 2)` or doubling the capacity of
    the line `(0, 1)` makes the flows overload a line, and synchrony is lost, which is Braess's paradox.

    Args:
        n_nodes:
            The number of nodes of the ring.

        K:
            The capacity of every line.

        power:
            The power of every generator, and its negative the power of every consumer.

        gamma:
            The damping.

        added_line:
            A pair of nodes to join by an extra line of capacity `K`, or `None`.

        doubled_line:
            A pair of neighbouring nodes whose line has capacity `2 K`, or `None`.

    Returns:
        system:
            The system, with the phases `node_<i>.phi` and frequency deviations `node_<i>.phi_d1`.

    Examples:
        >>> from skcomplex.systems.engineering.power_systems import power_grid
        >>> len(power_grid(added_line=(0, 2)).interactions)
        9
    """
    elements = [
        GridNodeElement(
            f'node_{index}',
            phi=0.0,
            phi_d1=0.0,
            P=power if index < n_nodes // 2 else -power,
            gamma=gamma,
        )
        for index in range(n_nodes)
    ]
    pairs = [(index, (index + 1) % n_nodes) for index in range(n_nodes)]
    doubled = None if doubled_line is None else frozenset(doubled_line)
    interactions = [
        TransmissionLineInteraction(
            f'line_{first}_{second}',
            ends=[f'node_{first}', f'node_{second}'],
            K=2 * K if frozenset((first, second)) == doubled else K,
        )
        for first, second in pairs
    ]
    if added_line is not None:
        first, second = added_line
        interactions.append(
            TransmissionLineInteraction(f'line_{first}_{second}', ends=[f'node_{first}', f'node_{second}'], K=K),
        )
    return System(elements=elements, interactions=interactions)

regenerative_braking(v=20.0, *, m=1000.0, A=1.2, L=0.01, R=0.05, C=680.0, alpha=50.0, beta=5.0)

Return a car braked by a Faraday disk generator on its wheel, which charges a supercapacitor.

The disk generates the voltage A v, with the structural constant A = |B| R^2 / (2 r_w) of the disk and the wheel, and drives the current through an inductor and a resistor into a supercapacitor, L I' = A v - R I - V_s and C V_s' = I. The current brakes the car with the force A I, together with dry and viscous friction, m v' = -A I - alpha sign(v) - β v, so the kinetic energy of the car ends up in the capacitor and in the losses.

Parameters:

Name Type Description Default
v float

The initial speed of the car.

 20.0
m float

The mass of the car.

 1000.0
A float

The structural constant of the generator, the voltage per unit speed.

 1.2
L float

The inductance of the charging circuit.

 0.01
R float

The resistance of the charging circuit.

 0.05
C float

The capacitance of the supercapacitor.

 680.0
alpha float

The dry friction force.

 50.0
beta float

The viscous friction coefficient.

 5.0

Returns:

Name Type Description
system System

The system, with the speed car.v, the current car.I, and the capacitor voltage car.Vs.

Examples:

>>> from skcomplex.systems.engineering.power_systems import regenerative_braking
>>> regenerative_braking(C=100.0).elements[0].C
100.0
Source code in src/skcomplex/systems/engineering/power_systems/_electromechanics.py 
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
@register_system(
    field='engineering/power_systems',
    tags=('energy', 'electromechanics', 'vehicles', 'ode'),
    references=('Trzaska 2023',),
    defaults={'until': 60.0, 'dt': 0.5},
)
def regenerative_braking(
    v: float = 20.0,
    *,
    m: float = 1000.0,
    A: float = 1.2,
    L: float = 0.01,
    R: float = 0.05,
    C: float = 680.0,
    alpha: float = 50.0,
    beta: float = 5.0,
) -> System:
    """Return a car braked by a Faraday disk generator on its wheel, which charges a supercapacitor.

    The disk generates the voltage `A v`, with the structural constant `A = |B| R^2 / (2 r_w)` of the disk and the
    wheel, and drives the current through an inductor and a resistor into a supercapacitor,
    `L I' = A v - R I - V_s` and `C V_s' = I`. The current brakes the car with the force `A I`, together with dry and
    viscous friction, `m v' = -A I - alpha sign(v) - β v`, so the kinetic energy of the car ends up in the capacitor
    and in the losses.

    Args:
        v:
            The initial speed of the car.

        m:
            The mass of the car.

        A:
            The structural constant of the generator, the voltage per unit speed.

        L:
            The inductance of the charging circuit.

        R:
            The resistance of the charging circuit.

        C:
            The capacitance of the supercapacitor.

        alpha:
            The dry friction force.

        beta:
            The viscous friction coefficient.

    Returns:
        system:
            The system, with the speed `car.v`, the current `car.I`, and the capacitor voltage `car.Vs`.

    Examples:
        >>> from skcomplex.systems.engineering.power_systems import regenerative_braking
        >>> regenerative_braking(C=100.0).elements[0].C
        100.0
    """
    car = RegenerativeCarElement('car', v=v, I=0.0, Vs=0.0, m=m, A=A, L=L, R=R, C=C, alpha=alpha, beta=beta)
    return System(elements=[car])

wheel_rail_contact(*, F=1000.0, T=10.0, width=1.0, spacing=3.0, k2=0.0, Mr=60.0, Ms=156.0)

Return a rail segment on a sleeper, loaded by the wheels of passing bogies.

The rail and the sleeper are lumped masses joined by the pad, a spring and a damper, and the sleeper rests on the ballast, Mr z'' + br (z' - zs') + kr (z - zs) = F(t) and Ms zs'' + br (zs' - z') + bs zs' + kr (zs - z) + ks zs + k2 (alpha zs^2 + β zs^3) = 0. The load is two rectangular pulses of force F per period T, one per wheel of a bogie, spacing apart. The nonlinear ballast term, with alpha = 1.75 and β = -1, is off for k2 = 0.

Parameters:

Name Type Description Default
F float

The force of a wheel.

 1000.0
T float

The period of the passing bogies.

 10.0
width float

The duration of the load of one wheel.

 1.0
spacing float

The time between the two wheels of a bogie.

 3.0
k2 float

The stiffness of the nonlinear ballast term.

 0.0
Mr float

The mass of the rail segment.

 60.0
Ms float

The mass of the sleeper.

 156.0

Returns:

Name Type Description
system System

The system, with the deflections rail.z of the rail and rail.zs of the sleeper, and their velocities.

Examples:

>>> from skcomplex.systems.engineering.power_systems import wheel_rail_contact
>>> wheel_rail_contact(k2=500.0).elements[0].k2
500.0
Source code in src/skcomplex/systems/engineering/power_systems/_electromechanics.py 
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
@register_system(
    field='engineering/power_systems',
    tags=('energy', 'vibration', 'railways', 'pulses', 'ode'),
    references=('Trzaska 2023',),
    defaults={'until': 60.0, 'dt': 0.05},
)
def wheel_rail_contact(
    *,
    F: float = 1000.0,
    T: float = 10.0,
    width: float = 1.0,
    spacing: float = 3.0,
    k2: float = 0.0,
    Mr: float = 60.0,
    Ms: float = 156.0,
) -> System:
    """Return a rail segment on a sleeper, loaded by the wheels of passing bogies.

    The rail and the sleeper are lumped masses joined by the pad, a spring and a damper, and the sleeper rests on the
    ballast, `Mr z'' + br (z' - zs') + kr (z - zs) = F(t)` and
    `Ms zs'' + br (zs' - z') + bs zs' + kr (zs - z) + ks zs + k2 (alpha zs^2 + β zs^3) = 0`. The load is two rectangular
    pulses of force `F` per period `T`, one per wheel of a bogie, `spacing` apart. The nonlinear ballast term, with
    `alpha = 1.75` and `β = -1`, is off for `k2 = 0`.

    Args:
        F:
            The force of a wheel.

        T:
            The period of the passing bogies.

        width:
            The duration of the load of one wheel.

        spacing:
            The time between the two wheels of a bogie.

        k2:
            The stiffness of the nonlinear ballast term.

        Mr:
            The mass of the rail segment.

        Ms:
            The mass of the sleeper.

    Returns:
        system:
            The system, with the deflections `rail.z` of the rail and `rail.zs` of the sleeper, and their velocities.

    Examples:
        >>> from skcomplex.systems.engineering.power_systems import wheel_rail_contact
        >>> wheel_rail_contact(k2=500.0).elements[0].k2
        500.0
    """
    rail = RailSegmentElement('rail', F=F, T=T, width=width, spacing=spacing, k2=k2, Mr=Mr, Ms=Ms)
    pad = SpringDamperInteraction(
        'pad',
        target='rail',
        states=('z', 'zs'),
        stiffness='kr',
        damping='br',
        divisor={'z': 'Mr', 'zs': 'Ms'},
    )
    return System(elements=[rail], interactions=[pad])