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