Measurement of hysteresis losses using LTspice

Question

I have modeled the B-H curve of a random transformer. I would like to calculate the hysteresis loss over several periods by calculating the area of the curve. Let's assume that the period T is 20 ms. I know that the power for one period is calculated by the formula P=V *f *Area , where Area is the area of the loop B(H), V is the volume of the core, f is the frequency. $$P = VfArea=Vf\int_0^{20ms}{b(t)h'(t)\mathrm{d}t}$$ But how to calculate the power for several periods? I assume that the power over three periods (T=60ms) is calculated by the formula: $$P = \frac{Vf\int_0^{60ms}{b(t)h'(t)\mathrm{d}t}}{3}$$ where 3 is the number of remagnetization cycles.

I check the result in the LTspice program, where it is possible to programmatically calculate the average power for any simulation interval. Please correct me if I'm wrong. A similar topic was discussed earlier: Link

--------

SHEET 1 3552 1384
WIRE 128 64 32 64
WIRE 160 64 128 64
WIRE 320 64 240 64
WIRE 400 64 320 64
WIRE 720 64 528 64
WIRE 848 64 720 64
WIRE 1200 64 1040 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 320 304 320 208
WIRE 320 304 32 304
WIRE 400 304 400 192
WIRE 400 304 320 304
WIRE 528 304 528 208
WIRE 528 304 400 304
WIRE 720 304 720 208
WIRE 720 304 528 304
WIRE 848 304 848 208
WIRE 848 304 720 304
WIRE 1040 304 1040 208
WIRE 1040 304 848 304
WIRE 1200 304 1200 208
WIRE 1200 304 1040 304
WIRE 32 320 32 304
WIRE 1264 416 1248 416
WIRE 96 432 32 432
WIRE 592 432 528 432
WIRE 1248 432 1248 416
WIRE 528 464 528 432
WIRE 32 480 32 432
WIRE 528 560 528 544
WIRE 32 576 32 560
WIRE 96 624 32 624
WIRE 544 640 528 640
WIRE 1280 640 1248 640
WIRE 32 656 32 624
WIRE 528 656 528 640
WIRE 1248 656 1248 640
WIRE 32 752 32 736
FLAG 848 64 B_Flux_density
FLAG 1200 64 H_magnetizing_force
FLAG 400 64 L_inductance_value
FLAG 32 576 0
FLAG 96 432 H_derivative
FLAG 528 560 0
FLAG 32 752 0
FLAG 96 624 Energy_loss
FLAG 592 432 Mult
FLAG 32 320 0
FLAG 528 736 0
FLAG 544 640 Loss_3_Period
FLAG 1248 736 0
FLAG 1280 640 Cycle
FLAG 1248 512 0
FLAG 1264 416 Volume
FLAG 128 64 in
SYMBOL ind 304 112 R0
WINDOW 3 -357 255 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 -218 -37 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 32 464 R0
SYMATTR InstName B2
SYMATTR Value V=ddt(V(H_magnetizing_force))
SYMBOL bv 528 448 R0
SYMATTR InstName B3
SYMATTR Value V=(V(B_Flux_density)*V(H_derivative))
SYMBOL bv 32 640 R0
SYMATTR InstName B4
SYMATTR Value V=idt(V(Mult))
SYMBOL res 256 48 R90
WINDOW 0 0 56 VBottom 2
WINDOW 3 32 56 VTop 2
SYMATTR InstName R4
SYMATTR Value 10m
SYMBOL voltage 32 48 R0
WINDOW 3 -325 66 Left 2
WINDOW 123 0 0 Left 0
WINDOW 39 0 0 Left 0
WINDOW 0 10 5 Left 0
SYMATTR Value SINE(0 1 {f} 0 0 0 1)
SYMATTR InstName V2
SYMBOL voltage 32 128 R0
WINDOW 3 -398 56 Left 2
WINDOW 123 0 0 Left 0
WINDOW 39 0 0 Left 0
WINDOW 0 10 5 Left 0
SYMATTR Value SINE(0 0.5 {f} {1/f} 0 0 1)
SYMATTR InstName V1
SYMBOL voltage 32 208 R0
WINDOW 3 -411 51 Left 2
WINDOW 123 0 0 Left 0
WINDOW 39 0 0 Left 0
WINDOW 0 10 5 Left 0
SYMATTR Value SINE(0 0.25 {f} {2/f} 0 0 1)
SYMATTR InstName V3
SYMBOL bv 528 640 R0
WINDOW 0 38 50 Left 2
SYMATTR InstName B19
SYMATTR Value V={(V(Energy_loss)*V(Volume)*f)/V(cycle)}
SYMBOL bv 1248 640 R0
WINDOW 0 38 50 Left 2
SYMATTR InstName B20
SYMATTR Value V={{stop_time}*f}
SYMBOL bv 1248 416 R0
WINDOW 0 38 50 Left 2
SYMATTR InstName B5
SYMATTR Value V={A*Lm}
TEXT 416 -56 Left 2 !.param Hc = 42 Bs = 1.896 Br = 1.4 A = 31u Lm = 0.235 Lg = 0.1u N = 160
TEXT 416 -16 Left 2 !.param Le = {(Lm+Lg)**2/(Lm+Lg)}
TEXT -240 24 Left 2 !.tran {stop_time}
TEXT -224 -48 Left 2 !.param f=50
TEXT 728 640 Left 2 ;[W]
TEXT 1376 416 Left 2 ;[m3]
TEXT -232 -16 Left 2 !.param stop_time={10/f}