顺抗磁性混合气体在梯度磁场中流动的蒙特卡罗直接模拟

D0I:10.13374/j.issn1001-053x.2006.05.036 第28卷第5期 北京科技大学学报 Vol.28 No.5 2006年5月 Journal of University of Science and Technology Beijing May 2006 顺抗磁性混合气体在梯度磁场中流动的 蒙特卡罗直接模拟 蔡军)王立1)吴平2) 尚晓航3) 1)北京科技大学机械工程学院，北京1000832)北京科技大学应用科学学院，北京100083 3)北京联合大学信息学院，北京100080 摘要以氧气和氮气的混合气体为研究对象，运用蒙特卡罗直接模拟方法模拟了顺抗磁性混合 气体在高梯度磁场中的流动情况，并把棋拟计算结果与已有的实验结果之间进行了对比分析，模 拟结果表明：磁场强度与梯度乘积由100T2·m1增加到1000T2·mˉ1时，氧体积分数增由 0.12%增加到1.89%；混合气体温度由283K升高到333K时，氧体积分数增量由0.4%减小到 0.1%:混合气体中氧气的初始体积分数由20%增加到30%时，氧体积分数增量由0.12%增加到 0.59%:压力的变化对氧体积分数的变化几乎没有影响. 关键词梯度磁场；蒙特卡罗直接模拟方法；顺磁性；仿真 分类号0552.32：V211.1*2 早在19世纪，人们就发现了氧气和氮气具有 时要求网格尺寸小于平均自由程，则碰撞只会在 不同的磁性：氧气为顺磁性磁介质，标准状态下体 同一网格中的分子间发生，这样就可以将分子的 积磁化率为146×106；而氮气则呈抗磁性，标准 运动和碰撞分开处理.因而蒙特卡罗直接模拟方 状态下体积磁化率为-0.58×106.根据这一特 法将分子的碰撞和运动在很短的时间步长内进行 性，氧气和氨气在梯度磁场中会受到不同的磁场 解耦，从而直接跟踪分子的运动情况.在模拟初 力作用，氧气所受磁场力指向磁场强度增大的方 始，根据物理条件在计算区域内布置模拟分子并 向，而氮气受到的磁场力指向磁场强度减小的方 给定其位置和速度，在每一时间步长内分别计算 向，朱重光就非均匀磁场对原子的作用力进行了 运动和碰撞.处理运动时还包括分子进入或逸出 分析和推导，并讨论了横向非均匀磁场对氧气和 计算区域（对应于进出口边界条件），与壁面的相 氨气的分离作用，论证了利用梯度磁场提取富氧 互作用；碰撞的处理主要是碰撞对的选择和碰撞 的可行性12].利用氧气和氮气不同的磁性制造 截面积的计算. 的氧分析仪器已经在工业上得到了广泛的应用. 本文以氧气和氮气的混合气体为研究对象， 一些学者还研究了利用梯度磁场来控制气体的流 运用蒙特卡罗直接模拟方法，在计算过程中加入 动以改善燃烧、加速呼吸和促进晶体生长等[36]. 了梯度磁场的影响作用，构造了顺抗磁性混合气 蒙特卡罗直接模拟方法是基于稀薄气体流动 体在高梯度磁场中流动的物理模型和数学模型， 数值模拟研究的需要发展起来的一种新的计算气 研究了混合气体在高梯度磁场中的流动情况，并 体流动的方法，由Bird等人在20世纪60年代提 对影响因素进行了分析 出的，其实质就是用适当数目的仿真分子代替大 1模拟方法 量真实气体分子，用计算机模拟气体分子的运动 过程.在实际情况中，分子的运动和碰撞是耦合 氧气在梯度磁场中所受的磁场力远大于氨气 在一起的，这在数值模拟中很难实现.如果限制 所受的磁场力，相差约200倍.作用在单位体积 时间步长小于平均碰撞时间，则在平均意义上一 顺磁性物质上的磁场力为： 个分子在一个时间步长内最多只有一次碰撞；同 F=XB dB (1) o dx 收稿日期：200503-21修回日期：200603-16 作者简介：菜军(1977一)，男，博士研究生；王立(1956一），男，教 式中，F为磁场力，40为真空磁导率，X为体积磁 授，博士 化率，B为磁感应强度，dB/dX为场强梯度.对

第 2 8 卷 第 5 期 2 0 06 年 5 月 北 京 科 技 大 学 学 报 J o u r n a l of U n ive 比 i ty o f Sc i e n Ce a n d T e c h n o l昭 y Be Ui gn V 0 1 . 2 8 N o . 5 M a y 2 0 0` 顺抗磁性混合气体在梯度磁场中流动的 蒙特卡罗直接模拟 蔡 军` ) 王 立 1) 吴 平 2) 尚 晓航 3) 1 ) 北京科技大学机械工程学院 , 北京 10 0 0 8 3 2 ) 北 京科技大学应用科学学院 , 北 京 10 0 0 5 3 3 ) 北京联合大学 信息学院 , 北京 1 0 0 0 8 0 摘 要 以氧气和氮 气的混 合气体 为研究对象 , 运用蒙特卡 罗直接模拟方 法模拟 了顺 抗磁性 混合 气体在高梯度磁场 中的流 动情况 , 并把模拟计算结 果与 已有 的实验结 果之 间进行 了对 比分析 . 模 拟结果 表 明 : 磁 场强 度 与 梯度 乘 积 由 10 0 T Z · m 一 ’ 增 加到 1 0 00 护 · m 一 ’ 时 , 氧体 积 分数 增 量 由 0 . 12 % 增 加到 1 . 89 % ; 混合 气体温 度 由 2 83 K 升高到 3 3 K 时 , 氧体积 分 数增量 由 0 . 4 % 减 小到 0 . 1 % ; 混合气体中氧气的初始体积分数 由 20 % 增 加到 30 % 时 , 氧体积 分数增 量 由 0 . 12 % 增 加到 0 . 5 9 % ; 压力的变化对氧体积分数的变化几乎没 有影响 . 关扭词 梯度磁场 ; 蒙特卡罗直接模拟方法 ; 顺 磁性 ; 仿真 分类号 0 5 5 2 . 3 + 2 ; V2 1 1 . l + 2 早在 19 世纪 , 人 们就发现 了氧气和 氮气具 有 不 同的磁性 : 氧气为顺磁性磁介质 , 标准状 态 下体 积磁 化率为 14 6 x 10 一 6 ; 而 氮气则 呈抗磁 性 , 标准 状态 下体积磁 化率为 一 0 . 5 8 x l0 ’ 6 . 根据这 一 特 性 , 氧气和 氮气在梯度 磁 场中会 受到 不 同的磁 场 力作用 , 氧气所 受磁 场力指 向磁场强 度增大的方 向 , 而 氮气受到 的磁 场 力 指 向磁 场强 度 减小 的方 向 . 朱重 光 就非均匀 磁场对原子的作用力进行了 分析和推导 , 并讨论 了横 向非均 匀磁 场对氧气和 氮气的分离作用 , 论 证 了利用 梯度磁 场提取富氧 的可行性【’一 利用 氧气和 氮气不 同的磁 性制造 的氧分析仪器 已 经 在 工 业 上 得到 了广 泛 的应 用 . 一些 学者还研究 了利用梯度磁场来控制气体的流 动以改善燃烧 、 加速呼 吸和 促进晶体生长等〔3司 . 蒙特卡罗直接模拟 方法 是基 于 稀薄气体流动 数值模拟 研究的需 要发展起来的一种新的计算气 体流 动的方法 , 由 iB dr 等人在 2 0 世纪 6 0 年代提 出的 . 其实 质就是用适 当数 目的仿真 分子代替大 量真实气体分子 , 用 计算机 模拟 气体分子 的运 动 过程 . 在实际情 况 中 , 分子 的运 动和 碰 撞 是 藕合 在一 起 的 , 这 在数值模拟 中很 难实现 . 如 果限制 时间 步长 小于 平 均碰 撞 时间 , 则 在平均 意 义上 一 个分子在 一个 时 间 步 长 内 最 多 只有 一 次 碰 撞 ; 同 收稿B 期 : 2 0 0 5习3 一 1 修 回 B 期 : 2 0 0 6一3 一 1 6 作者简介 : 蔡军 ( 1 9 7 一 ) , 男 . 博士研究生 ; 王 立( 19 56 一 ) . 男 , 教 授 , 博士 时要 求网格尺 寸 小于 平均 自由程 , 则 碰 撞 只 会 在 同一 网格中的分子 间发生 , 这 样就 可以 将分子 的 运 动和 碰撞分开 处理 . 因而蒙特卡罗 直接模拟 方 法 将分子的碰撞 和运 动在很 短的时间步长 内进行 解藕 , 从 而 直接跟 踪分子 的运 动情况 . 在模拟 初 始 , 根据物理 条件在计算 区域 内布置 模拟 分子 并 给定其位置和 速 度 , 在每一 时间步 长 内分别计算 运 动和 碰撞 . 处理运 动时还 包括 分子 进 入或逸 出 计算区域 (对 应 于进 出 口 边 界 条 件 ) , 与壁 面 的相 互 作用 ; 碰撞 的处理 主要 是碰 撞对 的选 择和 碰撞 截面积 的计算 . 本文 以氧气和 氮气的混 合气体 为研 究对象 , 运 用蒙特卡罗 直接模拟 方 法 , 在 计算过 程 中加 入 了梯度磁场的影 响作用 , 构 造 了顺抗 磁性 混 合气 体在高梯度 磁 场中流 动的物理 模型 和 数学 模型 , 研 究 了混合气体 在高梯度 磁 场中的流 动情况 , 并 对 影 响因 素进行了分析 . 1 模拟方法 氧气在梯度磁场中所 受的磁场力远 大于氮气 所受的磁 场力 , 相差 约 2 0 0 倍 . 作用 在单位体 积 顺 磁性物质上 的磁场 力为 : 一 X , 、 厂 = — 廿 ( l ) 一XB ` 一d . 刁U 式 中 , F 为磁 场力 , 产。 为真 空 磁导 率 , x 为体积 磁 化率 , B 为磁感应 强 度 , d B / d X 为场 强 梯度 . 对 DOI: 10. 13374 /j . issn1001 -053x. 2006. 05. 036

·476· 北京科技大学学报 2006年第5期 于特定的粒子，磁场力大小决定于B和dB/dX 确定了分子碰撞后的相对速度，可得到碰撞 在300K时，1T和120T·m1的磁场作用在1mol 后单个分子的运动速度为： 氧气上的磁场力约为4.2N,具体到每个分子上 Vi-G+m&.Vi-G-mI (6) 约为7×10-24N.氧分子的质量为1026数量级， 因此磁场力可以使氧气分子达到102~103数量 式中，g'为碰撞后分子的相对速度，G= 级的加速度 1.1程式化处理 m,y,+mbV为碰撞分子的质心速度，m,= m+mb (1)网格划分.网格用来选择碰撞对和统计 mnb为折合质量 流场，其线性尺度△r与流场宏观量梯度的标尺 m,+mb 长度比较应该是小量，一般取△r=入/3，其中入= 1.3仿真分子碰撞对的抽样方法 1/(2πd2n),为平均自由程. 分子碰撞对的抽样方法主要有时间计数器 (2)选取时间推进步长.蒙特卡罗直接模拟 (TC)方法、非时间计数器(NTC)方法、随机取样 方法的本质是在很小的时间步长△t内将分子的 频率方法以及Baganoff和McDonald抽样方 运动与分子间的碰撞解耦.因此△：应该比局部 法8].本文采用NTC方法，先给出在时间步长 平均碰撞时间小得多，以保证在这一时间间隔内 △t内单元网格中的仿真分子碰撞次数N,: 所有仿真分子是“自由”的，即应有： N,-7Nmn(oTR)mAt (7) △t<1/u (2) 式中，Nm为单元网格中仿真分子的数目，n为单 式中，”为气体分子的平均碰撞频率，由下式确 元网格中的分子数密度 定： 根据气体分子动力学理论，网格内仿真分子 =2 /2xT/m/(na) (3) 对的碰撞几率函数P可表示为： (3)选取仿真分子数.合理的仿真分子数的 P(g）=(oTg)x OT8 (8) 选取必须兼顾计算效率以及统计得到的宏观物理 量的真实性.考虑到在每一时间推进步长上每一 式中，cT=xd2,为分子的碰撞截面.本文采用的 网格内得到的宏观址都具有统计意义，Bird建议， 分子模型为变径硬球(VHS)模型，VHS模型中分 在每一个单元网格中布置20~30个模拟分子[] 子直径可以表示为相对速度的函数]： 1,2二元弹性碰撞模型 蒙特卡罗直接模拟方法只考虑二元碰撞的情 d-d(2KTl(m (9) T(2.5-w) 况，并且认为分子间的碰撞是弹性碰撞.由于碰 式中，下标ref表示参考值，m,为折合质址，g为 撞的时间非常短，在碰撞瞬间可以忽略磁力的影 相对速度，“为气体的粘性指数.由此可见，碰撞 响，因此碰撞的双分子体系必然遵守能量、动量以 截面T是分子间相对速度的函数. 及质量守恒定律，根据质量及动量守恒定律，无 把P与在区域上均匀分布的随机数R作 论气体分子间发生何种类型的碰撞，碰撞后分子 比较，如果R<P,认为它们形成了碰撞分子 体系的质心速度都将保持不变.根据能址守恒可 对，并根据二元碰撞理论，计算出气体分子碰撞后 以推导出碰撞前后分子间的相对速度大小不变， 的速度.如果R≥P,则重新选取分子对. 只是方向发生了变化.记n=(n1,n2,n3)为碰撞 1.4边界处理 后分子相对运动速度方向的单位矢量，由如下方 (1)分子与固体壁面的相互作用 法确定： 对于分子在固体壁面的反射模型，目前比较 n1=COSX1 常用的有两种，即镜面反射模型和完全漫反射模 n2=sinxicosx2 (4) 型.本文采用漫反射模型：假定离开壁面的分子 n3=sinxisinx2 以平衡的分布散射，气体分子与固体壁面碰撞后 式中，X1和X2为碰撞参数，由随机抽样确定，其 的速度与碰撞前的速度完全没有关系，反射分子 抽样方法为： 作为一个整体，在速度上满足半空间的Maxwell cosX1=2R1-1,X2=2πR2 (5) 速度分布1o),因此分子反射后法向速度V,的抽 式中，R1,R2为(0,1)区间上均匀分布的随机数 样几率函数为：

· 47 ` . 北 京 科 技 大 学 学 报 2 0 ` 年第 s 期 于特定的粒子 , 磁 场力大小决 定于 B 和 d B / d X . 在 3 0 0 K 时 , I T 和 12 0 T · m 一 ’ 的磁场作用 在 1 m o l 氧气上 的磁 场力约 为 4 . 2 N , 具 体到 每个分子上 约 为 7 X 10 ’ 24 N . 氧分子的质量 为 10 一 2 6数量级 , 因此 磁场力可以使氧 气分子 达 到 10 2 一 10 3 数量 级的 加速 度 . 1 . 1 程式化处理 ( l) 网格 划分 . 网格用 来选 择碰撞对和统计 流场 , 其线性 尺度 △: 与 流场宏观 量 梯度 的标尺 长度 比较应该是 小量 , 一 般取 △: = 几3/ , 其中 又= 1/ 盛耐 2 , ) , 为平均 自由程 . ( 2) 选取时间推进 步长 . 蒙特卡罗直接模拟 方法 的本质是 在很 小的时间步 长 △t 内将分子的 运 动与分子 间的碰撞 解祸 . 因此 △t 应 该 比局部 平均碰 撞时间 小得 多 , 以保证 在这 一 时间间隔 内 所有仿真分子是 “ 自由 ” 的 , 即应 有 : △t ( l / 。 ( 2 ) 式 中 , 。 为气体分 子 的平均碰 撞 频 率 , 由下 式 确 定 : 。 一 ( 2 丫五不丽 ) / ( ! 、 ) ( 3 ) ( 3) 选取仿真分子数 . 合理 的仿真分子数的 选取必须兼顾计算效率以及 统计得到 的宏观 物理 量的 真实性 . 考虑到 在每一时间推 进步长 上每一 网格内得到 的宏观量都具有统计意义 , iB dr 建议 , 在每一 个单元 网格中布置 20 一 30 个模拟 分子〔7 ] . 1 · 2 二元弹性碰摘模型 蒙特卡罗直接模拟 方法 只考虑 二元 碰撞 的情 况 , 并且 认为分 子 间的碰 撞是 弹性碰撞 . 由于 碰 撞 的时间非常短 , 在碰 撞瞬 间可以 忽略磁力 的影 响 , 因此 碰撞的双分子体 系必然遵守能量 、 动量以 及 质量守恒定律 . 根据质量及 动量守 恒定律 , 无 论气体分子间发生何种类型 的碰撞 , 碰撞 后 分子 体系的质心速度都将保持不变 . 根 据能量守 恒可 以推导出碰 撞前后分子 间的相 对速度大小不变 , 只是方 向发生 了变化 . 记 n = ( n , , n Z , , 3 )为碰撞 后分子 相对 运 动速 度 方 向的单位矢 量 , 由如 下 方 法确 定 : 确 定了分子 碰撞 后 的相对速度 , 可得 到碰 撞 后 单个分子 的运 动速度为 : V 二二 G + 贰 g V 云= G 一 二` g 产儿 b ( 6 ) 式 中 , g ` 为 碰 撞 后 分 子 的 相 对 速 度 , m . V 。 + m 卜 V 卜 、 , _ 、 * 一 , 、 ~ , , ~ 、 、 ~ 一 一不硕不戒一 ” ` , 分 , 。 。 心 。 二 , m r = m a m b m a + 1 。 3 为折合质量 . 阴 b 仿真分子碰抽对的抽样方法 分子碰 撞对 的抽样方 法 主 要 有时 间计数 器 ( T C )方法 、 非时间计数器 ( N T C ) 方 法 、 随 机取 样 频 率 方 法 以 及 B ag a on ft 和 M c oD an ld 抽 样 方 法 8[] . 本文 采用 N T c 方 法 , 先给 出 在时间步 长 △t 内单元 网格中的仿真分子碰 撞次数 N , : N ` 一 鲁N m , ( 。 T g ) anI x △, ` 式中 , N m 为单元 网格中仿真分子 的数 目 元网格 中的分子 数密度 . ( 7 ) n 为单 根据气体分子动力学理 论 , 网格 内仿真分子 对的碰 撞几率函数 尸 。 ol 可表示 为 : 尸col ( g ) = ` T g ( 。 T g ) ma x ( 8 ) 式中 , 。 T = 二 d ’ , 为分子 的碰 撞截面 . 本文 采用的 分子 模型为变径硬球 ( V H )S 模型 , v H S 模型 中分 子直径可以表示为相对速度的函数9[] : d = d er f 尸 ( 2 . 5 一 。 ) ( 9 ) ( 4 ) osxzixnz 义XX cos · sln · 一 式中 , x , 和 x : 为 碰撞 参数 , 由随机 抽样确 定 , 其 抽样方 法为 : co s X一 2 R I 一 l , X Z = 2 兀 R Z ( 5 ) 式 中 , R : , R : 为 ( 0 , l) 区 间上 均匀 分布 的随机 数 . 式中 , 下 标 er f 表示 参考值 , m r 为折合质量 , g 为 相对 速度 , 。 为气体的粘性指数 . 由此 可见 , 碰撞 截面 。 T 是分子间相对速度的函数 . 把 尸col 与在区域上 均 匀分布 的随 机 数 R 作 比较 , 如 果 R < 尸col , 认 为它 们形 成 了碰 撞分 子 对 , 并根据二元 碰撞理论 , 计算出气体分子碰撞后 的速 度 . 如果 R ) p col , 则重 新选取分子对 . 1 . 4 边界处 理 ( 1) 分子 与固体壁面 的相互 作用 . 对于分子在 固体壁 面 的反射模型 , 目前比较 常用的有两种 , 即镜面 反射模型 和 完全漫 反射模 型 . 本文采用 漫 反射模型 : 假定离开壁 面 的分子 以平衡的分布散 射 , 气体分子 与固体壁 面 碰撞 后 的速 度与碰撞 前的速 度 完 全 没有关 系 , 反 射分子 作为一 个整体 , 在速 度上 满足 半空 间的 M ax w el 速度 分布[ ’ “ ] , 因此 分 子反 射后 法 向速 度 v 。 的抽 样 几率 函数 为 :

Vol.28 No.5 蒸军等：顺抗磁性混合气体在梯度磁场中流动的蒙特卡罗直接模拟 ·477· exp(-B2v2) (10) 绝对值Bd的最大值出现在X=18.8mm和 √π dx 式中，=√m/(2kT).由于反射分子关于半空 X=-188mm的地方：在其余X值处， 间的散布是各向同性的，反射分子总是沿着壁面 的值急剧下降.由式(1)可知氧分子所受磁力与 外法向一侧的半空间运动，因此必须按照在半空 间内各向同性散布对反射分子的方位角进行抽 |成正比，因此选取为最大值附近 样，用中，中表示反射分子的方位角，其中中为子 的区域作为模拟时的梯度磁场区域，而把磁极间 午角(0≤p≤π/2)，中为水平角(0≤≤2π)，其抽 隙内的其他区域视为均匀磁场区域.本文选取的 样方法为： 梯度磁场区域沿X轴方向的模拟宽度为100μm. 少=2πR,cosp=1-R (11) 根据流动的对称性，只对一侧的梯度磁场区域进 式中，R为(0,1)区间上均匀分布的随机数 行模拟 (2)分子通过外边界的处理， (a)流入边界.单位时间内通过单位面积由 外边界进入到计算区域的分子数为： N,=-1_noe 2/元8eexp(-s2os29)+ √元scos0[1+erf(scos0)]} (12) 式中，Bm=√m/2kTo,s=|Vm|Bm为气体分子 图1磁极布置形式 速度比，n为自由来流分子数密度，日为自由来 Flg.I Layout of magnetic poles 流宏观速度V。与表面微元外法向向量的夹角， 1.2 600 ert(z)=二exp(-x2)dr为误差函数.时间 B 1.0 0气e 步长△：内从外边界进入到计算单元的分子先被 400 转换为仿真分子，然后再赋子其位置和速度. 300 (b)流出边界，当被跟踪的仿真分子的运动 L)/CXP/8P) 200 轨迹到达外边界以外区域后，就认为该仿真分子 02 100 Bx(dB/dx 已经逃逸到计算区域外，并在计算机中删除该仿 0 -20 -10 0 10 20 真分子的有关记录 距离Wmm 2 物理问题简述与计算条件设定 图2气隙碰场沿X轴方向上的碰场分布 Flg.2 Field intensity distribution of the alr-gap magnetic fleld 整个磁场区域由两块矩形永磁体形成的气隙 along the X direction 磁场构成，两块永磁体异极相对，其磁极布置形式 混合气体以流其Q,从磁场间隙左侧端面流 如图1所示，6=1mm为两块磁体之间的间隙宽 入，以流址Q2从磁场间隙右侧端面流出，且Q2 度.永磁体的结构尺寸为L×W×H为79mm× <Q1,其余部分气体从磁场间隙的两侧以流其 38mm×30mm.两块永磁体之间的间隙不仅构 Q3与Q4流出，如图3所示 成了气隙磁场，而且还是气体流动的通道，气体从 侧面沿Z轴进入磁体间隙. 以N35钕铁硼永磁体为例，按照图1所示的 磁极布置形式，用Ansys软件对其进行磁场计算. 图2所示为气隙磁场沿X轴方向的磁场分布结 果.从图2可以看出，气隙磁场在-15mm≤X≤ +15mm的范围内几乎是均匀的，当X>15mm 围3气体在磁场间隙内流动的宏观示意围 Flg.3 Gas flows in the magnetic gap 和X<15mm时，磁感应强度开始下降，约在磁场 间隙边界处处降为零，磁感应强度与梯度乘积的

V oj 。 2 8 N O 。 5 蔡军每 : 顺抗磁性混合气体在梯度磁场 中流动 的萦特卡罗直 接模拟 , , 、 , 、 _ { 理Y 丝 l _ _ _ _ , 。 2 二 r Z 、 f ( V , ) = } J 下 竺 ) e x p ( 一 月 ` V 盘) ( 10 ) 、 几 , 一“ r 、 。 · 。 , 式中 , 月= 了石兀丽而 . 由于 反 射分子关 于 半空 间的散布是各向同性 的 , 反 射分子总 是沿 着壁 面 外法 向一侧的半空 间运 动 , 因此 必 须按照 在半空 间 内各向同性散 布对 反 射分子 的方位角进 行抽 样 . 用 价 , 必表示反射分子的方位角 , 其中 协为子 午 角(0 簇 协簇 刁2) , 沪为水平角 (0 镇 沪簇 2动 , 其抽 样方法 为 : 沪= 2 二 R , co s 笋= l 一 R ( 1 1 ) 式 中 , R 为( 0 , l) 区间上 均匀分布的随机数 . ( 2) 分子通过 外边 界 的处理 . ( a) 流 入边 界 . 单位时间 内通 过 单位面 积 由 外边 界进入到计算区域 的分子数为 : 绝 对值 { 里哗 } 的最 大值出现 在 x 一 1 8 . 8 ~ 和 1 0 人 } 二 ` 二 。 一 二 、 人 , , 、 ` , } B d B } X = 一 18 . s m m 的地 方 ; 在 其余 x 值处 , _ 尸兰吕! _ . 。 二 . 二 ~ 二 , 、 一 , J , ~ , 、 , 二 、 ` · ~ 一 ’ } d X } 的值急剧 下降 . 由式 ( l) 可知氧分子所受磁 力 与 } 卫华 } 成正 比 , 因此 选 取 1 旦;华 { 为最大值附近 l d 入 } ! O 人 } 的区域作为模拟 时的梯度 磁场 区域 , 而把磁 极 间 隙内的其他区域视 为均匀 磁场 区域 . 本 文选取的 梯 度磁场区域沿 x 轴方 向的模拟 宽度为 10 拌m . 根据流动的对 称性 , 只对 一 侧的梯度 磁场 区域进 行模拟 . N , 一 众箭{ 二p` 一 ’ 。 osz 代’ ` 人 : 。 os 。 [ 1 + e fr ( : c o s 。 ) ] } ( 1 2 ) 式中 , 月co = 了石万丽呱 , : 二 } v 。 }夕co 为气体分子 速 度 比 , n co 为 自由来流 分子 数密度 , 口 为 自由来 流宏观 速 度 V co 与表面 微元 外法 向向量 的夹角 , - 一 一一 Z \ X 图 1 磁极布里形式 F ig . I L叮o u t o f m a g n e t l e OP Ies ǎ 一.日 · 件兰g弓ù与à 0 nU on ù 0 lnj on 七J ù 4 内l`, ` er f( 二 ) 一 鲁} ’ 。 x p ( 一 二 , ) d二 为误 差 函 数 . 时 间 了沉 J O 步长 △ t 内从外边 界进入 到计算单元 的分子先被 转换为仿真分子 , 然后 再赋 予其位置 和速 度 . ( b) 流 出边 界 . 当被跟 踪的仿真分子 的运 动 轨迹 到达 外边界 以 外区 域后 , 就认为该仿真分子 已经逃 逸到 计算区 域外 , 并在计算机 中删除该仿 真分子的有关记录 . 窝 。名 } B 火 ( d B厄 X ) 1 又X . _ _ 咱 1明 — 一 10 0 10 2 0 距离 x/ m m 图甩 64 00 八“,山 . 侧烈硕招 2 物理问题简述与计算条件设定 整个磁 场区域 由两块矩形 永磁 体形成的气隙 磁场构成 , 两块永磁体异极 相对 , 其磁极布置形 式 如 图 1 所示 , 占= l m m 为两 块磁 体之 间 的间隙宽 度 . 永磁 体的结构尺寸为 L x w x H 为 79 m m x 3 8 m m X 3 0 m m . 两 块永 磁体 之 间的间 隙不仅构 成了气隙磁场 , 而 且还 是气体流 动的通道 , 气体从 侧面 沿 Z 轴进入磁体间 隙 . 以 N 3 5 钱铁硼永磁 体为例 , 按照 图 1 所示 的 磁 极布置形 式 , 用 A sn ys 软件对其进行磁场 计算 . 图 2 所示为气隙磁 场沿 X 轴方 向的磁 场分布结 果 . 从图 2 可以看出 , 气隙磁 场在 一 巧 m m 镇 X 簇 + 1 5 m m 的范围内几乎是 均 匀的 , 当 X > 15 m m 和 X < 15 m m 时 , 磁感应强 度开 始下降 , 约在磁场 间 隙边界处 处降为零 . 磁 感应强 度与梯度 乘积 的 气隙磁场沿 X 轴方向上的磁场 分布 lF g . 2 a I O l g F i e ld i n t e sn i t y d i s t d b . t场 n o f t h e a i r · g a p ma 忍n . ti e n e ld t he X d i r e c ti o n 混 合气体以流量 Q ; 从磁 场间隙左 侧端面 流 入 , 以流 量 Q : 从 磁 场 间 隙右侧端 面 流 出 , 且 Q Z < Q l , 其 余 部分气体 从 磁 场 间 隙的两 侧 以流 量 Q 3 与 Q ; 流 出 , 如图 3 所 示 . 洲一十刊一刁一刊一号, / / / / / / / 应 、 、 、 入 、 曰 = 。 ’ 厂 ` 图 3 气体在磁 场间隙内流动 的史观示意 图 F lg . 3 G as fl o w s i n the am护 e ti e ga p

·478 北京科技大学学报 2006年第5期 根据对称性可知Q3=Q4·从磁场间隙沿X 25 一●一模拟值 轴方向两侧流出的气体将会穿过梯度磁场区域， ▲文献111实验值 由于存在场强梯度，氧气分子将受到指向Z轴方 2.0 向并与之垂直的磁场吸引力而留在磁场间隙内， 1.5h 而氮气分子则受到相反方向的作用力加速流出磁 1.0 场区域.可以看出，梯度磁场区域的作用类似于 的 一个具有选择性的“筛子”，氧气分子不易穿过该 0.5 筛子而氮气分子却相对更容易穿过，从而使得磁 2004006008001000 场间隙内气体的氧的体积分数升高 B*(dB/dX)/(T2.m-) 本文的计算工质为氧氮混合气体，网格尺寸 为108m,每个网格中的仿真分子总数为25，氧 围5氟的体积分数变化量随兴的变化情况 氨分子的个数按照各自组分的体积分数进行分 Flg.5 dx 配.网格单元中分子碰撞对的选取采用NTC方 法.根据式(2)与式(3)确定时间步长为10~Ⅱs, 3,2模拟计算结果 当计算区域内的仿真分子数围绕某个值上下波动 图6所示为磁场间隙内氧气的体积分数变化 时，认为流场稳定.计算中，进出口采用宏观上有 量随温度的变化情况.其参数条件为P=101325 流量通量的流体边界条件；分子与上下磁极壁面 Pa400 ml.min 的碰撞采用漫反射模型；分子之间的作用势采用 始氧的体积分数为24%.由图6可知，随着混合 变径硬球(VHS)模型.由于氨气分子的磁化率远 气体温度的不断升高，磁场间隙内的氧的体积分 远低于氧气分子的磁化率，因此不考虑梯度磁场 数变化量越来越小.当温度值从283K升高到 对氮气分子的作用.按照蒙特卡罗直接模拟方法 333K时，氧的体积分数变化量由0.4%减小到 的基本思想，编制出了计算流程图，如图4所示. 0.11%.这是因为氧气分子的磁化率与温度的平 扦始 方成反比，而由式(1)可知氧气分子所受磁力与其 磁化率成正比，当温度升高时，氧气分子受到的磁 设置模拟分子的初始条件及边界 条件并进行网格划分 力变小，从而削弱了梯度磁场对氧气分子运动的 阻碍作用，因此磁场间隙内气体的氧的体积分数 计算模拟分了的迁移运动及其与壁面的碰撞 ★ 也就越低 [重新标识分了 0.4 计算所有网格内的分了碰撞 计算磁场所的影啊作用 0.3 流场稳定？> TY 0.2 取样统计 输出结果 0.1 283293303313323333 (结束 气体温度K 图6氯的体积分数变化量随温度的变化情况 图4计算流程图 Flg.6 Change of oxygen concentration increment with tempera- Flg.4 Calculation flowchart tore 3计算结果及讨论 图7所示为磁场间隙内氧气的体积分数变化 量随压力的变化情况，其参数条件为T=239K, 3.1程序验证 图5所示为相同条件下，本文计算的氧的体 Bd5=100T2m1,Q1=400ml.min,初始氧 d 积分数变化量随磁场的变化情况与文献[11]的实 的体积分数为24%.由图7可知，随着压力的增 验结果的比较.由图可以看出，二者吻合的很好 大，氧的体积分数变化址的起伏不是很大，也就是

北 京 科 技 大 学 学 报 2 0 0 6 年 第 5 期 根据对称性可知 Q 3 = Q 4 . 从磁 场间隙沿 X 轴方 向两侧流 出 的气体将会 穿过 梯度磁 场区域 , 由于 存在场强梯度 , 氧气分子将受到 指向 Z 轴方 向并与之垂 直的磁 场吸引力而 留 在磁 场间隙内 , 而氮气分子则 受到相反方 向的作用力加速 流 出磁 场区域 . 可以 看出 , 梯度磁场 区域的作用 类似于 一个具有选择性 的 “ 筛子 ” , 氧气分 子 不 易穿过该 筛子而氮气分子却相对 更 容易 穿过 , 从 而 使得 磁 场间隙内气体的氧的体积 分数升高 . 本文 的计算工 质 为氧氮混 合气体 , 网 格尺寸 为 10 一 s m , 每个网格中的仿真分子 总数 为 25 , 氧 氮分子的个数按照 各 自组 分的体积 分数进 行分 配 . 网格单元 中分子 碰 撞对的选 取采用 N T C 方 法 . 根据式 ( 2) 与式 ( 3) 确定时间步长 为 10 一 ” s , 当计算区域内的仿真分子数围绕某个值上下波动 时 , 认为流场稳定 . 计算中 , 进 出 口 采用宏观上 有 流量通 量的流体边 界条件 ; 分子 与上 下磁 极 壁面 的碰撞 采用漫反 射模型 ; 分子 之 间的作用势采用 变径硬球 ( v H )S 模型 . 由于氮气分子的磁 化率远 远低 于氧气分子的磁化率 , 因此 不 考虑 梯度磁场 对氮气分子的作用 . 按照蒙特卡 罗 直接模拟 方法 的基 本思想 , 编制出了计算流 程图 , 如图 4 所示 . 2 . 5 2 . 0 一一 模拟值 ▲ 文献 1 1 实验值 势级佘彩挂酬屏芝g 20 0 4 0 0 6 0() 8 0 0 1 0 0 0 B * (胡闭 X ) /( T Z · m 一 , ) .卜 I L n 七J . O 八U 。 , 的体积分 , 变化 t 随箫 的变化 , 况 。 : · , e ` 。 o f o: y , n co n . o t . t . ou . n。 , .m 。 : w: 。 箫 设置模拟分 子的初始条件及边 界 条件并进行网格划分 计算模拟分子的迁移运 动及 其 与壁面的碰 撞 计算所有网格内的分了碰撞 计算磁场所的影响作用 3 . 2 模拟计算结果 图 6 所示 为磁场间隙内氧气的体积 分数变化 量随温度 的变化情况 . 其参数条件为 尸 = 10 1 3 2 5 _ B d B J 、 _ _ , _ : ~ 二 * , , . _ : P a , 一 吕长笋= 1 0 0 T Z · m 一 ’ , Q , = 4 0 0 m L · m i n 一 ` , ` “ , 初 dX ` 一 ” 一 ` ~ ’ 凭 且 ’ 一 ` ~ 一 “ ` . ~ 始氧 的体积分数 为 2 4 % . 由图 6 可知 , 随 着混 合 气体温 度的不断 升高 , 磁 场间 隙内的氧的体积分 数变 化量 越 来越 小 . 当温 度值从 2 83 K 升高到 3 3 3 K 时 , 氧的体积 分数变化量 由 0 . 4 % 减 小到 0 . 1 1% . 这 是 因为氧气分子的 磁化率与温 度的平 方成反 比 , 而 由式 ( l) 可知氧气分子 所受磁 力 与其 磁 化率成正 比 , 当温度升高时 , 氧气分子受到 的磁 力变小 , 从而 削弱 了梯度磁场 对氧气分子 运 动的 阻碍作用 , 因此磁 场 间隙内气 体的氧 的体积 分数 也就越 低 . 30 3 3 13 气体温度 服 32 3 3 3 3 知止9 ,`j, o0 - n 址名屏训势翻余彩芝 图 4 计算流程图 图 ` 级的体积分擞变化 t 随沮度的变化情况 F lg . ` C h a n乎 o f o x y ge n co n ce n t ar tl o n i n e 件me n t w i th t e m碑r a · lF g . 4 C a l e . l a ti o n fl ow e 加a d 3 计算结果及讨论 3 . 1 程序验证 图 5 所示 为相 同条 件下 , 本文计算的氧的体 积分数 变化量 随磁场的变化情况 与文 献 〔l lJ 的 实 验结果的 比较 . 由图可以看出 , 二者吻 合的很好 . 图 7 所示 为磁场间隙内氧气的体积 分数变化 量随 压力 的变化情况 , 其 参数 条件为 T = 2 39 K , 掣 一 10 。 T Z . m 一 , , Q 一 ; 。。 m , , . m i n 一 , , 初 始氧 d X - - 一 ` ~ ’ 气 1 ’ “ ” “ “ “ ~ · ~ , . “ ~ ~ 的体积分 数为 2 4 % . 由图 7 可知 , 随 着压 力的增 大 , 氧的体积 分数变化量 的起伏不是 很大 , 也就是

Vol.28 No.5 蔡军等：顺抗磁性混合气体在梯度磁场中流动的蒙特卡罗直接摸拟 ·479· 说压力的变化对氧的体积分数变化量的影响不是 磁场中的流动进行了计算机模拟，并把模拟结果 很大.虽然压力增大会导致氧气体积磁化率的增 与已有的实验结果进行了对比和分析，得出如下 加，但是其分子数密度也同时在增加，因此具体到 结论： 每个氧分子上的磁力并没有发生变化.所以磁场 (1)相同条件下，模拟结果与文款实验结果 间隙内氧的体积分数变化量也就不会发生变化. 吻合较好 0.5 (2)氧的体积分数增量与磁场强度和梯度的 e 乘积的绝对值成正比，当设由10下m增加 到1000T2m1时，氧的体积分数增量由0.12% 0.4 增加到1.89%. (3)氧的体积分数增量与混合气体温度成反 比，当温度由273K升高到333K时，氧的体积分 030 2000 4000 6000 数增量由0.4%减小到0.1%. 表k力，PIPa (4)压力的变化对氧的体积分数增量的影响 图7氯的体积分数变化量随压力的变化情况 不是很大. Fig.7 Change of oxygen concentration increment with pressure (5)氧的体积分数增量与混合气体中氧的初 图8所示为磁场间隙内氧气的体积分数变化 始体积分数成正比，当混合气体中氧气的初始体 量随初始氧体积分数的变化情况，其参数条件为 积分数由20%增加到30%时，氧的体积分数增址 P=101325Pa,T=239K,Bd9=1002.m1, 由0.12%增加到0.59%. dx Q1=400 mLmin-.由图8可知，随着氧气初始参考文献 体积分数的增加，磁场间隙内氧的体积分数变化 [1】朱重光.横向非均匀磁场对原子的作用，准北煤师陕学报： 量越来越大.当氧气的初始体积分数由20%增加 自然版，1996,17(1)：30 到30%时，氧的体积分数变化量由0.12%增加到 [2]朱重光.磁分制氧机理.准北煤师院学报，1997,18(1)： 0.59%.氧气初始体积分数越大时，氧分子的体 22 积磁化率也就越大，由式(1)可知它所受到的磁力 [3]Ueno S.Iwasaka M,Eguchi H,et al.Dynamics behavior of gas flow in gradient magnetic fields.IEEE Trans Magn, 就越强，因此留在磁场间隙内的氧气分子个数就 1993,29(6):3264 越多 [4]Wakayama N I,Sugie M.Magnetic promotion of combustion 0,6 in diffusion flames.Phys B,1996.216(3-4):403 [5]Wakayama N I,Wakayama M.Magnetic acceleration of in- haled and exhaled flows in breathing.Jpn J Appl Phys,2000. 0.4 39(3A/B):L262 [6]Sheibani H,Liu Y C.Sakai S.et al.The effect of applied magnetic field on the growth mechanisms of liquid phase elec troepitaxy.Int J Eng Sel,2003.41(3/5):401 [7]Bird G A.Molecular Gas Dynamics.Oxford:Clarendon 20 2224262830 Pres5,1976 初始氧体积分数% [8]Baganoff D,McDonald J D.A collision-section rule for a parti- 图8氧的体积分数变化量随初始氯体积分数的变化情况 cle simulation method suited to vector computers.Phys Flulds A,1990,2(7):1248 Flg.8 Change of oxygen concentration increment with initial [9]沈青.稀薄气体动力学，北京：国防工业出版社，2003 oxygen concentration [10]Vincenti W G.Kruger C H.Introduction to physical gas dy- namics.New York:Wiley,1965 4 结论 [11】周斌.内燃机富氧节能的研究一富氧空气的获取.西南 以氧气和氮气的混合气体为例，运用蒙特卡 交通大学学报，1995,30(3)：312 罗直接模拟方法，对顺抗磁性混合气体在高梯度

V o l 。 2 8 N 0 . 5 蔡军等 : 顺抗磁性 混合气体在梯度磁场 中流动的 , 特卡 罗 I 接模拟 说 压力的变化对 氧的体积分数变化量的影响不是 很大 . 虽然 压 力增大会 导致氧气体积磁 化率的增 加 , 但是 其分子数密度也 同时在增加 , 因此具体到 每个氧分子上 的磁力并没 有发生 变化 . 所以磁 场 间隙内氧的体积分数变化量也 就不 会发生变化 . 一 . — 一一\ 一 . / . 月呀 .0 喇契耸撅尔彩名屏芝 0 3 占 0 00 4 0() 0 表帐 力 , P/ aP 6 0 0 0 图 , 级的体积 分数变 化 t 随压力的变化情况 F ig . 7 C h a n ge o f o x y ge n e o n e e o t ar t i o n lcn 代 m e l t w lt h p r e, s吮 图 8 所示 为磁 场间隙内氧气的体积分数变化 量随初 始氧体积分数 的变化情况 , 其参数条 件为 尸 = 1 0 1 3 2 5 P a . T 一 2 3 9 K . 粤华 = l o o T , . m 一 1 . 一 v 一 “ ~ “ 一 , ` ~ “ “ ` ’ d X ` “ ” ` Q l = 40 0 m .L m in 一 ’ . 由图 8 可知 , 随着氧气初 始 体积分数的增加 , 磁 场间隙内氧的体积 分数变 化 量越来越大 . 当氧气的初 始体积分数由 2 0 % 增加 到 3 0 % 时 , 氧 的体积分数变化量 由 0 . 12 % 增加到 0 . 59 % . 氧气初 始体积 分数越 大时 , 氧分子 的体 积 磁 化率也 就越 大 , 由式 ( l) 可知它 所受到 的磁力 就越 强 , 因此 留在磁场 间隙内的氧气分子 个数就 越 多 . 磁场中的流动进 行了计算机 模拟 , 并把模拟 结果 与已有的实验结 果进 行 了对 比和 分析 , 得 出 如下 结论 : ( 1) 相 同条件下 , 模拟 结果与文 献实验 结果 吻合较好 . ( 2) 氧的体积 分数增量 与磁 场强度和 梯度的 乘积 的绝对 值成正 比 , 当毕 由 1。。 护 . m 一 , ` , 增加 、 ’ ` ’ ” J 一 ~ ’ ` ~ ~ 一 ~ ’ 一 d X ~ ` ” ” “ ~ 一 曰 ~ 到 1 0 0 0 T 2 · m 一 ’ 时 , 氧的体积 分数增量 由 0 . 12 % 增加到 1 . 8 9 % . ( 3) 氧的体积分数增量与混 合气体温度成反 比 , 当温 度由 2 7 3 K 升高到 3 3 K 时 , 氧的体积 分 数增 量 由 0 . 4 % 减小到 0 . 1% . ( 4 ) 压力 的变化对氧的体积 分数增量的影 响 不是 很大 . ( 5) 氧的体积分数增量与混 合气体 中氧的初 始体积分数成正 比 , 当混 合气体中氧气的初 始体 积分数由 2 0 % 增加到 3 0 % 时 , 氧的体积分 数增 量 由 0 . 1 2 % 增 加 到 0 . 5 9 % . 今 考 文 献 20 22 2 4 2 6 2 8 3 0 初 始氧体积分数 /% 圈 8 级 的体积 分数 变化. 随初始级体积分橄的变化情况 F lg . 8 Ch a gn e o f o x y ge n co n 沈n tar t l o n I n e er 吹 n t w i t h i n i t i a l o x y ge n o o n c e l t , t i o n 4 结论 以氧气和 氮气的混 合气体为例 , 运 用蒙特 卡 罗直接模拟 方 法 , 对 顺抗磁性 混 合气体在高梯度 【l] 朱重 光 . 横 向非均匀磁场对原子的作用 . 淮北煤师院学报 : 自然版 . 1 99 6 . 17 ( l ) : 30 〔2] 朱重光 . 磁分 制氧机理 . 淮北煤师院学 报 , 19 97 , 18 ( 1 ) : 2 2 [ 3 ] U e on S , lw a as k a M , Eg u e hi H , e t a l D y n a m i e s be h a v i o r o f g as fl ow i n g r a d i e n t m a , e t i e f i e ld s . I E E E T, 朋 M . g l , 1 99 3 , 2 9 ( 6 ) : 3 2 6 4 [ 4 ] Wa k a ” m a N I , S 呀i e M . M吧 n e t i e p or om t i o n o f co m b us t io n i n d if f u s i o n fl a m e s . P h y s B , 19 9 6 , 2 1 6 ( 3 一 4 ) : 4 0 3 [ 5 ] W a k a y a m a N l , w a k a y a am M . M 吧 n e t i 。 。 e e e l e ar t ion o f i n · h a l e d a n d e x h a l e d n o w , i n ber a t h i n g . J p n J A P PI P h y s , 2 0 0 0 . 39 ( 3 A / B ) : L 2 62 [ 6 ] S h e ib a n i H , I 碑 i u Y C , aS k a i S , e t a l . T h e e ff ce t o f a p p lide m a g n e t i e f i e ld o n t h e g or w t h m ce h a n i s m s o f Iiq u id p h a s e e lce - t ore p i t a x y . I nt J E雌 cS i , 2 0 0 3 , 4 1 ( 3 / 5 ) : 4 0 1 〔7 ] Bi记 G A . M o lec u l a : G a s D y n a m i e s . O x fo r d : C l a r e n d o n P r e s , 19 7 6 [ 8 〕 a3I g a n o f f D , M e oD n a ld J D . A co lli s i o n 一 s e e t i o n r u l e of r a 钾rt i - e l e s im u l a t i o n m e t h浏 s u i t e d t o v ce t o r e o m P u t e rs . P h y s n 川ds A , 19 9 0 , 2 ( 7 ) : 12 4 8 【9] 沈青 稀薄气体动力学 . 北京 : 国防工业 出版社 , 20 03 【10 ] V i n e e n t i w G , K r u g e r C H . I n t rod u e t i o n t o p h y s i e a l g a s d y - n a m i e s . N e w Y o r k : W il e y , 19 6 5 【1 1] 周斌 . 内燃机富 氧节能 的研究 — 富氧空气的获取 . 西南 交通大学学报 , 19 9 5 . 3 0 ( 3 ) : 3 12 4 nU ó , 训脚概求彩长名屏碧