Nonlinear capacitive transducer driven by a current source

A mass hangs on a spring and damper between the plates of a capacitor. The capacitor's charge depends nonlinearly on its voltage and on the position of the mass: q = a·λ̇^1/3 + b·λ̇·(x − c)⁴. A current source and a resistor feed the capacitor, and gravity acts on the mass. The paper uses the flux (node-voltage) formulation with the electric coenergy. It derives correct equations of motion, and its simulations reproduce exactly. Three things are missing. It skips the steady state. It does not explain why b = 0.005 "diverges": the electric force grows with y³, so above a critical voltage no equilibrium exists. And it does not explain why the simulation would not start at zero voltage: the differential capacitance is infinite there. This page adds all three, a robust charge-based formulation, and an interactive simulator.

Problem statement

The transducer in the figure is described by a characteristic equation. It relates the charge $q$ on the capacitor to the voltage $\dot\lambda$ and the position coordinate $x$:

$$q = a\,\dot\lambda^{1/3} + b\,\dot\lambda\,(x-c)^4$$

When $x = c$, the mass $M$ is centred between the capacitor plates and the spring is relaxed. Find the equations of motion.

The mass $M$ moves vertically between the plates. It is held by a spring $k$ and a damper $B$ to the ground, and gravity $g$ acts on it.
The capacitor and the resistor $R$ are connected in parallel to a current source $i(t)$. The voltage of the upper node is $\dot\lambda$, where $\lambda$ is the flux linkage.
The paper uses $y = c - x$, the downward displacement from the centred position, so the term $(x-c)^4$ becomes $y^4$.

The charge law is a model written for this exercise, not a real geometry. Its first term is a nonlinear capacitor that does not depend on position. Its second term is a linear capacitance $by^4$ that grows steeply as the mass leaves the centre.

Figure 1. The mass M between the capacitor plates rests on a spring k and damper B. The capacitor and resistor R are in parallel with the current source i(t). The node voltage is e = dλ/dt. The coordinate y = c − x points down.

Generalized coordinates and energies

With a current source and node voltages, the natural electrical coordinate is the flux $\lambda$, whose derivative is the node voltage $e = \dot\lambda$. The circuit has two nodes, and one is the reference, so there is one flux coordinate. The mechanical coordinate is $y$. In this formulation the capacitor contributes its coenergy, the integral of the charge over the voltage:

$$W_e^* = \int_0^{\dot\lambda} q\,d\dot\lambda = \tfrac34\,a\,\dot\lambda^{4/3} + \tfrac12\,b\,y^4\dot\lambda^2$$

TerminologyAs in the course: repulsive coenergy $\mathcal{R}^* = \tfrac12M\dot y^2$, attractive energy $\mathcal{A} = \tfrac12ky^2 - Mgy$ (spring plus gravity, since $y$ points down), electric coenergy $W_e^*$, magnetic energy $W_m = 0$ (no coil). The Lagrangian is $\mathcal{L} = \mathcal{R}^* + W_e^* - \mathcal{A} - W_m$.

$$\mathcal{L} = \tfrac12M\dot y^2 + \tfrac34a\dot\lambda^{4/3} + \tfrac12by^4\dot\lambda^2 - \tfrac12ky^2 + Mgy,\qquad \mathcal{F} = \tfrac12B\dot y^2 + \frac{\dot\lambda^2}{2R}$$

The current source is the generalized force of $\lambda$. There is no external force on $y$, because gravity is already in $\mathcal{A}$. The paper sets up all of this correctly. The electric co-Rayleigh function is written $\int d\dot\lambda/R$, which should read $\int(\dot\lambda/R)\,d\dot\lambda$, but its result is right.

Equations of motion

(1)$$M\ddot y + B\dot y + ky - 2b\,\dot\lambda^2y^3 = Mg$$

(2)$$\Big(\frac a3\dot\lambda^{-2/3} + by^4\Big)\ddot\lambda + \Big(4by^3\dot y + \frac1R\Big)\dot\lambda = i(t)$$

Equation (2) is Kirchhoff's current law at the node: $\dot q + e/R = i$. The first bracket is the differential capacitance $\partial q/\partial e$. The term $4by^3\dot y\,e$ is the motional current: moving the mass changes the capacitance at a given voltage.
The electric force $F_e = \partial W_e^*/\partial y = 2b\,e^2y^3$ has the sign of $y$. It pushes the mass away from the centre, because the capacitance $by^4$ grows in both directions. It therefore works against the spring, a negative stiffness that grows with $y^2$.
$\lambda$ itself never appears, only $\dot\lambda$ and $\ddot\lambda$. It is a cyclic coordinate, so the electrical part is really first order in $e = \dot\lambda$. The paper's second integrator ($\dot\lambda\to\lambda$) only produces an output that nothing uses.

Both equations agree with the paper's (*) and (**). One intermediate step, (2.21), has a stray "$2yk/s$", but the final result is correct.

Steady state and its limit (missing in the paper)

The table of contents lists section 2.7, "Steady-state equations", with a question mark instead of a page number, and the section is missing. For a constant current $i_0$, all derivatives vanish. Equation (2) then gives the node voltage and (1) the position:

(3)$$e^* = R\,i_0,\qquad ky^* - \beta\,y^{*3} = Mg,\qquad \beta = 2b\,e^{*2}$$

For the base parameters ($R = 10\ \Omega$, $i_0 = 1$ A, $b = 5\cdot10^{-5}$) this gives $e^* = 10$ V, $\beta = 0.01$ and $y^* = 5.03$ m. Without the electric force it would be $Mg/k = 4.905$ m. This explains the paper's observation 3°: a larger $R$ means a larger voltage, a larger outward force and a larger $y^*$. As $R\to0$, the position tends to the purely mechanical value.

When does an equilibrium exist?

The left side $ky - \beta y^3$ of (3) rises to a maximum and then falls. The maximum is at $y_m = \sqrt{k/(3\beta)}$ and equals $\tfrac23ky_m$. An equilibrium exists only if gravity does not exceed it:

(4)$$Mg \le \frac{2k}{3}\sqrt{\frac{k}{3\beta}}\quad\Longleftrightarrow\quad 2b\,(Ri_0)^2 \le \frac{4k^3}{27\,(Mg)^2}$$

Critical resistance. For the base parameters, (4) gives $R_{cr} = 24.8\ \Omega$. $R = 24\ \Omega$ settles towards $y^* = 6.45$ m (slowly, since $k_{ef}$ has dropped to 2.8 N/m), but $R = 26\ \Omega$ runs away.
Critical b. At $R = 10\ \Omega$, $b_{cr} = 3.1\cdot10^{-4}$. The paper's $b = 0.005$ is 16 times larger, so no equilibrium exists. That is the true reason for its "divergence" (1°, 7°), and why a smaller $R$ "mitigates" it. The motion is slow, because the capacitance $by^4$ becomes enormous and the voltage collapses to about 2 V. After 300 s the mass is at 12 m and still moving.
Stiffness. The stable equilibrium has an effective stiffness $k_{ef} = k - 3\beta y^{*2}$, which falls to zero at the limit (4). This is the same pull-in mechanism as in problems 17, 18 and 20.

Simulink model and results

The Simulink model integrates (1) and (2) with three subsystems for the nonlinear terms. A unit step drives the current source. The base parameters are $g = 9.81$, $k = 10$, $M = 5$, $b = 5\cdot10^{-5}$, $a = 0.05$, $B = 5$ and $R = 10$. Our simulation, with the robust charge formulation of section 6, reproduces every plot:

Case	Paper	Simulation
Base	y peaks ≈ 6.4 m at 2.4 s, settles ≈ 5; voltage dips to ≈ 6 V, peaks ≈ 13.5 V, settles at 10 V	6.45 m at 2.40 s → 5.03 m; 5.6 / 13.5 V → 10 V ✓
b = 0.005	"Diverges"; no steady state even after 100 000 s	No equilibrium exists (section 4); slow runaway ✓
R = 1 / 20 Ω	y* ≈ 4.9 / 5.6 m; voltage settles at 1 / 20 V	4.91 / 5.61 m; 1 / 20 V ✓
M = 0.5	y ≈ 0.49 m without overshoot; voltage flat at 10 V	0.491 m, overdamped (ζ = 1.1) ✓
k = 50	y ≈ 1 m, faster oscillation	Peak 1.58 m → 0.98 m ✓
B = 10	Smaller overshoot, faster settling	5.19 m peak → 5.03 m ✓
a = 0.9 / 0.005	No visible effect on y; a large a slows the voltage rise, a small a gives a jump to 10 V	Reproduced; e(t) = (t/a)³ at the start ✓

Observation 10° is also correct: with $M = 0.5$ the mass settles without oscillating, so the voltage shows no dip. The voltage dips and peaks because the mass moves. The capacitance $by^4$ changes, and with it the motional current $4by^3\dot y\,e$.

The zero-voltage trap, and a better formulation

The paper notes that the first simulation attempt failed. It attributes this to the equation containing only the derivatives of $\lambda$, and fixes it by starting the voltage integrator at 0.001. The real cause is the first term of the differential capacitance:

$$C_d(e) = \frac{\partial q}{\partial e} = \frac a3\,e^{-2/3} + by^4 \;\to\; \infty\quad\text{as } e\to0$$

At $e = 0$ the capacitor looks infinitely large. Equation (2) gives $\dot e = i/C_d = 0$, and the voltage stays at zero forever. Physically the current still charges the capacitor, $q = a\,e^{1/3} = i_0t$, so $e = (i_0t/a)^3$. Both $e \equiv 0$ and this solution satisfy (2), because the right-hand side $\propto e^{2/3}$ is not Lipschitz at zero. The solver picks the wrong one. The preset "Paper form, e(0) = 0" shows the mass simply sagging to $Mg/k$ with no voltage at all.
Negative voltages (for example a current reversal) give complex numbers in Simulink's $u^v$ block. The real cube root has to be used.

The fix: use the charge as the state. The conserved quantity of the cyclic coordinate $\lambda$ is $\partial\mathcal{L}/\partial\dot\lambda = q$, so (2) is simply

(5)$$\dot q = i(t) - \frac eR,\qquad e \text{ from } q = a\,e^{1/3} + by^4e$$

With $s = e^{1/3}$ the constitutive relation is the cubic $by^4s^3 + as - q = 0$. It is strictly increasing in $s$, so it has exactly one real root, found in a few Newton steps. The state $(y, \dot y, q)$ starts from zero without any trick, handles negative voltages, and never divides by an infinite capacitance. The simulator uses (5) by default and lets you switch to the paper's form for comparison.

Critical review

Missing steady state and stability analysis

Section 2.7 is announced but missing. Equation (3) would have predicted all the settled values in the plots. The existence condition (4) would have explained the "divergence" for $b = 0.005$ and the role of $R$ quantitatively, instead of leaving them as observations.

The direction of the electric force is not discussed

The term $-2b\dot\lambda^2y^3$ in (1) is a destabilizing, outward force, a negative stiffness. The paper never says so, although this one fact explains observations 1°, 3°, 7° and 9°.

The numerical workaround is misdiagnosed

The start-up failure is not caused by $\lambda$ appearing only through its derivatives. It is caused by the infinite differential capacitance at zero voltage, which makes $e \equiv 0$ a spurious solution. Setting $e(0) = 0.001$ works, but it hides the problem. A charge-based formulation (5) removes it, and the unused integrator for $\lambda$ can be dropped.

Physical scale and smaller points

With $M = 5$ kg and $k = 10$ N/m the mass sags almost 5 m. It would leave any real capacitor long before that, and the polynomial charge law cannot hold there. Lighter or stiffer parameters ($M = 0.5$ or $k = 50$) give more sensible motions. Several smaller points:

Observation 8° contradicts itself: "decreasing a gives a slow rise … decreasing a gives a jump". It should say that increasing $a$ slows the rise.
One axis is labelled "time t [t]".
The plot legends are missing, as the paper's footnote admits, so the two curves in each comparison cannot be told apart.
The conclusion is generic and quotes none of the results.

What is goodThe flux formulation with the electric coenergy is the right choice for a current source and a voltage-controlled charge law. The coenergy integral, the Lagrangian with gravity in the attractive energy, the Rayleigh functions and both equations of motion are correct. The Simulink model is reproduced to the plotted accuracy. The parameter study is systematic and its observations are correct, even though they are not explained.

Interactive simulator

The simulator integrates (1) with either the charge form (5) or the paper's voltage form (2), using 4th-order Runge–Kutta. The current is a step $i_0$ at $t = 0$. The last chart is the static picture behind (3)–(4): equilibria are where the curve $ky - \beta y^3$ meets the line $Mg$.

Cases

Electrical equation

Parameters

k [N/m]

M [kg]

B [N·s/m]

g [m/s²]

a [C·V⁻¹ᐟ³]

b [F/m⁴]

R [Ω]

i₀ [A] (step)

e(0) [V] (paper form)

Duration [s]

Displacement y(t)

Voltage e = dλ/dt

Forces on the mass

Static equilibrium at e* = R·i₀

Equilibria lie where the curve crosses Mg.

Animation

Results table

Review questions

Why is the electric coenergy, not the energy, used here?

The coordinate is the flux $\lambda$, so the voltage $\dot\lambda$ is a generalized velocity. The charge law gives $q$ as a function of the voltage. The capacitor is then described by the coenergy $W_e^*(\dot\lambda, y) = \int q\,d\dot\lambda$, which enters the Lagrangian like a kinetic term.

Why does the electric force push the mass away from the centre?

At constant voltage the force is $\partial W_e^*/\partial y = 2be^2y^3$. The capacitance $by^4$ grows when the mass moves away from the centre in either direction, so the system is pulled towards larger capacitance, which is outward.

Why does b = 0.005 never settle?

Then $2b(Ri_0)^2 = 1$, and the maximum of $ky - \beta y^3$ is only 12.2 N, far below $Mg = 49$ N. No equilibrium exists, and the mass keeps sinking. It sinks slowly, because the huge capacitance keeps the voltage low.

Why can a simulation started at exactly zero voltage stay at zero?

The differential capacitance $\tfrac a3e^{-2/3}$ is infinite at $e = 0$, so (2) gives $\dot e = 0$. The equation $\dot e \propto e^{2/3}$ has two solutions from zero, and the trivial one is chosen. Using the charge as the state avoids this.

What does the resistor do in the steady state?

All of the source current flows through it, because a capacitor carries no DC current. The voltage is therefore $e^* = Ri_0$. It sets the electric force and with it $y^*$, and through the limit (4) whether an equilibrium exists at all.