I Modeling a ferromagnetic core moving inside a solenoid with scipy+FEMM

[mkaluza]

TL;DR Summary

My model does not conserve energy and I have no idea what's wrong...

Hi,
I know, not another one, but I'm a more practical/engineering than scientific type - I need some tangible thing to learn stuff, so bear with me :)

Let me say hello first :)

The general idea was to model movement and electrical stuff with diff equations and do the hard part with femm. The thing is that since femm can't do transient calculations, I decided to run multiple static simulation for various projectile positions and various currents and store force and inductance to interpolate them later. This way I got F(x, i) and L(x, i) that I'm later using in scipy.
I get them from FEMM like so:
fz=femm.mo_blockintegral(19)
val1, val2, val3 = femm.mo_getcircuitproperties(self.name)
CoilResistance = val2/val1 # Calculate the ohmic resistance of the coil
CoilInductance = val3/val1 # Calculate the inductance of the coil


The electrical part of the model is just a charged capacitor C discharged into a coil and a serial diode to prevent the current from oscilating.

The problem I have is that the model does not conserve energy, so sth is obviously wrong, but I don't know what - either my di/dt equation or the interpolation, or... ???

vars: x (distance between coil and core centers), v (speed), u (cap voltage), i (current)
equations:
dx/dt = v
dv/dt = F(x, i)/m
du/dt = -i/C

To get di/dt I did:
u = iR + (d/dt)[L(x, i) * i]
u = iR + L(x, i)*di/dt + i * (d/dt)L(x, i)
but since both x and i are time dependend, I then did:
u = iR + L(x, i)*di/dt + i * [(δ/δx)L(x, i) * dx/dt + (δ/δi)L(x, i) * di/dt]
u = iR + L(x, i)*di/dt + i * (δ/δx)L(x, i) * dx/dt + i * (δ/δi)L(x, i) * di/dt
i * (δ/δi)L(x, i) * di/dt + L(x, i)*di/dt = u - iR - i * (δ/δx)L(x, i) * dx/dt
finally:
di/dt = (u - iR - i * (δ/δx)L(x, i) * dx/dt)/(i * (δ/δi)L(x, i) + L(x, i))
is that correct?

Resulting python func then fed to solve_ivp is:

def Fsim(t, X):
x, v, u, i = X
#core movement
dx_dt = v
dv_dt = A_i_x(i, x)[0][0] #A - just precalculated F(i, x)/m
#electrical circuit
du_dt = -i/C
di_dt = (u - i*R - i*v*dL_dx(i, x)[0][0])/(L_i_x(i, x)[0][0] + i*dL_di(i, x)[0][0])
if u <= 0 and i <=0: di_dt, du_dt = 0, 0 #serial diode
Y = [dx_dt, dv_dt, du_dt, di_dt]
return Y

the interpolation I use is by scipy.interpolate (yes, x and i are reversed):
F_i_x = RectBivariateSpline(I, Z, F)
A_i_x = RectBivariateSpline(I, Z, F/m)
L_i_x = RectBivariateSpline(I, Z, L)
dL_dx = L_i_x.partial_derivative(0, 1)
dL_di = L_i_x.partial_derivative(1, 0)

Link is to a zip file with precalculated data (those take a while...)
https://drive.google.com/drive/folders/1ACgUDOIAxobmZePrDnr4ZCaKOkDH52sO?usp=sharing

Complete code and sample charts attached.

Any suggestions will be very welcome :)

 

[Juanda]

Are you solving the differential equation numerically?
Are you using your own code to solve the equation or are you using the predefined libraries?
Numerical methods provide approximations that violate the conservation of energy. This violation is more or less evident depending on how good the approximation is.

Here is a related video that explains what I'm trying to say.


Your problem could be another though. I'm not experienced enough to judge only from what you wrote.