CN111709150B - 任意形状磁体磁场空间分布仿真模拟方法 - Google Patents

Abstract

任意形状磁体磁场空间分布仿真模拟方法，涉及计算机仿真技术。本发明包括下述步骤：步骤一、在直角坐标系下，将假想磁荷模型提出的位函数利用拉普拉斯方程进行有限差分，获得三维空间上磁位的迭代公式；步骤二、分别按笛卡尔直角坐标系和Yee网格模型将系统分为两套网格；第一套网格中，采用笛卡尔直角坐标系，根据磁体和添加材料的尺寸进行网格剖分，并建模：步骤三、第二套网格采用Yee网格模型，将步骤二求得的静态磁感应强度B_s利用线性插值的方法，插值到Yee网格模型节点处；步骤四、将Yee网格节点上的电磁场信息线性插值到宏粒子在笛卡尔坐标系中的实际位置处；循环执行步骤三、四，直到模拟结束。本发明的模拟更加符合实际。

Description

任意形状磁体磁场空间分布仿真模拟方法

技术领域

本发明涉及计算机仿真技术。

背景技术

粒子模拟(PIC)随着多年来的不断优化和改进已经渗入到生物医学、磁约束聚变、等离子体物理、微波传输及气体放电等各个先进的科学领域，应运而生出许多简单、高效、计算精准的PIC算法，但仍有边界条件、算法稳定性和数值色散特性等重要问题没有完全解决。因此，对于粒子模拟的研究，大批学者提出各类粒子模型，试图解决关键科学问题。

粒子模拟主要针对的是粒子的运动受力及电荷、电流密度更新等模块，而电磁场更新模块一般根据实际的需求采用不同的算法，例如利用泊松方程求解静电型问题，利用时域有限差分方法(FDTD)求解电磁型问题，而某些准中性或电荷密度极小的问题则不考虑带电粒子对整个电磁场的作用等等。带电粒子的运动和受力除了与电流、电荷激发的库伦电场有关，也与磁场相关，其运动受力服从牛顿-洛伦兹运动方程，粒子的推进基本是采用蛙跳法(LeapFrogMethod)来推进。

基于粒子模拟在各个领域中的丰富应用，将计算任意形状永磁体的空间磁场与粒子模拟相结合，不仅可以使设置空间磁场更加方便快捷，还可以仿真在空间中添加其他磁材料后，永磁体外部磁场的场型变化，以便得到更精确的三维磁场分布，从而反映更真实的粒子受力运动过程。

永磁体属于铁磁性材料，在仿真永磁体的磁场时，需要考虑如何设置永磁体参数才能使仿真符合实际。比较通用方法是：利用毕奥-萨伐尔定律，直接设置磁体的矫顽力、剩磁、磁导率等参数，通过有限积分的方法计算磁场。此种方法仅能实现某些磁体在真空或空气中的模拟，具有一定局限性。而当在空间中添加如顺磁物质、抗磁物质或其他磁材料时，磁场在空间中的场型将或大或小地发生变化，这时此种方法就不再适用。另外，有的方法直接用经验公式求解磁场，即根据经验对某一类磁体直接进行解析计算。此种方法仅能针对特定某一类磁体进行较准确的仿真，即一种经验公式对应某一类磁体材料，当需改变永磁体形状和材料属性时，经验公式也就不再适用，因此，此种方法局限性更大。还有的方法根据实验所得磁场数据将其拟合成函数，通过设置磁场函数来实现。这种方法的局限性在于：对于一些磁场数据，有时很难拟合出合适的函数，又或者拟合出的函数过于复杂和存在较大的误差。总之，上述方法存在以下几个方面的不足：

(1)有限积分方法相较于有限差分方法来计算磁场，前者无法仿真出添加其他材料对磁场的场型所造成的影响，不如后者简单灵活；

(2)经验公式局限性大，仅仅符合“一对一”(即一个经验公式对应某一种磁体)磁体选择上；

(3)利用函数来设置磁场，计算费时且精度不够，况且难以根据不规则磁场的场型来设定相应函数。

发明内容

本发明所要解决的技术问题是，结合粒子模拟，提供一种对任意形状磁体磁场空间分布进行仿真模拟的方法，以精确计算磁场和粒子的运动受力，减小目前粒子模拟的局限性。

本发明解决所述技术问题采用的技术方案是，任意形状磁体磁场空间分布仿真模拟方法，其特征在于，包括下述步骤：

步骤一、在直角坐标系下，将假想磁荷模型提出的位函数利用拉普拉斯方程进行有限差分，获得三维空间上磁位的迭代公式；

步骤二、分别按笛卡尔直角坐标系和Yee网格模型将系统分为两套网格；

第一套网格中，采用笛卡尔直角坐标系，根据磁体和添加材料的尺寸进行网格剖分，并建模：

对于任意形状磁体模型，其网格面心处赋有假想面磁荷密度σ，网格节点处赋有磁势

网格体心处赋有永磁体和其他添加材料的相对磁导率μ_r，静态磁感应强度B_s放置于网格节点处；设置误差精度ε和边界条件，利用磁位迭代公式进行迭代，当满足预设误差精度时迭代停止，得到模拟空间范围内所有的磁位分布；运用磁势求偏导数公式计算静态磁场强度H_s，静态磁场强度H_s乘以其所在空间处权重后的相对磁导率求出此网格节点处的静态磁感应强度B_s；

步骤三、第二套网格采用Yee网格模型，将步骤二求得的静态磁感应强度B_s利用线性插值的方法，插值到Yee网格模型节点处，对于电场平均权重到Yee网格线两端点上，对于磁场则平均权重到Yee网格面四个角点上，并对位于Yee网格节点处的场进行叠加；

步骤四、将Yee网格节点上的电磁场信息线性插值到宏粒子在笛卡尔坐标系中的实际位置处，用牛顿-洛伦兹运动方程计算受力，推动粒子动量的更新，继而推动粒子的坐标位置的更新；

步骤五、按顺序循环执行步骤三、四，直到模拟结束。

所述步骤三中采用下述权重方式，将位于原Yee网格棱心上的电场矢量E及网格面心上的磁场矢量B经过权重，通过线性插值的方法插值到笛卡尔坐标系网格节点上：

对于整网格点(i,j,k)x方向的场分量，

i为笛卡尔坐标系中x方向坐标值，j为笛卡尔坐标系中y方向坐标值，k为笛卡尔坐标系中z方向坐标值。

本发明基于假想磁荷模型，利用有限差分法，根据安培环路定律的微分方程提出的位函数，将拉普拉斯方程进行差分，并将其与粒子模拟方法(PIC)相结合，以使磁体的磁场计算更加方便灵活、粒子的受力更加精准。

在三维建模时，只需要定义假想面磁荷密度的分布及相对磁导率的分布就可以完成永磁体的建模和其相关铁磁材料的添加。然后通过磁场所满足的边界条件和拉普拉斯方程得到标量磁位的差分公式，设置误差精度，反复迭代直到满足精度，从而得到精准的有限空间内所有节点上磁位的值。当磁体或其他材料尺寸过大时，采用非均匀网格来计算和迭代，能有效减少内存占用和计算时间。

本发明方法，可以实现粒子模拟中任意形状永磁体磁场的计算，同时它也可以顾及当在模拟区域中添加其他磁性材料时对磁场分布的影响。

(1)与粒子模拟方法相结合，在仿真带有永磁体材料的有源器件(如负氢离子源、MRI设备等)时，模拟更加符合实际。

(2)等效磁荷面密度的弹性赋值为设计各类永磁体的多种形状和材料属性提供更有效的办法。

附图说明

图1为直角坐标系下永磁体的磁场差分格式示意图。

图2为Yee元胞示意图。

图3为Yee网格中电磁场信息权重到粒子位置示意图。

图4为本发明实现方式的流程图。

图5为本发明方法的模拟结果图。

图6为解析公式模拟结果。

具体实施方式

下面结合附图，对本发明进行详细描述。

本发明提供了一种添加任意形状永磁体空间磁场的方法并与粒子模拟相结合，使得模拟具有不规则永磁体边界的有源器件结果更加符合物理实际。由于永磁体表面赋有假想面磁荷密度，所以形式上类似于电荷密度的处理。

在永磁体产生磁场的模拟区域内，根据安培环路定律的微分形式

引入位函数

与磁场强度建立关系式，

由于磁感应强度的散度恒为0，即

对B＝μH+μ₀M式子左右同时取散度，可得

类似于电位公式的处理，假设

即为假想磁荷密度。在铁磁性物质之外的磁介质中，M可取为0，固ρ_m仅存在于铁磁性物质中，则有

类似于静电场的边界条件处理，推导假想磁荷模型的磁场边界条件为

其中H₁、H₂分别时非铁磁物质一侧和铁磁物质一侧的磁场强度分界面法向分量；

是铁磁物质一侧的磁化强度；

为分界面法向分量；σ是假想面磁荷密度；μ₀是介质空间的磁导率。

将永磁体视作均匀充磁，即ρ_m＝0，于是由

将此拉普拉斯方程离散，并建立直角坐标模型，如图1所示，

位于网格节点上，μ位于网格体心上，σ位于磁体面心上。不同介质分界面上，若分界面法线是沿着x正方向，则磁场满足的边界条件为：

μ₁、μ₂为1、2这两个单元内的相对磁导率，用μ₁，μ₂，...，μ_m来划定材料属性，用来添加其他材料，或者定义永磁体本身参数。同理，另外两个方向处理方法类似。

运用7点差分格式，对8个单元按照拉普拉斯方程进行差分，结合边界条件，可以得到计算标量磁势核心公式：

式中，

Q₀＝(σ₁+σ₂+σ₃+σ₄)/h₃(假设Z方向磁化)

h₁、h₂、h₃分别为x，y，z三个方向上的空间步长，μ₁、μ₂、μ₃、μ₄、μ₅、μ₆、μ₇、μ₈分别为8个单元的磁导率，σ₁、σ₂、σ₃和σ₄为分界面上的假想面磁荷密度，大小为

此式中Q₀的值与₃相关，是因为模型假设沿Z方向磁化。若模型是有多个磁化方向的永磁体时，Q₀的值应该是相应节点四周的面磁荷密度之和除以该σ沿法线方向上的空间步长。Br为剩磁。因此，模型可以采用非均匀网格来计算，可以有效减少总的迭代时间，扩大模拟范围。

满足误差精度后跳出循环，利用求出的标量磁势分布求出静态磁场强度H_s，静态磁场强度H_s乘以其所在空间处权重后的相对磁导率求出网格节点处的静态磁感应强度B_s，然后结合粒子模拟主程序。在Yee网格模型中(如图2)，电场分量是在棱心上的，而磁场分量是在面心上的，粒子则可能存在于网格中的任一位置，所以对于粒子的运动方程式并不能直接使用。因此需要将半网格上的电磁场分量权重到整网格上(即步骤三)，然后再插值到粒子的所在位置，进而用牛顿-洛伦兹运动方程算得受力，从而推动粒子动量的更新，继而推动粒子的坐标位置的更新推动电磁场在空间中的传播，实现粒子模拟的整个过程(即步骤四)。对于步骤四，可利用图3进行阐释说明。首先假定粒子A距离三个坐标面的距离分别为s_x，s_y，s_z。不妨设定：

于是对应的权重因子w有如下计算式：

w_i,j,k＝(1-λ_x)(1-λ_y)(1-λ_z)

w_i+1,j,k＝λ_x(1-λ_y)(1-λ_z)

w_i,j+1,k＝(1-λ_x)λ_y(1-λ_z)

w_i,j,k+1＝(1-λ_x)(1-λ_y)λ_z

w_i+1,j+1,k＝λ_xλ_y(1-λ_z)

w_i+1,j,k+1＝λ_x(1-λ_y)λ_z

w_i,j+1,k+1＝(1-λ_x)λ_yλ_z

w_i+1,j+1,k+1＝λ_xλ_yλ_z

其中，w下标是对应的网格索引。以粒子A在x方向电磁场分量权重后求解为例：

(E_x)_A＝w_i,j,k(E_x)_i,j,k+w_i+1,j,k(E_x)_i+1,j,k+w_i,j+1,k(E_x)_i,j+1,k+w_i,j,k+1(E_x)_i,j,k+1+w_i+1,j+1,k(E_x)_i+1,j+1,k+w_i+1,j,k+1(E_x)_i+1,j,k+1+w_i,j+1,k+1(E_x)_i,j+1,k+1+w_i+1,j+1,k+1(E_x)_i+1,j+1,k+1

(H_x)_A＝w_i,j,k(H_x)_i,j,k+w_i+1,j,k(H_x)_i+1,j,k+w_i,j+1,k(H_x)_i,j+1,k+w_i,j,k+1(H_x)_i,j,k+1+w_i+1,j+1,k(H_x)_i+1,j+1,k+w_i+1,j,k+1(H_x)_i+1,j,k+1+w_i,j+1,k+1(H_x)_i,j+1,k+1+w_i+1,j+1,k+1(H_x)_i+1,j+1,k+1

至此对应于每一个宏粒子的电磁场分量，经8个整网格节点上权重后的电磁场信息插值到粒子在网格中的实际位置，可以明确得出。同理可以推出其他方向的场分量。再由相对论牛顿-洛伦兹力方程即可求出对应宏粒子的位移和动量变化。

作为一个模拟多块矩形永磁体静态磁场分布的例子，其模拟结果将和解析公式的模拟结果进行对比。具体步骤如下：

1、确定模型的基本参数。(即矩形永磁体的尺寸，磁化方向，剩磁大小，磁体个数，相对磁导率等)。

2、网格剖分。具体操作如上述步骤二所述，这里不再赘述，根据等效磁荷模型迭代公式得到多块矩形永磁体的三维磁场分布。

3、根据解析公式，输入磁体的尺寸大小和剩磁大小，根据磁体的磁化方向，然后通过多次磁场叠加的方式来求解多块矩形永磁体的磁场。

4、将获得的解析结果与本发明方法模拟所得结果进行比较，以便验证本模拟方法的正确性，例如，模拟结果对比图，如图5、6所示。

Claims (2)

1.任意形状磁体磁场空间分布仿真模拟方法，其特征在于，包括下述步骤：

步骤一、在直角坐标系下，将假想磁荷模型提出的位函数利用拉普拉斯方程进行有限差分，获得三维空间上磁位的迭代公式；

步骤二、分别按笛卡尔直角坐标系和Yee网格模型将系统分为两套网格；

第一套网格中，采用笛卡尔直角坐标系，根据磁体和添加材料的尺寸进行网格剖分，并建模：

对于任意形状磁体模型，其网格面心处赋有假想面磁荷密度σ，网格节点处赋有磁势

网格体心处赋有永磁体和其他添加材料的相对磁导率μ_r，静态磁感应强度B_s放置于网格节点处；设置误差精度ε和边界条件，利用磁位迭代公式进行迭代，当满足预设误差精度时迭代停止，得到模拟空间范围内所有的磁位分布；运用磁势求偏导数公式计算静态磁场强度H_s，静态磁场强度H_s乘以其所在空间处权重后的相对磁导率求出此网格节点处的静态磁感应强度B_s；

步骤三、第二套网格采用Yee网格模型，将步骤二求得的静态磁感应强度B_s利用线性插值的方法，插值到Yee网格模型节点处，对于电场平均权重到Yee网格线两端点上，对于磁场则平均权重到Yee网格面四个角点上，并对位于Yee网格节点处的场进行叠加；

步骤四、将Yee网格节点上的电磁场信息线性插值到宏粒子在笛卡尔坐标系中的实际位置处，用牛顿-洛伦兹运动方程计算受力，推动粒子动量的更新，继而推动粒子的坐标位置的更新；

循环执行步骤三、四，直到模拟结束。

2.如权利要求1所述的任意形状磁体磁场空间分布仿真模拟方法，其特征在于，所述步骤三中采用下述权重方式，将位于原Yee网格棱心上的电场矢量E及网格面心上的磁场矢量B经过权重，通过线性插值的方法插值到笛卡尔坐标系网格节点上：

对于整网格点(i,j,k)x方向的场分量，

i为笛卡尔坐标系中x方向坐标值，j为笛卡尔坐标系中y方向坐标值，k为笛卡尔坐标系中z方向坐标值。

Images