Clear Filters
Clear Filters

How to avoid the overshoot in Jiles Atherton hysteresis loop

11 views (last 30 days)
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?

Answers (1)

Mathieu NOE
Mathieu NOE on 25 Apr 2024
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 Comments
Md Tanbir Siddik
Md Tanbir Siddik on 6 May 2024
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;

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!