How to avoid the overshoot in Jiles Atherton hysteresis loop

Question

Md Tanbir Siddik 2024-4-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2111566-how-to-avoid-the-overshoot-in-jiles-atherton-hysteresis-loop

评论： Mathieu NOE 2024-5-6

Trying to plot a hysteresis model using jiles atherton hysteresis model in matlab. All the parameters are defined in the beginning.

Here is the code I have written,

clear

%% Given parameters

a = 512;

c = 0.75;

k = 570;

M_s = 1.7*10^6;

i_s = 10^-6;

alpha = 0.032;

u_0 = 4*pi*10^-7;

f = 1;

h_max = 8000;

%% Generate time vector

t_start = 0; %start time

t_end = 60; %end time

dt = 0.001; %time steps

t = t_start:dt:t_end;

%% Initial values

M_irev(1) = 0; %Initial irreversible magnetization

h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength

%% Time stepping loop

for i=2:length(t)

% Compute effective field

h_e = h(i); % Applied field and effective field considered equal

% Compute anhysteric magnetization

M_an = M_s*(coth(h_e/a) - (a/(h_e))); % Langevin function

% Compute derivative of magnetic field with respect to time

dh_e = (h(i) - h(i-1)) / dt;

delta = sign(dh_e);

% Compute irreversible magnetization using Euler method

M_irev(i) = ((M_an - M_irev(i-1)) / ((k * delta) - (alpha * (M_an - M_irev(i-1))))) * dh_e;

% Update magnetization

M(i) = c*M_an + (1-c)*M_irev(i);

end

%% Magnetic flux density

B = u_0 * (h + M);

%% Plot results

figure;

plot(h, B);

xlabel('Magnetic Field Strength');

ylabel('Magnetic Flux Density');

title('HB Graph');

grid on;

The problem I am facing is when it plots the curve, it shows overshoot at three points. Suspecting the problem is in the 'Langevin function'. Lower values of h at the end part creates bigger values thus it shows overshoot at those points. I think I need to use approximation (Taylor series or similar method) at those points so that when overshoot happens the code will avoid those overshoots and take the approximation. Can anyone show me how can I do that approximation?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2024-4-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2111566-how-to-avoid-the-overshoot-in-jiles-atherton-hysteresis-loop#answer_1447926

在 MATLAB Online 中打开

hello

I am just trying to remove the spikes (with filloutiers) and further smooth the curve with a touch of smoothdata

look at these 2 lines I added at the plot section

Bs= filloutliers(B,'linear','movmedian',15, 'ThresholdFactor', 1.1);

Bs= smoothdata(Bs,'sgolay',50);

code :

clear

%% Given parameters

a = 512;

c = 0.75;

k = 570;

M_s = 1.7*10^6;

i_s = 10^-6;

alpha = 0.032;

u_0 = 4*pi*10^-7;

f = 1;

h_max = 8000;

%% Generate time vector

t_start = 0; %start time

t_end = 60; %end time

dt = 0.001; %time steps

t = t_start:dt:t_end;

%% Initial values

M_irev(1) = 0; %Initial irreversible magnetization

h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength

%% Time stepping loop

for i=2:length(t)

% Compute effective field

h_e = h(i); % Applied field and effective field considered equal

% Compute anhysteric magnetization

M_an = M_s*(coth(h_e/a) - (a/(h_e))); % Langevin function

% Compute derivative of magnetic field with respect to time

dh_e = (h(i) - h(i-1)) / dt;

delta = sign(dh_e);

% Compute irreversible magnetization using Euler method

M_irev(i) = ((M_an - M_irev(i-1)) / ((k * delta) - (alpha * (M_an - M_irev(i-1))))) * dh_e;

% Update magnetization

M(i) = c*M_an + (1-c)*M_irev(i);

end

%% Magnetic flux density

B = u_0 * (h + M);

%% Plot results

figure;

Bs= filloutliers(B,'linear','movmedian',15, 'ThresholdFactor', 1.1);

Bs= smoothdata(Bs,'sgolay',50);

plot(h, B, h, Bs);

xlabel('Magnetic Field Strength');

ylabel('Magnetic Flux Density');

title('HB Graph');

grid on;

7 个评论
显示 5更早的评论隐藏 5更早的评论

Mathieu NOE 2024-4-25

编辑：Mathieu NOE 2024-4-25

在 MATLAB Online 中打开

we could also do a cycle by cycle averaging and get the one cycle averaged curve , but at the end the result may look quite the same (almost , stil there is a small discontinuity, we can put a bit of smoothing there)

for the time being, this is the result without smoothing :

%% Given parameters

a = 512;

c = 0.75;

k = 570;

M_s = 1.7*10^6;

i_s = 10^-6;

alpha = 0.032;

u_0 = 4*pi*10^-7;

f = 1;

h_max = 8000;

%% Generate time vector

t_start = 0; %start time

t_end = 60; %end time

dt = 0.001; %time steps

t = t_start:dt:t_end;

%% Initial values

M_irev(1) = 0; %Initial irreversible magnetization

h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength

%% Time stepping loop

for i=2:length(t)

% Compute effective field

h_e = h(i); % Applied field and effective field considered equal

% Compute anhysteric magnetization

M_an = M_s*(coth(h_e/a) - (a/(h_e))); % Langevin function

% Compute derivative of magnetic field with respect to time

dh_e = (h(i) - h(i-1)) / dt;

delta = sign(dh_e);

% Compute irreversible magnetization using Euler method

M_irev(i) = ((M_an - M_irev(i-1)) / ((k * delta) - (alpha * (M_an - M_irev(i-1))))) * dh_e;

% Update magnetization

M(i) = c*M_an + (1-c)*M_irev(i);

end

%% Magnetic flux density

B = u_0 * (h + M);

%% cycle by cycle averaging

dtzc = 1/f;

tzc = t(1):dtzc:t(end);

NN = numel(tzc)-1;

hh = 0;

BB = 0;

for k = 1:NN

tstart = tzc(k);

tstop = tzc(k+1);

tt = linspace(tstart,tstop,360); % 360 => 1 degree resolution

hi = interp1(t,h,tt);

Bi = interp1(t,B,tt);

hh = hh + hi;

BB = BB + Bi;

end

% divide by N to obtain the averaged data

hh = hh/NN;

BB = BB/NN;

%% Plot results

figure;

plot(h, B, hh, BB);

xlabel('Magnetic Field Strength');

ylabel('Magnetic Flux Density');

title('HB Graph');

grid on;

Mathieu NOE 2024-4-25

在 MATLAB Online 中打开

just added now a bit of smoothing on top

see the effect here : smoothed curve is yellow

%% Given parameters

a = 512;

c = 0.75;

k = 570;

M_s = 1.7*10^6;

i_s = 10^-6;

alpha = 0.032;

u_0 = 4*pi*10^-7;

f = 1;

h_max = 8000;

%% Generate time vector

t_start = 0; %start time

t_end = 60; %end time

dt = 0.001; %time steps

t = t_start:dt:t_end;

%% Initial values

M_irev(1) = 0; %Initial irreversible magnetization

h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength

%% Time stepping loop

for i=2:length(t)

% Compute effective field

h_e = h(i); % Applied field and effective field considered equal

% Compute anhysteric magnetization

M_an = M_s*(coth(h_e/a) - (a/(h_e))); % Langevin function

% Compute derivative of magnetic field with respect to time

dh_e = (h(i) - h(i-1)) / dt;

delta = sign(dh_e);

% Compute irreversible magnetization using Euler method

M_irev(i) = ((M_an - M_irev(i-1)) / ((k * delta) - (alpha * (M_an - M_irev(i-1))))) * dh_e;

% Update magnetization

M(i) = c*M_an + (1-c)*M_irev(i);

end

%% Magnetic flux density

B = u_0 * (h + M);

%% cycle by cycle averaging

dtzc = 1/f;

tzc = t(1):dtzc:t(end);

NN = numel(tzc)-1;

hh = 0;

BB = 0;

for k = 1:NN

tstart = tzc(k);

tstop = tzc(k+1);

tt = linspace(tstart,tstop,360); % 360 => 1 degree resolution

hi = interp1(t,h,tt);

Bi = interp1(t,B,tt);

hh = hh + hi;

BB = BB + Bi;

end

% divide by N to obtain the averaged data

hh = hh/NN;

BB = BB/NN;

Bs= smoothdata(BB,'loess',20);

%% Plot results

figure;

plot(h, B, hh, BB, hh, Bs);

xlabel('Magnetic Field Strength');

ylabel('Magnetic Flux Density');

title('HB Graph');

grid on;

Mathieu NOE 2024-5-6

在 MATLAB Online 中打开

I am not sure it's the Langevin function which causes the issue

I traced down a few quantities (so I had to mdify a bit your code

if I plot Langevin function output (M_an), there is not much issue to be seen , except one single peak)

it seems to me there is a bigger issue in the Euler equation , if you trace down the numerator and denominator, it shows a lot of spikes , so there is some discontinuities in your process that must be corrected

%% Given parameters

a = 512;

c = 0.75;

k = 570;

M_s = 1.7*10^6;

i_s = 10^-6;

alpha = 0.032;

u_0 = 4*pi*10^-7;

f = 1;

h_max = 8000;

%% Generate time vector

t_start = 0; %start time

t_end = 60; %end time

dt = 0.001; %time steps

t = t_start:dt:t_end;

%% Initial values

M_irev(1) = 0; %Initial irreversible magnetization

h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength

%% Time stepping loop

for i=2:length(t)

% Compute effective field

h_e = h(i); % Applied field and effective field considered equal

% Compute anhysteric magnetization

M_an(i) = M_s*(coth(h_e/a) - (a/h_e)); % Langevin function

% Compute derivative of magnetic field with respect to time

dh_e(i) = (h(i) - h(i-1)) / dt;

delta(i) = sign(dh_e(i));

% Compute irreversible magnetization using Euler method

num(i) = (M_an(i) - M_irev(i-1))* dh_e(i);

den(i) = ((k * delta(i)) - (alpha * (M_an(i) - M_irev(i-1))));

% M_irev(i) = ((M_an(i) - M_irev(i-1)) / ((k * delta) - (alpha * (M_an(i) - M_irev(i-1))))) * dh_e(i);

M_irev(i) = num(i) / den(i);

% Update magnetization

M(i) = c*M_an(i) + (1-c)*M_irev(i);

end

%% Magnetic flux density

B = u_0 * (h + M);

%% Plot results

figure;

plot(h, B);

% plot(h, B, hh, BB, hh, Bs);

xlabel('Magnetic Field Strength');

ylabel('Magnetic Flux Density');

title('HB Graph');

grid on;

%% more plots (for debug purposes)

figure;

plot(M_an);

title('M-an');

figure;

plot(delta);

title('delta');

figure;

plot(num);

title('M-irev num');

figure;

plot(den);

title('M-irev den');

Md Tanbir Siddik 2024-5-6

Hi,

Thank you for your continuous support. The code with interpolation and smoothing solves the problem for now. I somehow overlooked that response when I was replying. It seems when I further decrease the step size from 0.0001, the curve shows some wave at those overshooting points. which later creates larger wave in the magnetostriction curve as well. But 0.0001 step size is fine for me till now. Here is the code that I am using for now,

clear

%% Given parameters

a = 512;

c = 0.75;

k = 570;

M_s = 1.7*10^6;

l_s = 10*10^-6;

alpha = 0.032;

u_0 = 4*pi*10^-7;

f = 1;

h_max = 8000;

%% Generate time vector

t_start = 0; %start time

t_end = 60; %end time

dt = 0.0001; %time steps

t = t_start:dt:t_end;

%% Initial values

M_irev(1) = 0; %Initial irreversible magnetization

h = h_max*sin(2*pi*f*t); % Sinusoidal magnetic field strength

%% Time stepping loop

for i=2:length(t)

% Compute effective field

h_e = h(i); % Applied field and effective field considered equal

% Compute anhysteric magnetization

M_an = M_s*(coth(h_e/a) - (a/(h_e))); % Langevin function

% Compute derivative of magnetic field with respect to time

dh_e = (h(i) - h(i-1)) / dt;

delta = sign(dh_e);

% Compute irreversible magnetization using Euler method

M_irev(i) = ((M_an - M_irev(i-1)) / ((k * delta) - (alpha * (M_an - M_irev(i-1))))) * dh_e;

% Update magnetization

M(i) = c*M_an + (1-c)*M_irev(i);

end

%% Magnetic flux density

B = u_0 * (h + M);

%% cycle by cycle averaging

dtzc = 1/f;

tzc = t(1):dtzc:t(end);

NN = numel(tzc)-1;

hh = 0;

BB = 0;

MM = 0;

for k = 1:NN

tstart = tzc(k);

tstop = tzc(k+1);

tt = linspace(tstart,tstop,360); %360 => 1 degree resolution

hi = interp1(t,h,tt);

Bi = interp1(t,B,tt);

Mi = interp1(t,M,tt);

hh = hh + hi;

BB = BB + Bi;

MM = MM + Mi;

end

% divide by N to obtain the averaged data

hh = hh/NN;

BB = BB/NN;

MM = MM/NN;

Bs= smoothdata(BB,'loess',20);

hs= smoothdata(hh,'loess',20);

Ms= smoothdata(MM,'loess',20);

%% Compute magnetostriction

l_ms = l_s*((Ms/M_s).^2);

%% Plot results

%Magnetic Field Strength vs Magnetic Flux Density

subplot(3,1,1);

plot(hh, Bs, "k");hold on;

%plot(h, B, hh, BB, hh, Bs); %Validate the estimation

xlabel('Magnetic Field Strength (H)');

ylabel('Magnetic Flux Density (B)');

title('HB Graph');

grid on;

%Magnetostriction As Function of Magnetic Flux Density

subplot(3,1,2);

plot(Bs,l_ms, "g");

xlabel('Magnetic Flux Density (B)');

ylabel('Magnetostriction (l_ms)');

title('Magnetostriction As Function of Magnetic Flux Density');

grid on;

%Magnetostriction As Function of Magnetic Field Strength

subplot(3,1,3);

plot(hs,l_ms, "b");

xlabel('Magnetic Field Strength (H)');

ylabel('Magnetostriction (l_ms)');

title('Magnetostriction As Function of Magnetic Field Strength');

grid on;

Mathieu NOE 2024-5-6

ok, glad you have now a work around !

请先登录，再进行评论。

How to avoid the overshoot in Jiles Atherton hysteresis loop

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

7 个评论
显示 5更早的评论隐藏 5更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How to avoid the overshoot in Jiles Atherton hysteresis loop

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

7 个评论 显示 5更早的评论隐藏 5更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论