As you found out yourself, parametric plots can't be integrated, but there is a way to do it. Mathematically, the Y axis would be \$b(t)\$ and the X axis \$h(t)\$, and they're both functions of time. To find out the area uner a curve, you take the integral over an interval, but since \$\mathrm{d}t\$ is now a function, itself a function of \$t\$, it all becomes:
$$A=\int_0^T{b(t)\mathrm{d}h(t)}=\int_0^T{b(t)h'(t)\mathrm{d}t}\tag{1}$$
Therefore you have to take the derivative of H and multiply that with B in order to have in integrated (without the need for a parametric plot):
In the plot, one square is 0.3*10=3 units, and a crude approximation count is about 18*3*2=90. It's clear that I'm underselling, but whether there are 4 more... However it is, B1 does what's described above, with V(h') being the derivative, and the rest is a way to apply definite integration. The .meas script shows -111.5, which is greater than the crude approximation by about 20% (which is better than 50%, I suppose), and it's a negative value because the product turned out this way, so it's just a matter of sign.
Word of caution: because you are using a derivative, if H is wild then the derivative will be even wilder, so don't forget to check the signals before pointing a finger.
Or you could try to use the formulas described in the help, under LTspice > Circuit Elements > L. Inductor (though they'll need to be used as a parametric plot).
Because it's difficult to determine the accuracy of the results above, it's never a bad idea to use a simpler approach. For example [Digital]/schmitt* for the hysteresis, and for the input a sine, since it needs to be differentiated. The plot shows the hysteresis window as being from -1 V to +1 V on the Y-axis, thus 2 V tall, and from -0.5 V to +0.5 V on the X-axis, thus 1 V wide, making its area 2. Just like above, G1+L1 differentiates the input giving V(z), which is then multiplied by the input, V(x), to be integrated (continuously) over a period (T1 with Td=1, since V(x) is a sine with period of 1), giving V(A) as the result, shown with a .meas script.
The output is still negative, because the multiplication starts with a negative number. The value is what matters.
I started this as a comment, but it looks like it will take a bit more.
Don't you think that the problem becomes from the mathematical definition that we took at the beginning of your post ?
Not at all. As I said in a previous comment, by not imposing initial conditions you are relying on LTspice to find the appropriate operating point for all the elements, which happens to be non-zero for the initial time. This is because when there are no initial conditions, the SPICE engine considers that the circuit has been running since the dawn of time, and when the user presses "run" what is seen is the steady-state as determined by the solver.
Given your circuit: I1 starts with a 100 A. This, alone, is an abomination if you consider that everything should have started at zero, but let's follow this thread. The solver considers that the circuit has been running since the Big Bang so when "run" is clicked, the whole circuit must behave as if I1 and the rest of the circuit has reached the 100 A value in a decreasing ramp. This ramp feeds an inductor, whose voltage feeds an integrator and its current a voltage source, and other mathematical operations based on these quantities. Therefore the solver must figure out what happened until the user clicks "run" such that at t=0 everything should be as close as possible to what exists in the schematic.
Because of the above, the integrator doesn't start at zero, but that's no reason why the formula is wrong. Hmm, looking over the whole answer and the coments, I see I forgot to mention this. I'm sorry. When performing indefinite integration and you're interested in the result over one period, you have to treat it the same way you would treat an integration with initial conditions: subtract the value at the beginning of the interval from the one at the end of it.
The same thing needs to be done here: V(Energy_loss) starts from 7.0443931 kV (as shown by the cursors) and wobbles to the final value of 7.0444878 kV, which means the difference is 94.726563 mV. It makes sense since your hysteresis has non-zero width. If you duplicate B4 and modify its expression to contain my other suggestion in the comments, idt(V(b), delay(V(b), 400)), you'll see that it settles at the final value of 184.24845 mV. Which is acceptable given the imposed timestep that the waveform compression is on.
I also hinted at behavioural sources vs primitives in the comments, so if you'll use the circuit shown below then this is what you'll get:
Take a bit of time to see what is shown:
- Top right there's the error log showing the results of the
.meas scripts seen in the bottom-right. The difference between what B4 (not shown but present in the schematic as giving V(Energy_loss)) is showing as ve=ve2-ve1=-0.41 (with the same timestep and .opt), vs V(test2), which uses the derivative as given by G2+L2, as vt=vt2-vt1=0.02, is more than an order of magnitude. And this is considering that the quantities are kV (don't forget that this is a SPICE engine, which deals in tolerances, not precision).
- Plotted are
V(energy_loss_2), which is my suggestion of a definite integral with a behavioural source (the idt(x-delay(y)) thing), vs. V(test) which is the version made out with primitives (and a B source). The readings are with the cursors shown in the top left, and they show 492 μV for the first and 6.8 nV for the second. Now use the parametric plot for the B-H curve and zoom in around zero. A simple click-drag rectangle will show about 9.6 mV on the X-axis and about 0.8 μV on the Y-axis. Multiplied they should give nV.
I conclude that the formula is very much valid, but also that it matters how it's implemented. If you want to see for yourself, here is the new source (sorry about the way they're thrown in, I just separated them for the screenshot):
Version 4
SHEET 1 2460 1380
WIRE 320 64 144 64
WIRE 400 64 320 64
WIRE 720 64 528 64
WIRE 848 64 720 64
WIRE 1200 64 1040 64
WIRE 144 128 144 64
WIRE 320 128 320 64
WIRE 528 128 528 64
WIRE 720 128 720 64
WIRE 1040 128 1040 64
WIRE 1200 128 1200 64
WIRE 400 144 400 64
WIRE 480 144 400 144
WIRE 848 144 848 64
WIRE 480 192 400 192
WIRE 144 272 144 208
WIRE 320 272 320 208
WIRE 320 272 144 272
WIRE 400 272 400 192
WIRE 400 272 320 272
WIRE 528 272 528 208
WIRE 528 272 400 272
WIRE 720 272 720 208
WIRE 720 272 528 272
WIRE 848 272 848 208
WIRE 848 272 720 272
WIRE 1040 272 1040 208
WIRE 1040 272 848 272
WIRE 1200 272 1200 208
WIRE 1200 272 1040 272
WIRE 144 288 144 272
WIRE 320 464 144 464
WIRE 864 464 720 464
WIRE 1632 464 1360 464
WIRE 144 512 144 464
WIRE 1360 512 1360 464
WIRE 720 528 720 464
WIRE 144 656 144 592
WIRE 720 656 720 608
WIRE 1360 656 1360 592
WIRE 624 800 560 800
WIRE 560 848 560 800
WIRE 240 864 176 864
WIRE 336 864 240 864
WIRE 384 864 336 864
WIRE 560 960 560 928
WIRE 592 1008 512 1008
WIRE 672 1008 592 1008
WIRE 464 1024 128 1024
WIRE 128 1072 128 1024
WIRE 224 1072 128 1072
WIRE 384 1072 320 1072
WIRE 464 1072 384 1072
WIRE 192 1248 128 1248
FLAG 144 288 0
FLAG 848 64 B_Flux_density
FLAG 1200 64 H_magnetizing_force
FLAG 400 64 L_inductance_value
FLAG 144 656 0
FLAG 320 464 H_derivative
FLAG 720 656 0
FLAG 1360 656 0
FLAG 1632 464 Energy_loss
FLAG 864 464 b
FLAG 560 960 0
FLAG 624 800 Energy_loss_2
FLAG 128 880 H_magnetizing_force
FLAG 128 928 0
FLAG 176 944 0
FLAG 240 944 0
FLAG 384 1152 0
FLAG 512 1088 0
FLAG 224 1104 0
FLAG 320 1104 0
FLAG 592 1072 0
FLAG 336 864 H`
FLAG 128 1152 0
FLAG 672 1008 test
FLAG 128 1328 0
FLAG 192 1248 test2
SYMBOL current 144 208 R180
WINDOW 0 24 80 Left 2
WINDOW 3 24 0 Left 2
WINDOW 123 0 0 Left 0
WINDOW 39 0 0 Left 0
SYMATTR InstName I1
SYMATTR Value PWL(0 0 {1*Is} {Is} {2*Is} 0 {3*Is} {-Is} {4*Is} 0 {5*Is} {Is})
SYMBOL ind 304 112 R0
WINDOW 3 -357 234 Left 2
WINDOW 39 36 108 Left 2
SYMATTR Value Hc = {Hc} Bs = {Bs} Br = {Br} A = {A} Lm = {Lm} Lg = {Lg} N = {N}
SYMATTR SpiceLine Rser = 0
SYMATTR InstName L1
SYMBOL res 704 112 R0
SYMATTR InstName R1
SYMATTR Value 100Meg
SYMBOL cap 832 144 R0
SYMATTR InstName C1
SYMATTR Value 1
SYMBOL bi 1040 208 R180
WINDOW 0 24 80 Left 2
WINDOW 3 -203 -23 Left 2
SYMATTR InstName B1
SYMATTR Value I=I(L1)*{N/Le}
SYMBOL res 1184 112 R0
SYMATTR InstName R2
SYMATTR Value 1
SYMBOL g 528 112 R0
SYMATTR InstName G1
SYMATTR Value {1/(A*N)}
SYMBOL bv 144 496 R0
SYMATTR InstName B2
SYMATTR Value V=ddt(V(H_magnetizing_force))
SYMBOL bv 720 512 R0
SYMATTR InstName B3
SYMATTR Value V=(V(B_Flux_density)*V(H_derivative))
SYMBOL bv 1360 496 R0
SYMATTR InstName B4
SYMATTR Value V=idt(V(b))
SYMBOL bv 560 832 R0
SYMATTR InstName B5
SYMATTR Value V=idt(V(b)-delay(v(b),400))/400
SYMBOL g 176 848 R0
WINDOW 0 12 -10 Left 2
SYMATTR InstName G2
SYMATTR Value 1
SYMBOL ind 224 848 R0
SYMATTR InstName L2
SYMATTR Value 1 rser=0
SYMBOL tline 272 1088 R0
WINDOW 0 -5 46 Bottom 2
WINDOW 3 -1 -45 Top 2
SYMATTR InstName T1
SYMATTR Value Td=400 Z0=1
SYMBOL res 368 1056 R0
SYMATTR InstName R3
SYMATTR Value 1
SYMBOL g 512 992 R0
WINDOW 0 27 90 Left 2
WINDOW 3 29 119 Left 2
SYMATTR InstName G4
SYMATTR Value 1
SYMBOL cap 576 1008 R0
SYMATTR InstName C2
SYMATTR Value 400 rpar=1g
SYMBOL bi 128 1152 M180
WINDOW 0 29 -9 Left 2
WINDOW 3 -13 -38 Left 2
SYMATTR InstName B6
SYMATTR Value I=V(B_Flux_density)*V(H`)
SYMBOL bi 128 1328 M180
WINDOW 0 29 -9 Left 2
WINDOW 3 -13 -38 Left 2
SYMATTR InstName B7
SYMATTR Value I=V(B_Flux_density)*V(H`) Rpar=1g Cpar=1
TEXT 136 -192 Left 2 !.param Hc = 0.01 Bs = 0.51 Br = 0.01 A = 31u Lm = 47m Lg = 0.725m N = 8
TEXT 136 -128 Left 2 !.param I_origine = Is
TEXT 136 -160 Left 2 !.param Is = 100
TEXT 136 -96 Left 2 !.param Cycle = 5
TEXT 136 -64 Left 2 !.param Le = {(Lm+Lg)**2/(Lm+Lg)}
TEXT -552 -136 Left 2 ;L : 1 volt = 1 Henry = Valeur de la self\nB : 1 volt = 1 Tesla = densité d'induction magnétique\nH : 1 volt = 1A/m = champ magnétique
TEXT -552 -184 Left 2 !.tran 0 {Is*Cycle} {I_origine} {Is/1000}
TEXT 1128 -232 Left 2 !.meas self0 find V(L_inductance_value) when I(L1) = 0\n.meas self418 find V(L_inductance_value) when I(L1) = 4.18 \n.meas self10 find V(L_inductance_value) when I(L1) = 10\n.meas self20 find V(L_inductance_value) when I(L1) = 20\n.meas self30 find V(L_inductance_value) when I(L1) = 30\n.meas self50 find V(L_inductance_value) when I(L1) = 50
TEXT 1120 432 Left 2 !.meas Energy_loss param V(Energy_loss)
TEXT 608 1208 Left 2 !.meas ve1 find v(energy_loss) at 0\n.meas ve2 find v(energy_loss) at 400\n.meas ve param ve2-ve1\n.meas vt1 find v(test2) at 0\n.meas vt2 find v(test2) at 400\n.meas vt param vt2-vt1
