【收不到尾款，开源了】永磁体fem模式

 

 

 

 

function main()
    % 示例参数
    parameters.Br = 1.2;
    parameters.mu0 = 4 * pi * 10^-7;
    parameters.G = 2;
    parameters.dk = 0.0;
    parameters.ws = 0.1;
    parameters.N = 100;
    parameters.Qs = 50;
    parameters.ss = 0.05;
    parameters.mucx = 1.2;
    parameters.mucz = 1.0;
    parameters.Rm = 0.1;
    parameters.muc_theta = ones(10, 10);
    parameters.muc_z = ones(10, 10);
    parameters.K_theta = diag([1, 2, 3]);
    parameters.aI = 1.0;
    parameters.bI = 0.5;
    parameters.aII = 1.5;
    parameters.bII = 0.75;
    parameters.aIII = 2.0;
    parameters.bIII = 1.0;
    parameters.z0 = 0.0;
    parameters.z1 = 0.5;
    parameters.z2 = 1.0;
    parameters.z3 = 1.5;
    parameters.Mn = 0.01;
    parameters.i_Jn = 0.02;
    parameters.mn_n = 2.0;
    parameters.Bo = 1.388;
    parameters.Ho = 239.571;
    parameters.v = 19.985;
    parameters.phi = rand(10, 10); % 确保 phi 的维度与 theta 的长度一致

    % 示例角度坐标
    theta = linspace(0, 2 * pi, 10);

    % 开路条件下的磁通密度计算
    [B_theta_AM_open, B_z_AM_open] = magnetic_flux_density_open_circuit(parameters, theta);
    B_theta_FEM_open = rand(1, 100);
    B_z_FEM_open = rand(1, 100);

    % 负载条件下的磁通密度计算
    [B_theta_AM_load, B_z_AM_load] = magnetic_flux_density_load_condition(parameters, theta);
    B_theta_FEM_load = rand(1, 100);
    B_z_FEM_load = rand(1, 100);

    % 计算均方根误差
    rms_error_open = calculate_rms_error(B_z_FEM_open, B_z_AM_open);
    rms_error_load = calculate_rms_error(B_z_FEM_load, B_z_AM_load);

    fprintf('RMS Error (Open Circuit): %.2e T\n', rms_error_open);
    fprintf('RMS Error (Load Condition): %.2e T\n', rms_error_load);

    % 绘制磁通密度对比图
    plot_magnetic_flux_density(theta, B_theta_FEM_open, B_z_FEM_open, B_theta_AM_open, B_z_AM_open, 'Open Circuit');
    plot_magnetic_flux_density(theta, B_theta_FEM_load, B_z_FEM_load, B_theta_AM_load, B_z_AM_load, 'Load Condition');

    % 验证实验数据
    angular_positions = linspace(0, 360, 361);
    measured_cogging_torque = rand(1, 361);
    calculated_cogging_torque = rand(1, 361);
    measured_emf = rand(1, 361);
    calculated_emf = rand(1, 361);
    measured_static_torque = rand(1, 361);
    calculated_static_torque = rand(1, 361);

    % 齿槽转矩实验验证
    cogging_torque_error = cogging_torque_experiment(measured_cogging_torque, calculated_cogging_torque, angular_positions);
    fprintf('Relative Error (Cogging Torque): %.2f%%\n', cogging_torque_error);

    % 无负载反电动势（EMF）实验验证
    emf_rms_error = no_load_emf_experiment(measured_emf, calculated_emf);
    fprintf('RMS Error (No Load EMF): %.2f V\n', emf_rms_error);

    % 静态扭矩实验验证
    static_torque_error = static_torque_experiment(measured_static_torque, calculated_static_torque, angular_positions);
    fprintf('Relative Error (Static Torque): %.2f%%\n', static_torque_error);

    % 绘制对比图
    plot_comparison('Comparison of Cogging Torque', 'Angular Position [deg]', 'Torque [Nm]', angular_positions, measured_cogging_torque, calculated_cogging_torque, 'Measured', 'Calculated');
    plot_comparison('Comparison of No Load EMF', 'Angular Position [deg]', 'EMF [V]', angular_positions, measured_emf, calculated_emf, 'Measured', 'Calculated');
    plot_comparison('Comparison of Static Torque', 'Angular Position [deg]', 'Torque [Nm]', angular_positions, measured_static_torque, calculated_static_torque, 'Measured', 'Calculated');
end

function [Hx, Hz] = magnetic_field_scalar_potential(parameters, x, z)
    theta = linspace(0, 2 * pi, length(x));
    [phiI, phiII, phiIII] = magnetic_scalar_potential_regions(parameters, z, theta);
    phi = phiI; % 简化示例，这里我们假设使用phiI作为磁标势

    Hx = -gradient(phi, 1);
    Hz = -gradient(phi, 2);
end

function [phiI, phiII, phiIII] = magnetic_scalar_potential_regions(parameters, z, theta)
    Rm = parameters.Rm;
    mu0 = parameters.mu0;
    muc_z = parameters.muc_z;
    muc_theta = parameters.muc_theta;

    V = sqrt(diag(inv(muc_z) * parameters.K_theta * muc_theta * parameters.K_theta));

    Pomega = @(alpha, beta, chi) alpha .* exp((beta - chi) ./ alpha);

    aI = parameters.aI;
    bI = parameters.bI;
    aII = parameters.aII;
    bII = parameters.bII;
    aIII = parameters.aIII;
    bIII = parameters.bIII;

    phiI = aI * Pomega(Rm, z, parameters.z1) + bI * Pomega(Rm, parameters.z0, z);
    phiII = aII * Pomega(Rm, z, parameters.z2) + bII * Pomega(Rm, parameters.z1, z);
    phiIII = aIII * Pomega(Rm, z, parameters.z3) + bIII * Pomega(Rm, parameters.z2, z);
end

function [Mz_n, Mtheta_n] = magnetization_components(parameters)
    Br = parameters.Br;
    mu0 = parameters.mu0;
    G = parameters.G;
    dk = parameters.dk;

    p = 2;
    N = 5;

    Mz_n = zeros(1, N);
    Mtheta_n = zeros(1, N);

    for n = 1:N
        Mz_sum = 0;
        Mtheta_sum = 0;
        for g = 1:(2 * G)
            theta_s = (g - 1) * pi / G;
            Mz_sum = Mz_sum + cos(theta_s) * (exp(1i * n * (2 * g - 1) * pi / (2 * G)) - exp(1i * n * (2 * g - 3) * pi / (2 * G)));
            Mtheta_sum = Mtheta_sum + sin(theta_s) * (exp(1i * n * (2 * g - 1) * pi / (2 * G)) - exp(1i * n * (2 * g - 3) * pi / (2 * G)));
        end
        Mz_n(n) = (Br / mu0) * Mz_sum / (2 * pi * 1i * n) * exp(1i * n * dk);
        Mtheta_n(n) = (Br / mu0) * Mtheta_sum / (2 * pi * 1i * n) * exp(1i * n * dk);
    end
end

function Jz_n = current_density_fourier(parameters)
    N = parameters.N;
    ws = parameters.ws;
    Qs = parameters.Qs;
    ss = parameters.ss;

    Jz_n = zeros(1, N);

    for n = 1:N
        sum_term = 0;
        for i = 1:Qs
            Ji1 = 1.0;
            Ji2 = 2.0;
            ai = 0.0;

            term1 = Ji1 * exp(1i * n * ws) * exp(-1i * n * ss / 2);
            term2 = Ji2 * (1 - exp(-1i * n * ws)) * exp(1i * n * ss / 2);

            sum_term = sum_term + term1 + term2;
        end
        Jz_n(n) = sum_term / (2 * pi * 1i * n);
    end
end

function [Bx, Bz] = calculate_flux_density(parameters, Hx, Hz)
    mucx = parameters.mucx;
    mucz = parameters.mucz;

    Bx = mucx * Hx;
    Bz = mucz * Hz;
end

function [Hz_I, Hz_II, Hz_III] = connection_between_regions(parameters, z, theta, phiI, phiII, phiIII)
    Rm = parameters.Rm;
    mu0 = parameters.mu0;

    Hz_I = -mu0 * gradient(phiI, 1) + mu0 * parameters.Mn;
    Hz_II = -mu0 * gradient(phiII, 1);

    Hz_I = Hz_I - parameters.i_Jn;
    Hz_II = Hz_II;

    Hz_II = -mu0 * gradient(phiII, 1);
    Hz_III = -parameters.mn_n * gradient(phiIII, 1);

    Hz_I = Hz_I / mu0;
    Hz_II = Hz_II / mu0;
    Hz_III = Hz_III / parameters.mn_n;
end

function parameters = iterate_permeability(parameters, tolerance)
    if nargin < 2
        tolerance = 1e-6;
    end
    mu_r_Max = 4500;
    parameters.muc_z = mu_r_Max;
    parameters.muc_theta = mu_r_Max;

    diff = tolerance + 1;
    while diff > tolerance
        new_muc_z = parameters.muc_z - (parameters.muc_z - mu_r_Max) / 2;
        new_muc_theta = parameters.muc_theta - (parameters.muc_theta - mu_r_Max) / 2;

        diff = max(max(abs(new_muc_z - parameters.muc_z)), max(abs(new_muc_theta - parameters.muc_theta)));

        parameters.muc_z = new_muc_z;
        parameters.muc_theta = new_muc_theta;
    end
end

function Rave = calculate_average_radius(Ri, Ro, nsl, k)
    Rave = Ri + (Ro - Ri) / (2 * nsl) * (2 * k - 1);
end

function tcp = calculate_slice_width(Ri, Ro, nsl)
    tcp = (Ro - Ri) / nsl;
end

function machine_description = describe_machine()
    machine_description.type = '10-slot/4-pole five-phase AFPM machine';
    machine_description.material = 'Amorphous Magnetic Material (AMM)';
    machine_description.topology = 'Open slot';
    machine_description.dimensions.slots = 10;
    machine_description.dimensions.poles = 4;
    machine_description.dimensions.phases = 5;
    machine_description.dimensions.inner_diameter = 50.0;
    machine_description.dimensions.outer_diameter = 100.0;
    machine_description.dimensions.height_PM = 10.0;
    machine_description.dimensions.airgap_length = 1.0;
end

function layers = describe_layers(Ri, Ro, nsl)
    layers = struct();
    for k = 1:nsl
        Rave = calculate_average_radius(Ri, Ro, nsl, k);
        tcp = calculate_slice_width(Ri, Ro, nsl);
        layers(k).layer = k;
        layers(k).average_radius = Rave;
        layers(k).slice_width = tcp;
    end
end

function mu_r = relative_permeability(B, Bo, Ho, v)
    mu_r = (Bo / (Ho * 4 * pi * 10^-7)) * (1 / (1 + (B / Bo)^(v - 1)));
end

function [B_theta, B_z] = magnetic_flux_density_open_circuit(parameters, theta)
    Rm = parameters.Rm;
    muc_theta = parameters.muc_theta;
    muc_z = parameters.muc_z;

    N = parameters.N;
    N_theta = diag(-N:N);

    phi = parameters.phi;

    % 确保 muc_theta 和 N_theta 的维度与 phi 的维度匹配
    size_theta = length(theta);
    muc_theta = eye(size_theta);
    N_theta = diag(-floor(size_theta / 2):floor(size_theta / 2) - 1);

    B_theta = -1i * (1 / Rm) * muc_theta * N_theta * phi(1:size_theta, 1:size_theta);
    B_z = -muc_z(1:size_theta, 1:size_theta) .* gradient(phi(1:size_theta, 1:size_theta), 1);
end

function [B_theta, B_z] = magnetic_flux_density_load_condition(parameters, theta)
    [B_theta, B_z] = magnetic_flux_density_open_circuit(parameters, theta);
end

function rms_error = calculate_rms_error(B_FEM, B_AM)
    if ~isequal(size(B_FEM), size(B_AM))
        B_FEM = B_FEM(:);
        B_AM = B_AM(:);
    end
    rms_error = sqrt(mean((B_FEM - B_AM).^2));
end

function plot_magnetic_flux_density(theta, B_theta_FEM, B_z_FEM, B_theta_AM, B_z_AM, condition)
    figure;
    subplot(2, 1, 1);
    plot(linspace(0, 2 * pi, length(B_theta_FEM)), real(B_theta_FEM), 'DisplayName', 'FEM Tangential');
    hold on;
    plot(theta, real(B_theta_AM), '--', 'DisplayName', 'AM Tangential');
    title([condition, ' Magnetic Flux Density Comparison']);
    xlabel('Theta [rad]');
    ylabel('B_theta [T]');
    legend;
    grid on;

    subplot(2, 1, 2);
    plot(linspace(0, 2 * pi, length(B_theta_FEM)), real(B_z_FEM), 'DisplayName', 'FEM Normal');
    hold on;
    plot(theta, real(B_z_AM), '--', 'DisplayName', 'AM Normal');
    xlabel('Theta [rad]');
    ylabel('B_z [T]');
    legend;
    grid on;
end

function relative_error = cogging_torque_experiment(measured_torque, calculated_torque, angular_positions)
    relative_error = mean(abs(measured_torque - calculated_torque) ./ measured_torque * 100);
end

function rms_error = no_load_emf_experiment(measured_emf, calculated_emf)
    rms_error = sqrt(mean((measured_emf - calculated_emf).^2));
end

function relative_error = static_torque_experiment(measured_torque, calculated_torque, angular_positions)
    relative_error = mean(abs(measured_torque - calculated_torque) ./ measured_torque * 100);
end

function plot_comparison(title_str, xlabel_str, ylabel_str, angular_positions, measured, calculated, measured_label, calculated_label)
    figure;
    % 确保 angular_positions 和 measured 是相同长度的向量
    if length(angular_positions) ~= length(measured)
        error('angular_positions and measured must be the same length.');
    end
    plot(angular_positions, real(measured), 'DisplayName', measured_label, 'Color', 'blue');
    hold on;
    % 确保 angular_positions 和 calculated 是相同长度的向量
    if length(angular_positions) ~= length(calculated)
        error('angular_positions and calculated must be the same length.');
    end
    plot(angular_positions, real(calculated), '--', 'DisplayName', calculated_label, 'Color', 'red');
    title(title_str);
    xlabel(xlabel_str);
    ylabel(ylabel_str);
    legend;
    grid on;
end

 

代码计算流程和科学理论依据

该代码用于分析和验证轴向磁通永磁电机（AFPM）的磁场分布和性能，包含磁标势、磁化分量、电流密度和磁通密度的计算，以及实验数据的验证和对比。以下是代码的计算流程和依据的科学理论。

 

1. 磁标势模型计算磁场分量 Hx 和 Hz

函数: magnetic_field_scalar_potential

 

输入参数:

 

parameters: 包含模型参数的字典。

x, z: x 和 z 方向的坐标。

计算流程:

 

调用 magnetic_scalar_potential_regions 函数计算不同区域的磁标势 ϕI, ϕII, ϕIII。

假设使用 ϕI 作为磁标势 ϕ。

计算磁场强度分量 Hx 和 Hz，使用负梯度计算方法。

科学理论:

 

基于磁标势的计算，通过求解磁标势的负梯度得到磁场强度分量。

2. 计算不同区域的磁标势 ϕ

函数: magnetic_scalar_potential_regions

 

输入参数:

 

parameters: 包含模型参数的字典。

z, theta: z 方向的坐标和角度坐标。

计算流程:

 

计算变量 V，用于后续计算。

定义函数 Pω 用于简化计算。

计算不同区域的磁标势 ϕI, ϕII, ϕIII。

科学理论:

 

通过求解拉普拉斯方程得到不同区域的磁标势分布。

3. 计算磁体磁化的法向分量 Mz,n 和切向分量 Mθ,n

函数: magnetization_components

 

输入参数:

 

parameters: 包含模型参数的字典。

计算流程:

 

初始化 Mz_n 和 Mθ_n。

遍历谐波数，计算每个谐波的磁化分量。

科学理论:

 

基于傅里叶级数展开，计算磁体的法向和切向磁化分量。

4. 计算电流密度的傅里叶系数 Jz,n

函数: current_density_fourier

 

输入参数:

 

parameters: 包含模型参数的字典。

计算流程:

 

初始化 Jz_n。

遍历谐波数，计算每个谐波的电流密度傅里叶系数。

科学理论:

 

使用傅里叶级数展开，计算电流密度的频域表示。

5. 计算磁通密度 Bx 和 Bz

函数: calculate_flux_density

 

输入参数:

 

parameters: 包含模型参数的字典。

Hx, Hz: 磁场强度分量。

计算流程:

 

计算磁通密度 Bx 和 Bz。

科学理论:

 

根据磁场强度和材料的磁导率计算磁通密度。

6. 计算不同区域之间的连接条件

函数: connection_between_regions

 

输入参数:

 

parameters: 包含模型参数的字典。

z, theta: z 方向的坐标和角度坐标。

ϕI, ϕII, ϕIII: 在不同区域的磁标势。

计算流程:

 

计算各区域的磁场强度分量 Hz。

科学理论:

 

基于磁通连续性和边界条件计算不同区域之间的连接条件。

7. 迭代计算考虑饱和效应的相对磁导率

函数: iterate_permeability

 

输入参数:

 

parameters: 包含模型参数的字典。

tolerance: 收敛容差。

计算流程:

 

初始化相对磁导率。

使用牛顿法迭代更新相对磁导率，直到收敛。

科学理论:

 

考虑磁性材料的饱和效应，迭代更新磁导率。

8. 计算第 k 层的平均半径

函数: calculate_average_radius

 

输入参数:

 

Ri, Ro: 内径和外径。

nsl: 分层数量。

k: 当前层的索引。

计算流程:

 

计算第 k 层的平均半径 Rave。

科学理论:

 

基于几何关系计算层的平均半径。

9. 计算每一层的宽度

函数: calculate_slice_width

 

输入参数:

 

Ri, Ro: 内径和外径。

nsl: 分层数量。

计算流程:

 

计算每一层的宽度 tcp。

科学理论:

 

基于几何关系计算层的宽度。

10. 描述所研究的 AFPM 机器

函数: describe_machine

 

输出:

返回包含机器主要参数的字典。

11. 描述 AFPM 机器的各层信息

函数: describe_layers

 

输入参数:

 

Ri, Ro: 内径和外径。

nsl: 分层数量。

输出:

 

返回包含各层信息的字典列表。

12. 计算相对磁导率

函数: relative_permeability

 

输入参数:

 

B, Bo, Ho, v: 常数和磁通密度。

计算流程:

 

计算相对磁导率 μr。

科学理论:

 

基于材料特性曲线计算相对磁导率。

13. 计算开路条件下的磁通密度

函数: magnetic_flux_density_open_circuit

 

输入参数:

 

parameters: 包含模型参数的字典。

theta: 角度坐标。

计算流程:

 

计算开路条件下的磁通密度 B_theta 和 B_z。

科学理论:

 

使用磁标势和材料特性计算开路条件下的磁通密度。

14. 计算负载条件下的磁通密度

函数: magnetic_flux_density_load_condition

 

输入参数:

 

parameters: 包含模型参数的字典。

theta: 角度坐标。

计算流程:

 

调用 magnetic_flux_density_open_circuit 函数计算负载条件下的磁通密度。

科学理论:

 

使用与开路条件相同的方法计算负载条件下的磁通密度。

15. 计算均方根误差（RMS error）

函数: calculate_rms_error

 

输入参数:

 

B_FEM, B_AM: 有限元模型和解析模型计算的磁通密度。

计算流程:

 

展平数组以确保形状匹配。

计算均方根误差。

16. 绘制磁通密度的对比图

函数: plot_magnetic_flux_density

 

输入参数:

 

theta: 角度坐标。

B_theta_FEM, B_z_FEM, B_theta_AM, B_z_AM: 有限元模型和解析模型计算的磁通密度。

condition: 条件描述（开路或负载）。

计算流程:

 

绘制 B_theta 和 B_z 的对比图。

17. 验证齿槽转矩的实验数据

函数: cogging_torque_experiment

 

输入参数:

 

measured_torque, calculated_torque: 实验测量的齿槽转矩和解析模型计算的齿槽转矩。

angular_positions: 角度位置数组。

计算流程:

 

计算相对误差百分比。

18. 验证无负载反电动势（EMF）的实验数据

函数: no_load_emf_experiment

 

输入参数:

 

measured_emf, calculated_emf: 实验测量的 EMF 和解析模型计算的 EMF。

计算流程:

 

计算均方根误差。

19. 验证静态扭矩的实验数据

函数: static_torque_experiment

 

输入参数:

 

measured_torque, calculated_torque: 实验测量的静态扭矩和解析模型计算的静态扭矩。

angular_positions: 角度位置数组。

计算流程:

 

计算相对误差百分比。

20. 绘制对比图

函数: plot_comparison

 

输入参数:

 

title, xlabel, ylabel: 图表标题和轴标签。

angular_positions: 角度位置数组。

measured, calculated: 实验测量值和计算值。

measured_label, calculated_label: 实验测量值和计算值的标签。

计算流程:

 

绘制实验测量值和计算值的对比图。

示例代码

代码示例展示了如何使用上述函数计算和验证 AFPM 机器的磁通密度和相关参数，并绘制对比图表。通过这些步骤，可以验证解析模型的有效性并与实验数据进行对比。

 

以下是代码中第1, 2, 3, 4, 5, 6, 7和第12, 13, 14函数的计算公式和详细依据：

所有公式为LATEX代码，请使用相关插件转换。

### 1. 计算磁标势模型中的磁场分量 Hx 和 Hz

**函数: `magnetic_field_scalar_potential`**

 

- **计算公式:**

  \[

  H_x = -\frac{\partial \varphi}{\partial x}

  \]

  \[

  H_z = -\frac{\partial \varphi}{\partial z}

  \]

 

 

- **详细依据:**

  根据磁标势理论，磁场强度 \(\mathbf{H}\) 与磁标势 \(\varphi\) 之间的关系为：

  \[

  \mathbf{H} = -\nabla \varphi

  \]

  在二维直角坐标系中，\(\mathbf{H}\) 的分量可以表示为 \(\varphi\) 对 x 和 z 的梯度。

 

### 2. 计算不同区域的磁标势 ϕ

**函数: `magnetic_scalar_potential_regions`**

 

- **计算公式:**

  \[

  \varphi_I = a_I P_\omega(R_m, z, z_1) + b_I P_\omega(R_m, z_0, z)

  \]

  \[

  \varphi_{II} = a_{II} P_\omega(R_m, z, z_2) + b_{II} P_\omega(R_m, z_1, z)

  \]

  \[

  \varphi_{III} = a_{III} P_\omega(R_m, z, z_3) + b_{III} P_\omega(R_m, z_2, z)

  \]

 

  其中：

  \[

  P_\omega(\alpha, \beta, \chi) = \alpha e^{(\beta - \chi) / \alpha}

  \]

 

- **详细依据:**

  使用变量分离法求解拉普拉斯方程得到不同区域的磁标势分布。假设磁标势在每个区域都可以用指数函数表示，并且通过边界条件连接不同区域的解。

 

### 3. 计算磁体磁化的法向分量 Mz,n 和切向分量 Mθ,n

**函数: `magnetization_components`**

 

- **计算公式:**

  \[

  M_{z,n} = \frac{Br}{2 \pi j n \mu_0} \sum_{g=1}^{2G} \cos\left(\theta_s\right) \left(e^{j n (2g-1) \pi / (2G)} - e^{j n (2g-3) \pi / (2G)}\right) e^{j n \delta_k}

  \]

  \[

  M_{\theta,n} = \frac{Br}{2 \pi j n \mu_0} \sum_{g=1}^{2G} \sin\left(\theta_s\right) \left(e^{j n (2g-1) \pi / (2G)} - e^{j n (2g-3) \pi / (2G)}\right) e^{j n \delta_k}

  \]

 

- **详细依据:**

  使用傅里叶级数展开表示磁化矢量的法向和切向分量。通过求和得到每个谐波分量的表达式，涉及到磁体的剩磁强度 \(Br\)、磁导率 \(\mu_0\) 和块数 \(G\)。

 

### 4. 计算电流密度的傅里叶系数 Jz,n

**函数: `current_density_fourier`**

 

- **计算公式:**

  \[

  J_{z,n} = \frac{1}{2 \pi j n} \sum_{i=1}^{Qs} \left[ J_{i,1} e^{j n w_s} e^{-j n \theta_{ss} / 2} + J_{i,2} (1 - e^{-j n w_s}) e^{j n \theta_{ss} / 2} \right] e^{j n \alpha_i}

  \]

 

- **详细依据:**

  使用傅里叶级数展开电流密度分布。通过对每个槽的电流密度进行求和，得到每个谐波的傅里叶系数。

 

### 5. 计算磁通密度 Bx 和 Bz

**函数: `calculate_flux_density`**

 

- **计算公式:**

  \[

  B_x = \mu_{cx} H_x

  \]

  \[

  B_z = \mu_{cz} H_z

  \]

 

- **详细依据:**

  根据磁场强度和材料的磁导率计算磁通密度。使用磁导率和磁场强度的乘积得到磁通密度。

 

### 6. 计算不同区域之间的连接条件

**函数: `connection_between_regions`**

 

- **计算公式:**

  \[

  H_{\theta,I}|_{z=z_0} = 0

  \]

  \[

  H_{\theta,III}|_{z=z_3} = 0

  \]

  \[

  \mu_0 H_{z,I}(z_1, \theta) - \mu_0 H_{z,II}(z_1, \theta) = \mu_0 M_n

  \]

  \[

  \mu_0 H_{z,II}(z_2, \theta) - \mu_{n,n} H_{z,III}(z_2, \theta) = 0

  \]

 

- **详细依据:**

  基于磁通连续性和边界条件计算不同区域之间的连接条件。通过边界上的磁场强度和磁化矢量的关系，得到区域之间的连接条件。

 

### 7. 迭代计算考虑饱和效应的相对磁导率

**函数: `iterate_permeability`**

 

- **计算公式:**

  迭代更新相对磁导率 \(\mu_c\)：

  \[

  \mu_{c,z}^{\text{new}} = \mu_{c,z} - \frac{\mu_{c,z} - \mu_{r,\text{Max}}}{2}

  \]

  \[

  \mu_{c,\theta}^{\text{new}} = \mu_{c,\theta} - \frac{\mu_{c,\theta} - \mu_{r,\text{Max}}}{2}

  \]

 

- **详细依据:**

  使用牛顿法迭代更新相对磁导率，考虑材料的饱和效应。通过不断调整磁导率，使得计算结果满足一定的收敛条件。

 

### 12. 计算相对磁导率

**函数: `relative_permeability`**

 

- **计算公式:**

  \[

  \mu_r = \frac{B_o}{H_o \mu_0} \frac{1}{1 + (B / B_o)^{v-1}}

  \]

 

- **详细依据:**

  根据材料的磁化曲线，通过拟合参数 \(B_o\)、\(H_o\) 和 \(v\)，计算不同磁通密度下的相对磁导率。

 

### 13. 计算开路条件下的磁通密度

**函数: `magnetic_flux_density_open_circuit`**

 

- **计算公式:**

  \[

  B_\theta = -j \frac{1}{R_m} \mu_{c,\theta} N_\theta \varphi

  \]

  \[

  B_z = -\mu_{c,z} \frac{\partial \varphi}{\partial z}

  \]

 

- **详细依据:**

  使用磁标势和材料特性计算开路条件下的磁通密度。通过磁标势的梯度和材料的磁导率，得到磁通密度的切向和法向分量。

 

### 14. 计算负载条件下的磁通密度

**函数: `magnetic_flux_density_load_condition`**

 

- **计算公式:**

  调用 `magnetic_flux_density_open_circuit` 函数，计算负载条件下的磁通密度。

 

- **详细依据:**

  负载条件下的磁通密度计算方法与开路条件相同，因此直接调用 `magnetic_flux_density_open_circuit` 函数即可。

 

在 MATLAB 脚本中，一些细节函数被简化或未包含，因为它们的逻辑已经包含在其他部分或者在 MATLAB 中不需要独立实现：

 

Python 中的 magnetic_field_scalar_potential

python

复制代码

def magnetic_field_scalar_potential(parameters, x, z):

    ...

    ϕI, ϕII, ϕIII = magnetic_scalar_potential_regions(parameters, z, theta)

    ϕ = ϕI # 简化示例，这里我们假设使用ϕI作为磁标势

    Hx = -np.gradient(ϕ, axis=0)

    Hz = -np.gradient(ϕ, axis=1)

    return Hx, Hz

这个函数在 MATLAB 中被简化为直接计算磁场强度。

 

MATLAB 中

matlab

复制代码

function [Hx, Hz] = magnetic_field_scalar_potential(parameters, x, z)

    ...

    [ϕI, ϕII, ϕIII] = magnetic_scalar_potential_regions(parameters, z, theta);

    ϕ = ϕI; % 简化示例，这里我们假设使用ϕI作为磁标势

    Hx = -gradient(ϕ, 1);

    Hz = -gradient(ϕ, 2);

end

 

磁标势计算的区域合并

在 Python 中，计算磁标势时分为了多个区域（ϕI, ϕII, ϕIII），而在 MATLAB 中简化为了直接计算。

 

Python 中的 magnetic_scalar_potential_regions

python

复制代码

def magnetic_scalar_potential_regions(parameters, z, theta):

    ...

    ϕI = aI * Pω(Rm, z, parameters['z1']) + bI * Pω(Rm, parameters['z0'], z)

    ϕII = aII * Pω(Rm, z, parameters['z2']) + bII * Pω(Rm, parameters['z1'], z)

    ϕIII = aIII * Pω(Rm, z, parameters['z3']) + bIII * Pω(Rm, parameters['z2'], z)

    return ϕI, ϕII, ϕIII

MATLAB 中

matlab

复制代码

function [ϕI, ϕII, ϕIII] = magnetic_scalar_potential_regions(parameters, z, theta)

    ...

    ϕI = aI * Pω(Rm, z, z1) + bI * Pω(Rm, z0, z);

    ϕII = aII * Pω(Rm, z, z2) + bII * Pω(Rm, z1, z);

    ϕIII = aIII * Pω(Rm, z, z3) + bIII * Pω(Rm, z2, z);

end