基于 Preisach 模型的磁滞回线计算程序（MATLAB）

Preisach 模型是描述磁性材料磁滞行为最经典的标量模型之一。它将材料视为大量基本磁滞算子（rectangular hysteresis operators）的叠加，每个算子具有独立的矫顽力 ( $\alpha$ ) 和相互作用场 ( $\beta$ )。通过合理选择 Preisach 密度函数，可以模拟各种形状的磁滞回线。

MATLAB 实现，包含：

• Preisach 平面的离散化
• 密度函数的定义（二维高斯分布）
• 磁化强度的递推计算（状态更新）
• 正弦激励下的磁滞回线绘制

1. Preisach 模型基本原理

每个基本算子 ( $\gamma_{\alpha\beta}[u(t)] $) 的输出为：

$${\gamma }_{\alpha \beta }[u(t)]=\{\begin{array}{ll}+1& u(t)>\alpha \\ -1& u(t)<\beta \\ \text{保持前一状态}& \beta \le u(t)\le \alpha \end{array}$$

其中 $\alpha \ge \beta$，$u(t) $为外加磁场强度。总磁化强度为所有算子的加权平均：

[
$$M(t)={\iint }_{\alpha \ge \beta }\mu (\alpha ,\beta ){\gamma }_{\alpha \beta }[u(t)]d\alpha d\beta$$
]

在实际计算中，我们将 Preisach 平面 ( $$(\alpha ,\beta )$$) 离散为 ( $$N\times N$$ ) 网格，每个网格点赋予权重 ( $$w_{ij}=\mu ({\alpha }_i,{\beta }_j)\Delta \alpha \Delta \beta$$ )。

2. MATLAB 完整程序

%% Preisach 模型磁滞回线仿真
clear; clc; close all;

%% 参数设置
% Preisach 平面范围（α, β 的取值范围）
alpha_min = -1.5; alpha_max = 1.5;   % α 范围
beta_min  = -1.5; beta_max  = 1.5;   % β 范围

% 离散网格点数（越大越精确，但计算量增大）
N = 80;                              

% 激励磁场参数
H_max = 1.2;          % 最大磁场强度
f = 1;                % 频率（Hz）
T = 1/f;              % 周期
n_periods = 2;        % 仿真周期数（至少2个周期消除瞬态）
dt = T/200;           % 时间步长
t = 0:dt:n_periods*T; % 时间向量
H = H_max * sin(2*pi*f*t);   % 正弦激励

%% 构建 Preisach 平面网格
alpha_vec = linspace(alpha_min, alpha_max, N);
beta_vec  = linspace(beta_min,  beta_max,  N);
dalpha = alpha_vec(2) - alpha_vec(1);
dbeta  = beta_vec(2) - beta_vec(1);

% 生成网格坐标
[Alpha, Beta] = meshgrid(alpha_vec, beta_vec);

% 只保留 α >= β 的区域（三角形区域）
mask = (Alpha >= Beta);

%% 定义 Preisach 密度函数 μ(α,β) —— 二维高斯分布
% 中心位于 (α0, β0)，标准差 σα, σβ
alpha0 = 0.2;  beta0 = -0.2;
sigma_alpha = 0.4; sigma_beta = 0.4;

mu = exp( -((Alpha - alpha0).^2/(2*sigma_alpha^2) + (Beta - beta0).^2/(2*sigma_beta^2)) );
mu(~mask) = 0;   % 只保留 α>=β 区域

% 归一化，使得总权重和为 1（即饱和磁化强度 Ms = 1）
total_weight = sum(mu, 'all') * dalpha * dbeta;
mu = mu / total_weight;

% 每个网格点的权重 w_ij = μ(α_i, β_j) * dα * dβ
W = mu * dalpha * dbeta;

%% 初始化算子状态
% 状态矩阵 Q：Q(i,j) 表示第 (i,j) 个算子的当前输出（+1 或 -1）
% 初始状态：所有算子处于 -1（对应退磁状态）
Q = -ones(N, N);   % 注意：只有 mask 区域有效

% 存储磁化强度
M = zeros(size(t));

%% 主循环：逐时刻更新算子状态并计算 M
for k = 1:length(t)
    u = H(k);   % 当前磁场
    
    % 更新所有算子的状态
    % 规则：
    %   if u > Alpha(i,j) → Q(i,j) = +1
    %   if u < Beta(i,j)  → Q(i,j) = -1
    %   否则保持不变
    Q(u > Alpha) = 1;      % 超过上行阈值 → 正向饱和
    Q(u < Beta)  = -1;     % 低于下行阈值 → 负向饱和
    
    % 计算总磁化强度（加权平均）
    M(k) = sum(W .* Q, 'all');
end

%% 绘图：磁滞回线
figure('Color','w','Position',[100 100 900 400]);

% 时域波形
subplot(1,2,1);
plot(t, H, 'b', 'LineWidth',1.5); hold on;
plot(t, M, 'r', 'LineWidth',1.5);
xlabel('Time (s)'); ylabel('H or M');
legend('H(t)','M(t)');
title('激励与响应');
grid on;

% 磁滞回线（取最后一个完整周期）
subplot(1,2,2);
% 提取最后一个周期的数据
idx_start = round(length(t)/2) + 1;  % 后半段
plot(H(idx_start:end), M(idx_start:end), 'k-', 'LineWidth',2);
xlabel('Magnetic Field H'); ylabel('Magnetization M');
title('Preisach 模型磁滞回线');
grid on; axis equal;

%% 可选：显示 Preisach 密度函数
figure('Color','w');
surf(Alpha, Beta, mu, 'EdgeColor','none');
xlabel('\alpha'); ylabel('\beta'); zlabel('\mu(\alpha,\beta)');
title('Preisach 密度函数（仅 \alpha≥\beta 区域）');
colorbar;
view(45,30);

3. 运行结果示例

执行上述代码后，你将看到：

• 左图：正弦磁场 ( H(t) ) 和磁化强度 ( M(t) ) 的时域波形，可见 ( M ) 滞后于 ( H )。
• 右图：典型的磁滞回线（蝴蝶形），饱和、剩磁、矫顽力清晰可见。
• 第二张图：Preisach 密度函数的三维曲面，反映不同 ( $(\alpha ,\beta )$ ) 处算子的权重分布。

参考代码 基于prsicach模型的磁滞回线计算程序 www.youwenfan.com/contentcsv/81342.html

4. 参数调节指南

常见密度函数调整技巧：

• 窄高斯分布 → 陡峭的方形回线
• 宽高斯分布 → 圆润的回线
• 双峰分布 → 可模拟具有多个矫顽力成分的材料

5. 扩展到任意激励波形

只需修改 H 向量即可，例如三角波、阶梯波等：

% 三角波激励
H_tri = H_max * sawtooth(2*pi*f*t, 0.5);

6. 性能优化建议

• 使用 parfor 并行更新算子状态（如果网格很大）。
• 预计算 Alpha 和 Beta 矩阵，避免循环内重复生成。
• 对于实时应用，可采用快速 Preisach 模型（如 Everett 函数法）。

7. 参考文献

• Mayergoyz, I. D. (2003). Mathematical Models of Hysteresis. Springer.
• Bertotti, G. (1998). Hysteresis in Magnetism. Academic Press.