- #1

mkaluza

- 1

- 1

- 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 :)

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 :)