Magnetic hysteresis modelling using the Flatley differential equation model

Good day,

I am working on a Cube Satellite and am trying to correctly model magnetic hysteresis. Essentially, the hysteresis rods are ferromagnetic rods that can be easily magnetized and demagnetized by an external field (Earths Geomagnetic field in this case). The purpose of them is to dampen rotational oscillations while the satellites main magnet aligns the satellite with the local geomagnetic field vector in orbit.

I am trying to apply a magnetic field in the form of a sin wave (A) currently the amplitude is 8, but my goal is to apply various sin amplitudes in order to get multiple transient hysteresis loops. I have attached the flatley equations and what the output plots (with varyious sin amplitudes applied) should look like. I've been working on this forever with no luck, I thought it might be time to reach out. Thanks!

Bs = 0.73;         %Saturation
Br = 0.35;         %Remenance
Hc = 1.59;         %Coercivity
q = 0;
p = 2;

%%Time specifications:
SamplesPs = 1000;                   % samples per second
dt = 1/SamplesPs;                   % seconds per sample
StopTime = 2;             % sec
t = (0:dt:StopTime-dt)';     % sec
%%Sine wave:
Fc = 1;                     % Hz
H = 8*sin(2*pi*Fc*t);  %% **********(A)**********
dH = diff(H);
dHdt = dH/dt;

H(2000) = [];
t(2000) = [];

K = (1/Hc)*tan((pi/2)*(Br/Bs)) ;
i = 0;

while i <= length(H)-1
i = i + 1;  

if dHdt(i) >= 0 
B = (2*Bs/pi)*atan(K*(H(i)-Hc)); 
HL = (tan((pi*B)/(2*Bs))/K)-Hc;
Bp = (2*K*Bs/pi).*cos((pi*B)/(2*Bs)).^2;
f = (H(i)-HL)./(2*Hc); 
Q = q+(1-q)*f.^p;
dBdH = Q*Bp;
dBdt = dBdH*dHdt; 
else
B = (2*Bs/pi)*atan(K*(H(i)+Hc));   
HR = (tan((pi*B)/(2*Bs))/K)+Hc;
Bp = (2*K*Bs/pi).*cos((pi*B)/(2*Bs)).^2;
f = (H(i)-HR)./(2*Hc); 
Q = q+(1-q)*f.^p;
dBdH = Q*Bp;
dBdt = dBdH*dHdt;
B1 = trapz();
end
B_ = int(dBdt,t(i),t(i+1));
b(i) = B;
bb(i) = B_;
end
%%H is integration limit

%%Plotting Flatley Model
hold on
grid on
xlabel('Magnetizing Field [A/m]')
ylabel('Magnetic Flux Density [Tesla]')
plot(0,Br,'o')
text(0,Br,'Br','VerticalAlignment','top','HorizontalAlignment','left')
plot(0,Bs,'*')
text(0,Bs,'Bs','VerticalAlignment','top','HorizontalAlignment','left')
plot(Hc,0,'d')
text(Hc,0,'Hc','VerticalAlignment','top','HorizontalAlignment','left')
plot(0,-Br,'o')
text(0,-Br,'-Br','VerticalAlignment','top','HorizontalAlignment','left')
plot(0,-Bs,'*')
text(0,-Bs,'-Bs','VerticalAlignment','top','HorizontalAlignment','left')
plot(-Hc,0,'d')
text(-Hc,0,'-Hc','VerticalAlignment','top','HorizontalAlignment','left')

plot(H,b)
 
hold off

where constants, q = 0 and p = 2

[darova]

You should use ode45 since you have differential equation. Did you try?

[Alex]

I tried that in the past with no luck. It is possible I am using it wrong.

Answer 1

[darova]

Here is a start

function main
% code
% constants
[t,B] = ode45(f,time,B0);
plot(t,B)

    function dB = f(t,B)
        dH1 = interp1(time,dHdt,t);     % current dH        
        H1 = interp1(time,H,t);         % current H
        if dh > 0
            Hc1 = -Hc;
        else
            Hc1 = Hc;
        end
        C1 = 2*k*Bs/pi*cos(pi*B/2/Bs)^2;
        dB = 1/Hc*(H1 - 1/h*tan(pi*B/2/Bs)+Hc1)^p * C1 * dH1
    end
end

[Alex]

Thank you. What is the difference between t and time here?

[Alex]

in the dH1 and H1 calculations specifically

[darova]

I don't understand the question. time and t are different variable

time - is your time span ([0 2] i think)

t - is time point (current timestep of integration)

[Alex]

Sorry, I'll clarify. I am getting the following error within the 'f' function:

-

Not enough input arguments.

Error in main/f (line 23)

dH1 = interp1(time,dHdt,t); % current dH

Error in main (line 19)

[t,B] = ode45(f,time,B0);

-

Here is the adapted code from your response in which the error is occuring:

function main
Bs = 0.73;         %Saturation
Br = 0.35;         %Remenancewhere
Hc = 1.59;         %Coercivity
p = 2;
B0 = 0;
%%Time specifications:
SamplesPs = 1000;                   % samples per second
dt = 1/SamplesPs;                   % seconds per sample
time1 = 2;             % sec
t1 = (0:dt:time1);     % sec   to calculate magnetizing field sine wave
time = [1 2];
%%Sine wave:
Fc = 1;                     % Hz
H = 8*sin(2*pi*Fc*t1);      %Magnetizing field
dH = diff(H);
dHdt = dH/dt;

[t,B] = ode45(f,time,B0);

function dB = f(t,B)
        dH1 = interp1(time,dHdt,t);     % current dH        
        H1 = interp1(time,H,t);         % current H
        if dh > 0
            Hc1 = -Hc;
        else
            Hc1 = Hc;
        end
        C1 = 2*k*Bs/pi*cos(pi*B/2/Bs)^2;
        dB = 1/Hc*(H1 - 1/h*tan(pi*B/2/Bs)+Hc1)^p * C1 * dH1
end
end

[darova]

Try these modifications

[Alex]

Works! Thank you!

[Timothy Kuhamba]

The h on the equation must be a k dB = 1/(2*Hc)*(H1 - 1/k*tan(pi*B/2/Bs)+Hc1)^p * C1 * dH1.However the graph that am getting is different from the hystersis loop .if i make the changes proposed .

function main
Bs = 0.73;         %Saturation
Br = 0.35;         %Remenance
Hc = 1.59;         %Coercivity
p = 2;
B0 = 0;

k = (1/Hc)*tan((pi/2)*(Br/Bs))% shaping factor 
%%Time specifications:
SamplesPs = 1000;                   % samples per second
dt = 1/SamplesPs;                   % seconds per sample
time1 = 2;             % sec
t1 = (0:dt:time1);     % sec   to calculate magnetizing field sine wave
time = [1 2];
%%Sine wave:
Fc = 1;                     % Hz

H = 8*sin(2*pi*Fc*t1);      %Magnetizing field
dH = diff(H);
dHdt = dH/dt;
[t,B] = ode45(@f,time,B0);
plot(t,B)
function dB = f(t,B)
     
        dH1 = interp1(t1(2:end),dHdt,t);     % current dH   

        H1 = interp1(t1,H,t);   
        if dH1 > 0
            Hc1 = -Hc;
        else
            Hc1 = Hc;
        end
        C1 = 2*k*Bs/pi*cos(pi*B/2/Bs)^2;
        dB = 1/(2*Hc)*(H1 - 1/k*tan(pi*B/2/Bs)+Hc1)^p * C1 * dH1
end
end