一种微波场混沌搅拌器的数值仿真分析方法

摘要

本发明公开了一种微波场混沌搅拌器的数值仿真分析方法，首先建立微波场混沌搅拌器的有限元模型以及微波场混沌搅拌器中搅拌部件的随机运动函数模型；并定义微波场混沌搅拌器几何模型区域的控制方程和边界约束方程；再利用有限元法求解控制方程得到微波场混沌搅拌器内的电磁场分布。本发明将随机运动函数模型与有限元模型相结合，利用随机函数模型模拟微波场混沌搅拌器中非固定搅拌部件的随机运动，进而得到随搅拌部件变化的移动网格，再利用有限元法求解电磁场的偏微分方程，从而可以获得随时间变化的微波场混沌搅拌器的动态电磁场分布。

说明书

技术领域

本发明属于计算机数据处理技术领域，涉及一种电磁场的数值仿真分析，具体涉及一种应用于微波炉的微波场混沌搅拌器的数值仿真分析方法。

背景技术

随着社会进步和微波能利用技术的发展，微波炉的使用越来越普及。因微波具有良好的穿透性，能直接将能量送入被加热物体内部，实现对物体内部和外部的同时加热，从而极大地缩短了加热时间，提高了加热效率。但由于微波加热的不均匀性会使得食物内部出现过熟和欠熟区域，导致食物无法安全食用，从而制约微波能的有效应用。因此，改善微波加热过程中受热物料温度分布不均的现象显得意义重大。

微波炉加热不均匀的主要原因之一是由于微波炉腔体能量分布的不均匀性；为了改善能量在空间分布的不均匀性，目前的通用做法是添加微波场搅拌器；微波场搅拌器一般由多片共轴安装在在微波炉腔体的顶部的金属扇叶构成，在电机带动，金属扇叶的旋转可以改变微波炉腔体的边界条件，使微波炉内的电磁场分布随着叶片的旋转而跟着旋转或不断的变化，从而可以使微波场的均匀性得到一定的改善；但是微波场搅拌器的位置固定，搅拌具有周期性，这将导致而能量分布也存在周期性的变化，从而使微波场能量均匀性改善存在一定的局限性，改善效果不是很明显。

为克服上述已有技术缺陷，需要设计一种可以实现微波场能量均匀化的非固定式微波场搅拌器；为此，需要先通过数值仿真方法对非固定式微波场搅拌器的搅拌效果进行计算机分析，以得到微波腔体中能量分布情况及相关参数，进而获得微波腔体内电磁分布规律。

目前，人们已经提出了多种方法以实现对电磁环境的仿真计算，例如矩量法、时域有限差分法、有限元法等。其中矩量法是由R.F.Harrington首先应用于电磁问题的求解，其主要思想是将积分方程化为差分方程，或者是将积分方程中积分化为有限求和，从而建立代数方程组，最后输入计算机求解代数方程组；在求解代数方程组的过程中，矩阵规模的大小涉及到占用内存的多少，将会直接影响计算的速度；此外该方法仅适用于计算静止物体的电磁场问题，对于移动物体的实时计算不在适用。时域有限差分法的基本思想是用中心差商代替场量对时间和空间的一阶偏微商，通过在时域的递推模拟波的传播过程，从而得出场分布，该方法适合研究瞬态过程中的电磁波传播问题。有限元法则是将连续的求解域离散为一组单元的个组合体，用在每个单元内假设的近似函数来分片的表示求解域上待求的未知场函数；近似函数通常由未知场函数及其导数在单元各节点的数值插值函数来表达，从而使一个连续的无限自由度问题变成离散的有限自由度问题。

然而，对于非固定式微波场搅拌器，其中的搅拌部分是不固定的，使得微波腔体中不仅含有运动的物体，而且物体的运动具有不确定性，这已经超出上述矩量法和时域有限差分方法的适用范围；此外，由于网格是表征几何的计算区域，物体移动造成网格的变形，使得网格发生反转和畸变导致不能计算，传统有限元法也难以实现非固定式微波场搅拌器微波腔体内电磁场分布的数值模拟。

发明内容

本发明的目的旨在，针对非固定式微波场搅拌器的电磁场分布，缺少有效、可行的数值仿真技术的现状，提供一种微波场混沌搅拌器的数值仿真分析方法，以得到微波腔体中能量分布情况。

为了达到上述目的，本发明提供了一种微波场混沌搅拌器的数值仿真分析方法，包括以下步骤：

(1)在有限元分析软件中绘制或插入微波场混沌搅拌器的几何模型，对建立的微波场混沌搅拌器几何模型进行网格划分，得到微波场混沌搅拌器有限元模型；

(2)定义微波场混沌搅拌器几何模型网格区域的控制方程和边界约束方程；

(3)建立微波场混沌搅拌器中搅拌部件的随机运动函数模型；

(4)求解微波场混沌搅拌器内的动态电磁场分布，包括以下分步骤：

(41)初始化，对步骤(1)得到的微波场混沌搅拌器有限元模型进行初始化，并设定电磁场求解程序运行时间；

(42)函数更新，将当前时刻代入步骤(3)获得的搅拌部件随机运动函数模型，获取搅拌部件位移值，将其作为当前时刻搅拌部件网格节点位移值；并依据搅拌部件网格节点位移值获得搅拌部件之外网格节点位移值，完成一次函数更新；

(43)网格移动，将步骤(42)获取的当前时刻微波场混沌搅拌器网格节点位移值加载到网格上，完成网格移动；

(44)网格畸变判断，判断最小网格质量是否小于设定值，若小于设定值，网格发生畸变，进入步骤(47)；若不小于设定值，网格没有发生畸变，进入步骤(45)；

(45)电磁场分布计算，根据步骤(43)获得的移动之后的网格、边界约束方程以及上一时刻得到的微波场混沌搅拌器的电磁场结果，利用有限元法求解控制方程获得当前时刻微波场混沌搅拌器的电磁场分布；

(46)判断是否达到设定运行时间，若达到设定运行时间，得到微波场混沌搅拌器的动态电磁场分布，任务完成；若没有达到设定时间，在当前时刻基础上增加设定时间间隔作为下一个当前时刻，然后返回步骤(42)；

(47)网格重构，依据搅拌部件当前时刻的位移值，获取当前时刻搅拌部件的坐标值，再根据网格划分原理对搅拌器的几何模型进行剖分，得到重构后的几何模型网格，并以网格畸变时刻为当前时刻返回步骤(42)。

上述微波场混沌搅拌器的数值仿真分析方法，利用随机函数建立微波场混沌搅拌器中搅拌部件的随机运动函数模型，并将该模型与有限元模型相结合，得到不同时刻的移动网格，再利用有限元方法求解获得微波场混沌搅拌器的电磁场分布；为了避免因网格畸变，导致电场强度、磁场强度在一个几何坐标点上出现多个值，造成奇异矩阵，在对建立的微波场混沌搅拌器的有限元模型进行求解得到微波场混沌搅拌器内的电磁场分布过程中，在网格完成移动后，增加了对移动后网格是否发生畸变的判断，若网格发生畸变，就对网格进行重构，通过对网格进行重构，避免因网格畸变出现奇异矩阵而导致的结果不收敛或者结果错误，实现对具有运动物体的电磁场分布的计算。

上述微波场混沌搅拌器的数值仿真分析方法，所述步骤(1)的目的在于构建微波场混沌搅拌器的几何模型，并对其进行网格划分，以便于利用有限元方法求解对电磁场偏微分方程获得微波场混沌搅拌器内的电磁场分布。微波场混沌搅拌器的几何模型可以通过常规绘图软件(如AutoCAD、Solidworks)等，然后插入有限元分析软件中，也可以直接在有限元分析软件中绘制微波场混沌搅拌器的几何模型；所构建的搅拌器几何模型包括微波腔体，位于微波腔体内的搅拌部件以及位于微波腔体一侧、为微波腔体提供激励电磁波的波导；当插入有限元软件的几何模型为三维几何模型时，为了减少计算量，可将三维几何模型简化为二维几何模型，例如可沿微波场混沌搅拌器入射波导轴向取其横截面，得到简化的二维几何模型。当搅拌器几何模型构建好后，便可得知其几何区域包括几何坐标原点、搅拌部件的边长尺寸、微波腔体的边长尺寸、波导尺寸以及波导馈入端口等。

上述微波场混沌搅拌器的数值仿真分析方法，根据有限元求解偏微分方程的流程，在求解前处理中要对几何体进行网格划分，将几何体离散成一个一个小的区域，然后求解其节点上的值；网格的尺寸越小，求解所得到的结果就越精确，但是所耗费的求解资源也越大；对于有限元的网格剖分算法可以采用本领域已经披露的常规算法(参见文献《平面区域有限元三角网格剖分算法研究》)；为了求解电磁波动方程，最大网格单元尺寸应该小于激励电磁波波长的六分之一才能够精确解析电磁波；在优选实施方式中，对于二维几何模型，网格划分为三角形单元或四边形单元。

上述微波场混沌搅拌器的数值仿真分析方法，所述步骤(1)建立微波场混沌搅拌器有限元模型，进一步包括定义微波场混沌搅拌器中搅拌部件材料属性、微波场混沌搅拌器腔体内空间介质属性以及定义空间框架坐标系(固定不动)和材料框架坐标系(随搅拌部件运动而变化)，以便于后期电磁场的计算；若采用的软件已经有自定义的空间框架坐标系(固定不动)和材料框架坐标系(随搅拌部件运动而变化)，可以不需要再重复定义。

上述微波场混沌搅拌器的数值仿真分析方法，为了求解电磁场，需要进一步给出微波场混沌搅拌器几何区域的控制方程和边界约束方程，通过求解控制方程，便可得到电磁场分布，给出的控制方程可以是有关电场的，也可以是有关磁场的；本发明给出的微波场混沌搅拌器所在几何区域的控制方程为：

其中，E为电场强度，k₀自由空间的波矢、μ_r为相对磁导率，ε_r为相对介电常数，ε₀为真空中的介电常数，σ为电导率，ω为角频率，j为虚数

边界约束方程为：

其中，H为磁场强度，E为电场，n为法向矢量，μ₀为真空中的磁导率，μ_r为相对磁导率，ε_r为相对介电常数，ε₀为真空中的介电常数，σ为电导率，ω为角频率。

上述微波场混沌搅拌器的数值仿真分析方法，通过MATLAB数学分析软件的随机函数定义微波场搅拌器中搅拌部件的随机运动。本发明采用ALE(Arbitrary Lag range-Euler)方法对搅拌器中搅拌部件的随机运动进行描述，当搅拌部件为金属搅拌棒时，金属搅拌棒四条边在空间框架坐标系中x轴方向上的运动方程为：

x＝X·cos(2πt)-Y·sin(2πt-X)-int(t)，金属搅拌棒四条边在空间框架坐标系中y轴方向上的运动方程为：y＝Y·cos(2πt)+X·sin(2πt-Y)+int(t)；

当搅拌部件由金属搅拌棒以及设置于金属搅拌棒一端的金属柔性丝带构成时，金属搅拌棒三条边在空间框架坐标系中x轴方向上的运动方程为：

x＝X·cos(2πt)-Y·sin(2πt-X)-int(t)，金属搅拌棒三条边在空间框架坐标系中y轴方向上的运动方程为：y＝Y·cos(2πt)+X·sin(2πt-Y)+int(t)；金属柔性丝带三条边在空间框架坐标系中x轴方向上的运动方程为：

x＝int(t)·X·cos(2πt)-Y·sin(2πt-X)+int(t)，金属柔性丝带三条边在空间框架坐标系中y轴方向上的运动方程为：y＝int(t)·Y·cos(2πt)+X·sin(2πt-Y)+int(t)；其中X，Y为搅拌部件在材料框架坐标系下对应的坐标值，int(t)为随着时间变化的随机变量，将a对应的随机数b通过函数b＝randn(siza(a))得到，从而得到数组(a，b)，然后将数组(a，b)的值通过插值方法得到函数int(t)，其中t为时间变量，a在0到1之间取值。

上述微波场混沌搅拌器的数值仿真分析方法，步骤(4)的目的是利用有限元法通过求解电磁场的偏微分方程得到不同时刻的微波场混沌搅拌器内的电磁场分布，即动态电磁场分布。为了在时间统一尺度下计算搅拌微波场的过程，可以采用时域-有限元方法(TD-FEM)计算整个过程。

步骤(41)中，在计算动态电磁场分布之前,需要进行初始化，得到有关微波场混沌搅拌器有限元模型的初始化参数,并设定电磁场求解程序运行时间；设定电磁场求解程序运行时间，并在该时间段内设定一个时间间隔，当按照步骤(42)至步骤(45)完成一次电磁场计算时，对计算结果进行输出保存，并对运行时间进行判断【即步骤(46)】，若没有达到时间段的设定时间，在当前时刻基础上叠加一个时间间隔，再返回至步骤(42)继续计算下一时刻微波场搅拌器腔体内的电磁场分布；初始化参数包括网格节点上对应的电场或磁场分量和网格移动的位移值,进行初始化，就是对上述参数进行赋值,从而进行下一步计算；为了简化计算，可以对初始化时刻的网格节点上对应的电场或磁场的各个分量以及网格移动位移值赋予零值。

步骤(42)中，由于本发明针对的微波搅拌器中的搅拌部件是非固定式的，每一时刻的位置都不相同，为此，本发明利用随机运动函数模型来模拟搅拌部件的运动，将当前时刻代入搅拌部件的随机函数中，采用本领域常规计算方法(参见《成型充填过程中非等温非牛顿粘性流动的ALE有限元与无网格自适应耦合模拟》)得到当前时刻搅拌部件位移值；由于搅拌部件的运动会带动几何模型中所有网格节点运动，从而使所有网格节点产生位移，因此这里将搅拌部件的位移值作为该部件上网格节点的位移值，并且由该部件的网格主动变化引起了周围空间网格的被动变化，根据本领域已经披露的常规方法(《基于调和函数的移动网格有限元方法》)可以得到搅拌部件之外的其它网格节点被动变形后的网格节点位移值；这里搅拌部件位移值为搅拌部件当前时刻位置相对于初始化时刻位置的位移值或者相对于前一时刻位置的位移值；这里的当前时刻为完成前一次电磁场计算的下一时刻(由前一时刻加上时间间隔得到)，或者为发生网格畸变的时刻。

步骤(43)中，获得的是网格移动后的网格节点坐标值，当步骤(2)搅拌部件位移值为搅拌部件当前时刻相对于初始化时刻的位移值时，是将网格节点位移值加载到初始化时刻网格节点坐标值上，得到网格移动后的网格节点坐标值；当步骤(2)搅拌部件位移值为搅拌部件当前时刻相对于前一时刻的位移值时，是将网格节点位移值加载到前一时刻网格节点坐标值上，得到网格移动后的网格节点坐标值。

步骤(44)中，由于在网格瞬态运动过程中，网格的变化会导致网格发生反转、重叠等畸变问题，从而导致物理量(例如电场强度、磁场强度等)在一个几何坐标点上出现多个值，造成奇异矩阵，因此需要对移动后的网格进行判断，可以根据最小网格质量(可以以网格单元面积、长宽比等作为网格质量)来判断；根据经验数据，给定一个设定值，将得到的最小网格质量与设定值比较，若最小网格质量若不小于设定值，网格没有发生畸变，可以进一步计算电磁场分布；若最小网格质量小于设定值，网格发生畸变，此时，需要停止计算，对网格进行重构。

步骤(45)中，根据当前时刻的网格分布情况、边界约束方程以及上一时刻得到的微波场混沌搅拌器的电磁场分布(例如网格节点上的电场分量或磁场分量，激励电磁波频率、功率等)，利用有限元法(例如《时域有限元电磁计算方法的研究》)求解前面给出的控制方程便可获得搅拌器微波腔体内的电磁场分布。

步骤(47)中，网格重构是对网格进行重新划分，为了避免畸变后网格的影响，网格重构与发生畸变前的网格无关，是在当前时刻搅拌部件位置基础上，根据网格划分原理《平面区域有限元三角网格剖分算法研究》进行剖分的；当前时刻搅拌部件位置可以根据步骤(42)中获得的当前时刻搅拌部件的位移值；当搅拌部件位移值为搅拌部件当前时刻相对于初始化时刻的位移值时，当前时刻搅拌部件坐标值是将搅拌部件位移值加载到初始化时刻搅拌部件坐标值上得到；当搅拌部件位移值为搅拌部件当前时刻相对于前一时刻的位移值时，当前时刻搅拌部件坐标值是将搅拌部件位移值加载到前一时刻搅拌部件坐标值上得到；网格重构后，需要以网格发生畸变的时刻为当前时刻，返回步骤(42)继续求解微波场混沌搅拌器的电磁场分布。

传统电磁场的数值仿真方法，均是针对具有固定结构的微波腔体研究提出的，而对于含有非固定式搅拌部件的微波场搅拌器，采用传统电磁场的数值仿真方法，由于网格畸变等原因，程序容易被中断，难以得到微波场搅拌器的电磁场数值仿真结果。

与现有技术相比，本发明具有以下有益效果：

1、本发明将随机运动函数模型与有限元模型相结合，利用随机函数模型模拟微波场混沌搅拌器中非固定搅拌部件的随机运动，进而得到随搅拌部件变化的移动网格，再利用有限元法求解电磁场的偏微分方程，从而可以获得随时间变化的微波场混沌搅拌器的电磁场分布。

2、本发明针对微波场混沌搅拌器中非固定搅拌部件的随机运动，所引起的网格畸变，采用网格重构方式，避免了因网格畸变出现奇异矩阵而导致的结果不收敛或者结果错误，确保电磁场求解过程的顺利进行。

3、本发明微波场混沌搅拌器的数值仿真结果显示，含有非固定搅拌部件的微波场混沌搅拌器能够使微波场能量分布更均匀，有效避免因能量分布不均匀导致的过热效应。

附图说明

图1为本发明实施例1中微波场混沌搅拌器二维几何模型。

图2为本发明实施例2中微波场混沌搅拌器二维几何模型。

图3为本发明实施例1中使用三角形对图2微波场混沌搅拌器二维几何模型剖分的网格图。

图4为本发明实施例1中通过MATLAB软件绘制的随机图形。

图5为本发明求解微波场混沌搅拌器内电磁场分布的流程示意图。

图6为本发明实施例1(a)与实施例2(b)得到的微波场混沌搅拌器内电场模分布图。

图7为本发明实施例1(a)与实施例2(b)得到的微波场混沌搅拌器内电场模时间平均值对比图。

具体实施方式

以下将结合附图对本发明各实施例的技术方案进行清楚、完整的描述，显然，所描述实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施例，都属于本发明所保护的范围。

实施例1

本实施例以设置有非固定金属搅拌棒的搅拌器为微波场混沌搅拌器，其中非固定金属搅拌棒作为搅拌部件，本实施例需要获得的是搅拌器腔体内的动态电磁场分布，具体过程包括以下步骤：

(1)建立微波场混沌搅拌器有限元模型，具体包括以下分步骤，

1)直接在有限元分析软件中绘制微波场混沌搅拌器的二维几何模型，如图1所示，该几何模型包括微波腔体1，金属搅拌棒2设计在微波腔体1内，矩形波导3位于微波腔体1的一侧；以微波腔体1和金属搅拌棒2正中心为几何坐标原点(0，0)，金属搅拌棒2的边长尺寸为15cm×1cm，微波腔体1的边长尺寸为40cm×40cm，在微波腔体1的左边边上的中间，矩形波导3的尺寸为25cm×10cm，并且以最左边为波导馈入端口。

2)定义金属搅拌棒2以及微波腔体1内介质空气的电磁特性，其参数如表1所示，

表1：金属搅拌棒及空气的电磁特性参数

3)对建立的微波场混沌搅拌器几何模型进行网格划分，本实施例中使用最大单元尺寸为24[mm]，最小单元尺寸为0.0812[mm]的三角形对图1中几何模型进行剖分，得到的网格如图3所示。

(2)定义微波搅拌器所在几何区域的控制方程和边界约束方程，本实施例中几何区域的控制方程为：

其中，E为电场强度(矢量)，k₀自由空间的波矢、μ_r为相对磁导率，ε_r为相对介电常数，ε₀为真空中的介电常数，σ为电导率，ω为角频率，j为虚数

边界约束方程为：

其中，H为磁场强度，E为电场(矢量)，n为法向矢量，μ₀为真空中的磁导率，μ_r为相对磁导率，ε_r为相对介电常数，ε₀为真空中的介电常数，σ为电导率，ω为角频率；

本实施例中端口激励电磁波为TE₁₀波，其频率为2.45GHz、功率为1W。

(3)建立微波搅拌器中金属搅拌棒2的随机运动函数模型；通过MATLAB数学分析软件对金属搅拌棒2的随机运动进行定义，其中随机函数代码为

a＝0:0.01:1,b＝randn(size(a))；通过plot(a,b)绘图得到的随机图形如图4所示；金属搅拌棒四条边在空间框架坐标系x轴方向上的运动方程为：

x＝X·cos(2πt)-Y·sin(2πt-X)-int(t)，金属搅拌棒四条边在空间框架坐标系y轴方向上的运动方程为：y＝Y·cos(2πt)+X·sin(2πt-Y)+int(t)，其中X，Y为材料框架下对应的坐标值；int(t)为随着时间变化的随机变量，将a对应的随机数b通过函数b＝randn(siza(a))得到，从而得到数组(a，b)，然后将数组(a，b)的值通过插值方法得到函数int(t)，其中t为时间变量，a在0到1之间取值。

(4)求解得到微波场混沌搅拌器内的动态电磁场分布，包括以下分步骤：

(41)初始化，由前面网格划分得到的节点坐标，依次对初始化时刻的网格节点上对应的电场分量和网格移动的位移值赋零值，设定激励电磁波为TE₁₀波，其频率为2.45GHz、功率为1W，并设定动态电磁场求解程序运行时间为1s，每隔0.01s输出保存一次结果，即时间间隔为0.01s；

(42)函数更新，将当前时刻代入步骤(3)获得的金属搅拌棒2随机运动函数模型，依据x＝X·cos(2πt)-Y·sin(2πt-X)-int(t)和y＝Y·cos(2πt)+X·sin(2πt-Y)+int(t)确定金属搅拌棒2在空间框架坐标系下相对于初始化时刻位置的位移值，并将其作为当前时刻搅拌部件网格节点位移值，由于搅拌部件网格的移动导致周围空间网格的被动变形，根据文献《基于调和函数的移动网格有限元方法》中的方法可以得到周围网格被动变形后的网格节点位移值，完成一次函数更新；

(43)网格移动，将步骤(42)获取的当前时刻网格节点位移值加载到初始化网格节点坐标上，完成网格移动；

(44)网格畸变判断，本实施例中以三角形的最短边与最长边的比值作为网格质量，加入的判断条件(即设定值)为0.2，判断最小网格质量是否小于0.2，若小于0.2，网格发生畸变，进入步骤(47)；若不小于0.2，网格没有发生畸变，进入步骤(65)；

(45)电磁场分布计算，根据步骤(63)获得的移动之后的网格位置、前面给出的边界约束方程、激励电磁波数据等，利用有限元法求解前面给出的控制方程获得当前时刻微波腔体1内的电场分布；再依据电场与磁场的关系(麦克斯韦方程组)，获得当前时刻微波腔体1内的磁场分布；

(46)判断是否达到设定运行时间，若达到设定运行时间，得到微波腔体1内不同时刻的电磁场分布情况(即动态电磁场分布)，任务完成；若没有达到设定时间，在当前时刻基础上增加设定时间间隔0.01s作为下一个当前时刻，然后返回步骤(42)；

(47)网格重构，由初始化时刻金属搅拌棒2坐标值与当前时刻位移值相加得到当前时刻金属搅拌棒2在空间框架坐标系中的坐标值，再根据网格划分原理对搅拌器的二维几何模型进行剖分，得到重构后的几何模型网格，并以网格畸变时刻为当前时刻返回步骤(42)。

实施例2

本实施例以设置有非固定金属搅拌棒2以及连接在金属搅拌棒2上的金属柔性丝带4的搅拌器为微波场混沌搅拌器，其中非固定金属搅拌棒2以及连接在金属搅拌棒2上的金属柔性丝带4作为搅拌部件，本实施例需要获得的是搅拌器腔体内的动态电磁场分布，具体过程包括以下步骤：

(1)建立微波场混沌搅拌器有限元模型，具体包括以下分步骤，

1)直接在有限元分析软件中绘制微波场混沌搅拌器的二维几何模型，如图2所示，该几何模型包括微波腔体1，金属搅拌棒2以及连接在金属搅拌棒2上的金属柔性丝带4设计在微波腔体1内，矩形波导3位于微波腔体1的一侧，其中金属搅拌棒2和金属柔性丝带4构成搅拌部件；以微波腔体1、金属搅拌棒2以及金属柔性丝带4正中心为几何坐标原点(0，0)，金属搅拌棒2的边长尺寸为15cm×1cm，金属柔性丝带4的边长尺寸为15cm×1cm，微波腔体1的边长尺寸为40cm×40cm，在微波腔体1的左边边上的中间，矩形波导3的尺寸为25cm×10cm，并且以最左边为波导馈入端口。

2)定义金属搅拌棒2、金属柔性丝带4以及微波腔体1内介质空气的电磁特性，其参数如表2所示，

表2：金属搅拌棒、金属柔性丝带及空气的电磁特性参数

3)对建立的微波场混沌搅拌器几何模型进行网格划分，本实施例中使用最大单元尺寸为24[mm]，最小单元尺寸为0.0812[mm]的三角形对图2中几何模型进行剖分，得到的网格如图3所示。

(2)定义微波搅拌器所在几何区域的控制方程和边界约束方程，本实施例中几何区域的控制方程为：

其中，E为电场强度(矢量)，k₀自由空间的波矢、μ_r为相对磁导率，ε_r为相对介电常数，ε₀为真空中的介电常数，σ为电导率，ω为角频率，j为虚数

边界约束方程为：

其中，H为磁场强度，E为电场，n为法向矢量，μ₀为真空中的磁导率，μ_r为相对磁导率，ε_r为相对介电常数，ε₀为真空中的介电常数，σ为电导率，ω为角频率；

本实施例中端口激励电磁波为TE₁₀波，其频率为2.45GHz、功率为1W。

(3)建立微波搅拌器中金属搅拌棒2以及金属柔性丝带4的随机运动函数模型；通过MATLAB数学分析软件对金属搅拌棒2以及金属柔性丝带4的随机运动进行定义，其中随机函数代码为a＝0:0.01:1,b＝randn(size(a))；通过plot(a,b)绘图得到的随机图形如图4所示；金属搅拌棒三条边在空间框架坐标系中x轴方向上的运动方程为：x＝X·cos(2πt)-Y·sin(2πt-X)-int(t)，金属搅拌棒三条边在空间框架坐标系中y轴方向上的运动方程为：y＝Y·cos(2πt)+X·sin(2πt-Y)+int(t)；金属柔性丝带三条边在空间框架坐标系中x轴方向上的运动方程为：x＝int(t)·X·cos(2πt)-Y·sin(2πt-X)+int(t)，金属柔性丝带三条边在空间框架坐标系中y轴方向上的运动方程为：y＝int(t)·Y·cos(2πt)+X·sin(2πt-Y)+int(t)；其中X，Y为搅拌部件在材料框架坐标系下对应的坐标值，int(t)为随机变量，将随机函数(a，b)对应的值通过插值函数得到，t为时间变量，a在0到1之间取值，b＝randn(siza(a))。

(4)求解得到微波场混沌搅拌器内的动态电磁场分布，包括以下分步骤：

(41)初始化，由前面网格划分得到的节点坐标，依次对初始化时刻的网格节点上对应的电场分量和网格移动的位移值赋零值，设定激励电磁波为TE₁₀波，其频率为2.45GHz、功率为1W，并设定动态电磁场求解程序运行时间为1s，每隔0.01s输出保存一次结果，即时间间隔为0.01s；

(42)函数更新，将当前时刻代入步骤(3)获得的金属搅拌棒2以及金属柔性丝带4随机运动函数模型，依据x＝X·cos(2πt)-Y·sin(2πt-X)-int(t)和y＝Y·cos(2πt)+X·sin(2πt-Y)+int(t)确定金属搅拌棒2在空间框架坐标系下相对于初始化时刻位置的位移值，依据x＝int(t)·X·cos(2πt)-Y·sin(2πt-X)+int(t)和y＝int(t)·Y·cos(2πt)+X·sin(2πt-Y)+int(t)确定金属柔性丝带4在空间框架坐标系下相对于初始化时刻位置的位移值，并将金属搅拌棒2在空间框架坐标系下相对于初始化时刻位置的位移值与金属柔性丝带4在空间框架坐标系下相对于初始化时刻位置的位移值作为当前时刻搅拌部件网格节点位移值，由于搅拌部件网格的移动导致周围空间网格的被动变形，根据文献《基于调和函数的移动网格有限元方法》中的方法可以得到周围网格被动变形后的网格节点位移值，进而完成一次函数更新；

(43)网格移动，将步骤(42)获取的当前时刻网格节点位移值加载到初始化网格节点坐标上，完成网格移动；

(44)网格畸变判断，本实施例中以三角形的最短边与最长边的比值作为网格质量，加入的判断条件(即设定值)为0.2，判断最小网格质量是否小于0.2，若小于0.2，网格发生畸变，进入步骤(47)；若不小于0.2，网格没有发生畸变，进入步骤(45)；

(45)电磁场分布计算，根据步骤(43)获得的移动之后的网格位置、前面给出的边界约束方程、激励电磁波数据以及上一时刻得到的微波场混沌搅拌器的电磁场计算结果等，利用有限元法求解前面给出的控制方程获得当前时刻微波腔体1内的电场分布；再依据电场与磁场的关系(麦克斯韦方程组)，获得当前时刻微波腔体1内的磁场分布；

(46)判断是否达到设定运行时间，若达到设定运行时间，得到微波腔体1内不同时刻的电磁场分布情况(即动态电磁场分布)，任务完成；若没有达到设定时间，在当前时刻基础上增加设定时间间隔0.01s作为下一个当前时刻，然后返回步骤(42)；

(47)网格重构，由初始化时刻金属搅拌棒2以及金属柔性丝带4坐标值与当前时刻两者位移值分别相加得到当前时刻金属搅拌棒2以及金属柔性丝带4在空间框架坐标系中的坐标值，再根据网格划分原理对搅拌器的二维几何模型进行剖分，得到重构后的几何模型网格，并以网格畸变时刻为当前时刻返回步骤(42)。

通过实施例1和实施例2获得带有非固定金属搅拌棒2的搅拌器微波腔体以及带有非固定金属搅拌棒2和金属柔性丝带4的搅拌器微波腔体内的动态电磁场分布，并从中提取不同时刻的电场分布，得到如图6和图7所示【(a)为带有非固定金属搅拌棒的搅拌器，(b)为带有非固定金属搅拌棒和金属柔性丝带的搅拌器】，其中图6给出的搅拌器微波腔体和矩形波导内电场随时间变化的分布情况，图7给出的是微波腔体内同一位置电场随时间变化情况，通过比对可以明显看出带有非固定金属搅拌棒以及带有非固定金属搅拌棒和金属柔性丝带的搅拌效果均得到了改善，且带有非固定金属搅拌棒和金属柔性丝带的搅拌效果更加均匀。

通过上述数值仿真分析，带有非固定式搅拌部件具有运动随机性，从而使得能量分布更加均匀，能够有效避免微波加热食品过程中的过热效应，从而验证了设计非固定式搅拌部件的微波搅拌器的可行性，为进一步研发带有非固定式搅拌部件的微波搅拌器提供可靠依据和有效数据。

本领域的普通技术人员将会意识到，这里的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。