介绍
热喷涂和静电喷涂具有喷涂材料利用率高、施工速度快、安全等优点，近年来在各国得到广泛推广。喷涂技术在工程领域中发挥着至关重要的作用，并吸引了大量基于计算流体动力学（CFD）的研究，包括粉末传输行为、1-3粉末浓度分布4、5和喷涂过程的研究。6 , 7获取颗粒分布数据对于工件加工具有重要意义。此外，工作环境中空气颗粒物的传播也威胁着人类健康。因此，了解工作环境中空气颗粒物的扩散情况对于制定工业防护策略具有重要意义。CFD 已广泛应用于模拟空气中颗粒的传输和扩散。4 , 5 , 8然而，传统CFD方法计算成本较高，限制了其通过风扇旋转引起的无组织气流扰动来模拟高浓度条件下喷雾颗粒分布的应用。
近年来，人们对风扇旋转引起的气流扰动进行了研究，重点关注颗粒-风扇相互作用9 – 11和气流模式12 – 15然而，风扇引起的气流在不同的吹风方向下具有不同的流型。迄今为止，在房间尺度上捕获那些复杂的气流场仍然具有挑战性。一是吊扇和落地扇诱导的下吹和上吹气流不均匀，其核心区气流较强，周边区气流较弱，导致喷涂环境出现无组织的气流扰动。 。此外，高浓度可能会产生沉重的CFD计算负担。需要基于流体流动控制方程进行合理简化，以降低CFD模型的计算复杂度。
本研究旨在研究风扇旋转时无组织气流扰动的变化，并模拟高浓度条件下喷雾的颗粒分布。模型喷涂车间内的气流由喷涂气流建立并由风机再循环。理论计算、实验和CFD是研究密闭空间气流的主要方法。理论计算基于流体流动控制方程和斯托克斯数，将高浓度颗粒的分布计算简化为气体组分的分布状态。采用MRF模型、Fan模型建立CFD模型以加快仿真速度，并通过现场试验进行验证。据我们所知，本研究首次利用CFD模型预测高浓度条件和无组织气流扰动下的喷洒颗粒物浓度分布。相信这种Stokes-MRF-Fan-CFD模型将在流体力学领域得到广泛的应用。
方法
数值方法
在有气流的封闭空间（喷涂车间）中，本研究主要关注的是气流运动。将空气视为单一连续物质，不考虑空气的成分以及各成分分子水平的分散性。经计算（见2.4节），喷口处气流最大速度不大于71 m/s，声速为340 m/s。气流速度约为0.2Ma。在流体力学中，当流量小于0.3Ma时，可以不考虑流体的压缩性。空气可以被视为不可压缩的流体。16空气密度与温度的关系符合理想气体状态方程。并且喷涂车间配有恒温恒湿设备，无需计算温度场。
流量模型
CFD模拟分析离不开流体力学的三个守恒方程（流体流动的控制方程）17：连续性方程（方程（1））、动量方程（方程（2））和能量方程（方程（7）），其中对应三个基本物理原理：质量守恒定律、牛顿第二定律（动量守恒定律）、能量守恒定律，如图S1流体力学基本控制方程的补充材料。通过求解质量守恒、动量守恒、能量守恒、化学成分浓度守恒等一系列偏微分方程，可以实现对流体流动、热交换、物质输运等物理现象的分析。

动量方程是矢量方程，连续性方程和能量方程是标量方程(包括三个速度方向的方程)。包括流体流动控制方程、压力和能量（内能）计算公式等7个方程。方程包括密度、三个速度方向的速度分量、内能、温度和压力等七个未知量。上述方程是封闭的。在这项研究中，密度和温度是常数，其他未知量可以求解。
在高浓度条件下，粒子碰撞后运动轨迹会发生变化。为了减少计算量，粒子之间的参数设置为碰撞后不合并。碰撞后粒子的后续运动可以根据动量方程计算。
拖动模型
涂层材料（环氧树脂）直径大于10μm，不会像气溶胶一样长期漂浮在空气中。20当环氧树脂接触到车间的表面（例如天花板、地板和墙壁）时，它会粘附在表面上。并且环氧树脂不再随气流飘入车间。在计算颗粒物浓度分布的过程中，需要考虑颗粒物重力引起的沉降效应。对于颗粒的沉降，颗粒之间的阻力是阻力系数的线性函数。

移动参考系
欧拉网格通常用于流体动力学计算。22欧拉网格可以处理大变形（因为节点不运动），23但处理节点运动的能力不足。24然而，现实中网格边界运动的例子有很多，例如发动机气缸运动、阀门开闭、机翼运动等，以及本研究中喷涂车间风扇叶片的旋转。需要将网格运动引入CFD计算中来模拟风扇叶片的旋转。
CFD计算的基本物理量是速度、温度、压力和成分，不计算网格节点的位移。25因此，为了使网格移动，通常对节点施加的物理约束是速度。网格的移动速度是相对于参考坐标系定义的，参考坐标系可以是静态的（实验室坐标系，本研究中为喷涂车间的参考系），也可以是相对于参考坐标系旋转和平移。在不需要精确时间解的稳态仿真或瞬态仿真中，移动参考系 (MRF) 提供了一种将旋转和平移建模为稳态问题的方法，26同时保持网格静止。图1(a)显示了考虑恒速旋转和平移的运动参考坐标系。

革命模型
喷涂车间的搅拌风机不仅有自转运动，还有公转运动（摆动）。风扇的旋转运动方式和参数可以在建模时直接添加，但公转运动比较复杂，需要单独定义。喷涂车间的三台风机绕着旋转轴旋转，需要建立它们的局部坐标系。对于简单的 CFD 模拟应用，可以使用实验室坐标系。但喷涂车间内存在多台风机的自转和公转运动，气流运动复杂。有必要定义和使用局部坐标系。
实验室坐标系是定义其他局部坐标系的基准，实验室坐标系属于笛卡尔坐标系。局部坐标系可以定义为笛卡尔坐标系、柱坐标系或球坐标系，并且可以使用基本实验室坐标系或先前定义的笛卡尔局部坐标系作为参考系。27本研究定义了三个笛卡尔局部坐标系，如图2所示。风扇的转数可以通过以三个风扇的笛卡尔局部坐标系为参考系定义旋转运动来定义。

在实验流体动力学中，斯托克斯数是粒子图像测速 (PIV) 实验中流动示踪剂保真度的度量，其中非常小的颗粒被夹带在湍流中并通过光学观察来确定流体运动的速度和方向（也称为速度）流体场）。为了获得可接受的跟踪精度，粒子响应时间应快于流的最小时间尺度。斯托克斯数越小表示追踪精度越高。对于Stk  >> 1，颗粒将从流中分离，尤其是在流突然减速的情况下。对于Stk  << 1，颗粒紧密跟随流体流线。如果Stk  < 0.1，则跟踪精度误差低于 1%。
罗马诺等人。34研究了PIV实验中不同粒径的丙二醇颗粒（蓝点）在驻点气流场（灰色流线）中的平流运动。图 3显示了两种不同粒径的 PIV 跟踪精度的比较。模拟丙二醇颗粒（蓝点）在驻点流场（灰色流线）中平流。请注意，1 毫米的颗粒撞击到停滞板上，而 0.1 毫米的颗粒则沿着流线运动。在本研究中，环氧树脂粉末的粒径仅为~30μm。介质为空气，丙二醇的密度为1.038g·cm -3，环氧树脂粉末的密度为1.204g·cm -3。斯托克环氧树脂 ≈0.104 Stk丙二醇可由式（13）和（14 ）求得。可以认为环氧树脂粉末颗粒在空气中跟随流体流线运动。

本研究的重点是确定高浓度条件和无组织气流扰动（即气固两相流问题）下喷雾的颗粒分布。由于本研究中颗粒的体积分数较小（<< 1%，仅为 1/8.7564 × 10 −3，计算过程见2.4节），我们根据本节的讨论结果，采用合理的简化来求解问题。环氧树脂粉末颗粒跟随气流，并且注入的空气和环氧树脂粉末的混合物被设定为气体成分。该气体成分与喷涂车间内的原始气体一起随气流流动，模拟该气体成分在车间内的分布。由此，根据气体成分的分布状态来评价环氧树脂粉末的分布，并数值计算喷涂车间内环氧树脂粉末的浓度分布。
喷涂车间模型
喷涂车间3D模型是根据实际喷涂车间布局进行一定简化后建立的，包括一个喷口和三个风扇。喷涂车间由一个密封室、三台风机和颗粒浓度监测系统组成。表1列出了喷涂车间的尺寸，配合图2可以帮助构建喷涂车间的几何印象。车间长（L）、宽（W）、高（H）为3690毫米×1490毫米×3000毫米。该车间容积为16.49 m 3。车间表面（如天花板、地板、墙壁）均由不锈钢板制成。总共使用了三个编号为 1-3 的风扇，分别位于天花板下方和地板上方。它们被放置在车间的中心线（长度方向）上，形成气流再循环区，也是喷涂区。另外，16个操作孔分布在距喷涂车间底部1.1m处。

喷涂车间内有8个喷涂操作通道对工件进行操作。每个通道的中心线距两侧墙体和相邻通道中心线的距离为145、400、600、400、600、400、600、400、145毫米。操作通道中心线高度为1100mm。喷涂车间外部铺有保温材料。两端厚度均为305毫米。喷涂车间长度方向的一端设有用于喷涂涂料的喷口。喷嘴直径为12毫米，高度为1100毫米，与操作通道中心线的高度一致。图2为喷涂车间操作通道尺寸及喷口位置。图S2补充材料中是拓扑结构和流体域中操作通道的细节。
本研究选择了三台直径为31厘米、带有四个金属叶片的吊扇和落地扇（型号：KDS-310），固定在车间的中心线（长度方向）。风扇叶片的几何尺寸如图4(a)所示。风扇配备交直流有刷电机（额定工作电压为48V），能够以每分钟1450转（r/min）的速度进行不同方向旋转，转速为6.11r/min。


为了提高计算精度并防止计算发散，需要设置非常小的时间步长，这需要非常大的计算量。因此，Fan模型应进行简化。在对风机模型进行简化之前，需要通过MRF模型计算真实风机的流量特性，得到特定转速下风机流量与压力的关系。本研究对风扇进行建模并将其置于流体通道中，对真实风扇进行流体流量计算，得到转速为1450 r/min（正常转速）时叶片上下游压差风扇的工作状况）。图4(b)是根据KDS-310风机叶片数据建立的几何模型。风扇模型的叶片倾斜角为16.7°，轴承振动幅度为0.15 mm。
当能够获得风扇性能曲线并且风扇的精确几何形状对计算结果不产生影响时，可以使用简化风扇模型来简化风扇模型。KDS-310风机属于轴流风机。气流的流入和流出方向相同，即沿风机轴向流动。本研究利用流体软件提供的风扇动量源模型建立轴流风扇的风扇模型。将CFD计算出的1450 r/min转速下叶片上下游压差引入Fan模型作为风扇的动量源，可得到将风扇流量和压力结合起来的集总Fan模型已确立的。风扇被建模为一个体积区域，当气流经过该区域时，会引起风机上下游压力变化。计算时风扇不动，可以大大减少计算量。此外，Fan模型还可以模拟风扇引入的涡流，并计算风扇下游的涡流速度。
网格生成
使用多面体网格对喷涂车间模型的体积进行离散化。喷涂车间内有流线方向快速变化的扇区。为了保证计算精度，对包围风机转数范围的区域进行网格加密35，并通过重叠网格对风机区域和喷涂车间区域进行网格重叠处理。36
考虑到模型的计算速度，设置了风机区域和喷涂车间区域。包括风机的真实性能计算和喷涂车间的流量计算。对于本研究中使用的高雷诺数模型（k -Epsilon模型）， y +的取值应满足30≤y +  ≤300的要求 ， 一般研究中取y + =30。37 , 38边界层厚度和层数分别为2mm和5，y +的值为30。扇面网格数为113万，体积网格数为500万。喷涂车间表面网格数60万个，体网格数323万个。由于风扇区域附近的气流较为复杂，因此风扇区域的网格比其他区域更密集。喷涂车间最终网格如图5所示。

图5 . (a) 扇形网格细化与 (b 和 c) 周围网格尺寸之间的比较。喷涂车间的(d)表面网格和(e)体积网格。
模型验证
为了确保 CFD 模型的准确性，对颗粒浓度进行了验证。由于高浓度下的颗粒分布没有类似的验证数据，因此我们使用实测数据的实验数据集。喷涂车间内的气流是复杂的，它随着喷管和风机的运行条件而变化，在喷涂过程中不断变化。因此，选择喷淋和沉降的整个周期进行验证。
测试设置包括三个阶段，例如喷雾、分散和沉降。设置如下。1-15秒：将涂料以0.7MPa吹入喷涂车间。1-30秒：使用风扇使涂料分散均匀。31-90秒：静置以稳定包衣材料的浓度。喷涂车间涂料平均浓度为1000mg·m -3。表2显示了涂层材料的物理参数。

将空气压缩机的参数（储气罐容积V  =160L、初始压力p  =0.7MPa、终了压力p  =0.4MPa、室温T  =294.15K）带入理想气体方程pV  =  nRT。可以得到从空压机储气罐注入的气体总量Δn  =19.64mol，体积V  =120L。注入粒子的体积比为1/8.7564×10 -3。喷口半径r  = 6 mm，喷射时间t  = 15 s，喷射空气速度v  =  V /(πr 2  ×  t ) = 70.74 m/s ≈ 0.2 < 0.3 Ma，因此空气可视为不可压缩流体。16
颗粒物浓度监测
如今，颗粒物浓度监测的标准已经以生活环境为基础。这些标准的一个共同特点是检出限低、灵敏度高（颗粒物浓度检测标准补充材料表S1 ）。然而，该颗粒在喷涂车间空气中浓度较高（>1000 mg·m -3）。在线获取这些信息一直具有挑战性，因为它很容易超过传统颗粒浓度监测仪的最高检测限。
静电耦合和数字信号处理的最新进展使得实时监测高浓度车间颗粒物成为可能。本研究设置KJ-Z型颗粒物浓度监测仪（华瑞环境科技有限公司）对喷涂车间颗粒物浓度进行监测。KJ-Z的原理是监测颗粒物接触传感器时产生的交流感应电动势信号。感应电动势与仪器范围内的颗粒浓度成正比。它通过颗粒的感应电动势信号来区分颗粒浓度，从而对浓度信号的强度和数量进行分类。该方法测量颗粒物浓度上限较高，但时间分辨率较低。
KJ-Z的进风口位于喷涂车间的操作孔处，进风口方向与喷涂气流方向垂直。进风口距墙壁的距离为0.745 m。进风口周围没有屏蔽。
结果与讨论
气流模式
风机计算结果
风扇性能计算为稳态计算，采用MRF模型模拟风扇旋转。对MRF模型计算出的风机结果文件进行处理后，得到叶片上下游的压差。该压力差可用于测量风机动量源的特性。当风扇被视为Fan模型时，认为通过Fan模型的气流受到上游和下游压差的加压。根据计算，风机上游静压为-16.01 Pa，风机下游静压为0.88 Pa，风机引起的压力上升为0.88−(−16.01) = 16.89 Pa。根据制造商提供的性能参数，KDS-310风机正常工作条件下上下游压差为16Pa（25℃）。MRF模型在压力上升方面的结果显示，绝对误差为-0.89 Pa，相对误差为-5.56%，这表明MRF模型的模拟结果总体良好。图6（a）是风扇的静压及其上游和下游的轮廓线。图6(b)是风扇的速度云图。图6(c)是风扇的流线图。

图6 . (a) 上下游静压及其等值线，(b) 速度等值线，(c) 风扇流线。
喷涂车间计算结果
图7给出了喷射车间喷射、混合扩散和沉降三个时期的流线分布和速度矢量场。可以看出，在喷射工况下，由于涂料颗粒的喷射速度较高，扰动了风机流线，使其偏离垂直方向，形成集中在喷涂车间中部的三个涡流。使喷涂车间内的颗粒物无法有效分散。在混合扩散工况下，喷射工况下形成的涡流被风扇气流破坏。新形成的大型涡流位于喷涂车间上部，能有效分散颗粒而不沉降。沉降工况下，风速逐渐减小（<1 m·s−1 )，涡流逐渐稳定。颗粒开始沉降并积聚在远离喷嘴且未被涡流覆盖的角落区域。60 s后，风速稳定在0.1 m·s -1左右，涡流尺寸变大并逐渐覆盖整个喷涂车间。沉降在角落区域的部分颗粒开始随流线向喷涂车间上部上浮，此时喷涂车间内颗粒浓度趋于稳定。可以观察到，三种工况不同气流组织方式下，喷涂车间气流组织差异主要集中在喷涂车间中部区域、风机区域和远离喷口的角落区域。吊扇和落地扇产生的向下和向上吹送的气流不均匀，其核心区域产生较强的空气运动，但周边区域产生较弱的气流，导致喷涂环境中出现无组织的气流扰动。

图 7 . 喷涂车间流线分布及速度矢量场：(a) t  = 8 s, (b) t  = 23 s, (c) t  = 30 s, (d) t  = 40 s, (e) t  = 50 s ，(f) t  = 60 秒，(g) t  = 70 秒，(h) t  = 80 秒。
粒子的分散和命运
1-15 秒内喷雾分散颗粒
1~15s后，喷涂车间风机开启，利用空压机将涂料注入喷涂车间。根据上述讨论结果，本研究简化了模型中涂层材料的以下特性。注入了涂料颗粒的气体被视为气相。图8为喷涂车间1~15 s内的速度云图和颗粒质量分数。可以看出，由于喷射风速（71 m·s -1）远高于风机风速（4 m·s -1），颗粒物集中在中间的高浓度区域。喷涂车间内并没有有效分散在喷涂车间内。

图 8 . 喷涂车间1~15 s内的速度云图和颗粒质量分数：(a) t  = 2.5 s，(b) t  = 5 s，(c) t  = 7.5 s，(d) t  = 10 s，(e ) t  = 12.5 秒，(f) t  = 15 秒。
16-30秒内风扇气流使颗粒分散
16-30 秒后，停止注入干扰剂，用风扇搅拌，使颗粒分散均匀。图9为喷涂车间16~30  s内的速度云图和颗粒质量分数。可以看出，停止释放颗粒后，在风机的搅拌作用下，颗粒在喷涂车间内得到有效分散，没有沉淀。

图 9 . 喷涂车间16~30  s内的速度云图和颗粒质量分数：(a) t  = 17.5 s，(b) t  = 20 s，(c) t  = 22.5 s，(d) t  = 25 s，(e ) t  = 27.5 秒，(f) t  = 30 秒。
31-60 秒内颗粒沉降
31-90 秒，喷涂车间风机关闭，开始工件作业。其中，车间静置31-60 秒，以稳定颗粒物浓度。此时，颗粒的喷射和风扇搅动的气流扰动已经停止，颗粒开始靠重力沉降。图10为喷涂车间31~60  s内的速度云图和颗粒质量分数。可以看出，风机关闭后，喷雾车间内颗粒沉降，而携带气体的颗粒在远离喷口的区域开始沉降和积累。

图 10 . 喷涂车间31~60  s内的速度云图和颗粒质量分数：(a) t  = 32.5 s，(b) t  = 37.5 s，(c) t  = 42.5 s，(d) t  = 47.5 s，(e ) t  = 52.5 秒，(f) t  = 57.5 秒。
61-90 秒内颗粒沉降
61~90 s期间 ，喷涂车间的工况与31~60 s期间一致 。颗粒的喷射和风扇引起的气流扰动已经停止，颗粒靠重力沉降。图 11为喷涂车间 61~90  s 内的速度云图和颗粒质量分数。可以看出，70  s前颗粒就沉降在喷涂车间的下部区域。70s后 ，携带颗粒的气流开始向喷涂车间上部飘浮，说明喷涂车间内仍有风速。

图 11 . 喷涂车间61~90  s内的速度云图和颗粒质量分数：(a) t  = 62.5 s，(b) t  = 67.5 s，(c) t  = 72.5 s，(d) t  = 77.5 s，(e ) t  = 82.5 秒，(f) t  = 87.5 秒。
颗粒分散分析
图12为喷涂车间颗粒物浓度-t模拟曲线。图 13显示了运行通道中颗粒物浓度的模拟值。从图中可以直观地看到颗粒的喷洒、分散和沉降过程。

图 12 . 喷涂车间中心操作通道颗粒浓度-t模拟曲线。

图 13 . (a)第七和(b)第八操作通道中颗粒浓度的模拟值。
喷射颗粒时，喷射的颗粒集中分布在喷涂车间的中部。这反映在图12和图13 中0-15秒颗粒浓度的突然跳跃，且集中在喷涂车间的5~8个通道。颗粒喷射后，携带颗粒的空气在惯性作用下继续向前运动，从而在颗粒的喷射路径上形成低压区域。此时，受到搅拌风扇扰动的自由气流将充满低压区，导致颗粒浓度下降。最后，颗粒物的喷射和气流扰动停止，颗粒物浓度趋于稳定。在3-8个通道中，前20秒内颗粒浓度是标准工作浓度的3倍。工作人员应做好呼吸和身体防护，加强通风。
根据3.2.4节的结果，在70 s之前，颗粒物沉降在喷涂车间的下部区域。70s后，由于喷漆车间内的流动，携带颗粒的气流开始飘向喷漆车间上部。工作人员应采取呼吸和眼睛防护措施。图 13显示了第 7 和第 8 个操作通道在 31-60 和 61-90 s 内颗粒物浓度的标准差（σ  ） 。61~90 s的颗粒浓度模拟值的分散程度低于31~60 s，这表明喷涂车间颗粒的分散更加稳定，适合工件操作。
定量验证
通道7和通道8是喷涂工作中最常用的位置，因此我们选择这两个通道进行数值模拟验证。图 14为运行通道颗粒物浓度实验值与模拟值对比。图14中，实验值是真实观测颗粒浓度数据，模拟值是通过CFD模型模拟的。计算出仿真值并与实验值进行比较，评价模型的可用性。通过现场试验验证了模型的仿真结果。可以看出，图 14中两条曲线几乎彼此重叠。

图 14 . (a)第七和(b)第八操作通道中颗粒浓度的实验值和模拟值比较。
图15显示了模拟的操作通道中颗粒浓度的绝对误差率和相对误差率。对模拟结果误差分析表明，绝对误差和相对误差大部分集中在±200 mg·m -3和±20%区间内，模拟结果总体较好。由于研究过程中模拟值会受到很多方面的影响，图15中的几个相对误差较大。

图 15 . 模拟（a）第七和（b）第八操作通道中颗粒浓度的绝对误差率和相对误差率。
本研究还计算了 CFD 模型的准确性。20模拟精度计算公式如下：
 
其中，S i表示模拟值，E i表示实验值。CFD模型模拟结果的准确度等级为88.5747。成绩表明，该CFD模型能够有效模拟无组织气流下高浓度喷涂车间的颗粒浓度。从绝对误差、相对误差和模拟精度等级来看，总体结果较好，误差率在可接受的范围内，表明该CFD模型是本研究的理想模型。初步解决了高浓度工况、无组织气流扰动下CFD预测喷洒颗粒物浓度分布计算量大的问题。
结论
由于密闭空间内气流无组织、颗粒物浓度高，传统的平均浓度计算方法无法真实反映颗粒物的浓度分布。Stokes-MRF-Fan-CFD方法的应用是模拟密闭空间颗粒浓度分布的一次必要的创新尝试。本研究的主要结论如下。
(1)根据流体流动控制方程、Stokes数、MRF模型和Fan模型对问题进行合理化简求解。从而将高浓度粒子的分布计算简化为气体成分的分布状态。喷口处携带颗粒的气流速度约为0.2Ma，可视为不可压缩流体。可以认为，环氧树脂粉末颗粒（Stk环氧树脂 <0.1）在空气中跟随流体流线运动，注入的空气与环氧树脂粉末的混合物被设定为气体成分。
(2)颗粒物集中在喷涂车间中部高浓度区域，1~15s内未有效分散。在风机的搅拌作用下，颗粒在16-30秒内有效分散，不产生沉淀。风机关闭后，颗粒在喷涂车间内沉降，31~60s内，携带气体的颗粒开始在远离喷口的区域沉降和积累。61~90 s，携带颗粒的气流开始向喷涂车间上部飘浮，说明喷涂车间内仍有风速。61~90 s的颗粒浓度模拟值的分散程度低于31~60 s，这表明喷涂车间颗粒的分散更加稳定，适合工件操作。
(3)CFD模型的应用取得了良好的效果。绝对误差和相对误差大部分集中在±200 mg·m -3和±20%区间，模拟精度为88.5747，统计结果表明总体结果较好，误差率在可接受的范围内。