流体机械的内流仿真在现代透平机械设计中起着至关重要的作用，非均匀介质复杂流动现象（如相变和空化）为传统数值模拟方法带来挑战. 传统网格法（如有限体积法）在计算叶轮机械均质流动方面已较为成熟且有广泛工业应用，但前处理过程须分模块划分网格在计算时要引入混合平面法、旋转坐标系法或转子冻结法等假设或简化方法，捕捉多相界面须基于体积分数进行界面重构，上述设置易引入误差，丢失部分流动细节.

无网格粒子法在拉格朗日框架下，采用无固定拓扑关系的流体粒子代替网格和节点，适用于计算非定常大变形流动. 采用整体离散和整体求解，无须将流体域分块建模和拼接、无须额外处理动静交界，易于捕捉多相界面运动，有望为叶轮机械复杂内流模拟，尤其是针对非均质多相流动，提供一种新的解决方案.

已有许多学者将无网格法应用于流体机械内流数值模拟，Sun等^[1]用移动粒子半隐式法(moving particle semi-implicit，MPS)模拟带有6个直叶片与6个导叶的三维搅拌器. 王锋等^[2]利用MPS方法研究Y型微混合器内的流动混合机理. Rahim^[3]采用改进的不可压缩光滑粒子流体动力学方法(incompressible smoothed particle hydrodynamic，ISPH)模拟微泵混合器内的流动和传质. Min等^[4]建立了一个高压斜齿轮泵模型，采用MPS方法获得了较FVM更优的结果. Kakuda等^[5]用MPS方法研究液环泵内水环的形成. 此外，由于在齿轮交界处难以建立高质量的动网格，无网格法也常用于模拟齿轮箱或减速器内的油液流动^[6-8].

以上研究针对的几何模型大多较为规则且不存在尖角或薄壁结构，使用的壁面边界均是布置三层壁面粒子^[9-10]并使内壁面粒子参与压力计算. 对于存在复杂回转面、尖角以及细长弯扭叶片的流体机械，如离心泵，较难采用传统方法直接布置. 布置精度不足引起的粒子数密度缺失也会导致粒子穿透.

基于距离函数的单层壁面模型^[11]常用于离散复杂几何壁面，但直接引入距离函数容易导致近壁面液体粒子压力剧烈波动^[12]，进而影响计算的稳定性. 此外，壁面函数在三维情况下较难准确补偿粒子数密度缺失，Zhang等^[13]提出基于几何关系的边界粒子自适应布置法，但参与压力计算的边界粒子仍增加了计算量.

本研究基于MPS法，针对流体机械内复杂壁面离散问题，结合光滑壁面边界模型^[14]与虚拟粒子法^[15]，建立了一种复杂曲面离散方法. 该方法中壁面粒子仅有单层且不参与压力计算，边界截断通过虚拟粒子进行补偿，壁面作用力由接触模型来提供. 利用局部加密来提高复杂型面的离散精度，定义法向光滑角度来描述模型离散精度. 通过计算2个三维复杂静水算例以及诱导轮水中旋转算例，验证该离散建模方法的有效性. 本研究所提模型是无网格法应用于流体机械内流模拟的重要前处理工作之一.

MPS方法用于模拟不可压缩流动，质量守恒方程与动量守恒方程如下：

(1) $ \nabla \cdot {{{{\boldsymbol{u}}}}} = 0, $

(2) $ \rho \frac{{{\text{d}}{{\boldsymbol{u}}}}}{{{\text{d}}t}} = - \nabla p+\mu {\nabla ^2}{{\boldsymbol{u}}}+\rho {{\boldsymbol{f}}}. $

式中： $ \rho $为流体密度； $ {\boldsymbol{u}} $为流体速度； $ p $为压力； $ \mu $为动力黏度； ${\boldsymbol{ f}} $为体积力.

MPS法将计算区域离散为带有物性、流动以及热力学参数的粒子，通过粒子间相互作用，离散控制方程中的各项微分算子，导出梯度算子和拉普拉斯算子的光滑近似式，计算各时间层粒子的流动参数，随时间推进动态获取整场流动信息.

粒子间作用力的强弱可以用核函数描述，本研究采用立方样条核函数. 流体机械流动模拟中通常存在较强的动静相干现象，压力波动剧烈，粒子常受到较大的挤压和瞬时冲击. 立方样条核函数曲线相对光滑，两端变化趋势缓慢且存在有限值，当两粒子发生挤压而导致粒子间距急剧缩小时，不容易出现瞬时核函数值过大而导致压力振荡，适用于计算剧烈流动问题. 该核函数已经被应用于水轮机模拟^[14]和齿轮箱内油液流动^[16]仿真并经过验证. 核函数表达式如下：

(3) $ w(R)=\alpha \times \left\{\begin{array}{l}\dfrac{2}{3}-{R}^{2}+0.5{R}^{3}\text{，}\text{   }R\leqslant 1\text{；}\\ \dfrac{1}{6}{\left(2-R\right)}^{3}\text{，}\text{       }1 < R\leqslant 2\text{；}\\ \text{ }\stackrel{}{{\displaystyle }}\text{  0}\text{，}\text{               }R > 2.\end{array}\right. $

(4) $ R = {r \mathord{\left/ {\vphantom {r h}} \right. } h}. $

式中： $w(R)$为核函数值； $ r $为两粒子 $i$、 $j$之间的距离， $ r = \left| {{{{\boldsymbol{r}}}_i} - {{{\boldsymbol{r}}}_{j}}} \right| $； $h$为光滑长度； $\alpha $为归一化系数.

梯度算子与拉普拉斯算子可以通过空间泰勒展开结合光滑近似公式求得^[17]. 为了保证动量守恒，采用改进后梯度公式^[18]. 梯度算子与拉普拉斯算子计算公式如下：

(5) $ {\left\langle {\nabla \varphi } \right\rangle _i} = \frac{d}{{{n^0}}}\sum\limits_{j \ne i} {\frac{{{\varphi _j}+{\varphi _i}}}{{{{\left| {{{{\boldsymbol{r}}}_j} - {{{\boldsymbol{r}}}_i}} \right|}^2}}}} ({{{\boldsymbol{r}}}_j} - {{{\boldsymbol{r}}}_i})w(\left| {{{{\boldsymbol{r}}}_j} - {{{\boldsymbol{r}}}_i}} \right|), $

(6) $ {\left\langle {{\nabla ^2}\varphi } \right\rangle _i} = \frac{{2d}}{{\lambda {n^0}}}\sum\limits_{j \ne i} {({\varphi _j} - {\varphi _i})} w(\left| {{{{\boldsymbol{r}}}_j} - {{{\boldsymbol{r}}}_i}} \right|). $

式中： $\varphi $为任意标量函数， $d$为维度， ${n^0}$为粒子数密度常数， $\lambda $为修正系数.

MPS法基于粒子数密度的偏移与修正来构造不可压缩模型，粒子数密度定义如下：

(7) $ {\left\langle n \right\rangle _i} = \sum\limits_{j \ne i} {w(\left| {{{{\boldsymbol{r}}}_j} - {{{\boldsymbol{r}}}_i}} \right|)} . $

在拉格朗日框架下，将1个时间步分为显、隐2步来解耦合速度与压力. 显式计算黏性力以及体积力，得到粒子位移、速度的估算值，由此估算粒子数密度 ${n^*}$，利用该估算值结合质量守恒、动量守恒方程推导出压力泊松方程，由此方程隐式求解压力，求得压力梯度力并修正各项参数.

采取散度自由条件法^[19]构造压力泊松方程，引入伪可压系数 $\gamma $，压力泊松方程如下：

(8) $ {\left\langle {{\nabla ^2}{p^{k+1}}} \right\rangle _i} = (1 - \gamma )\frac{\rho }{{\Delta t}}\nabla \; \cdot \;{{{\boldsymbol{u}}}_{i}} - \gamma \frac{\rho }{{\Delta {t^2}}}\frac{{{n^*} - {n^0}}}{{{n^0}}}. $

如图1所示为通用光滑壁面边界(generic smooth wall, GSW)模型，壁面形状由单层壁面粒子表示且壁面粒子中心精确位于几何表面. 定义距壁面0.5 ${l_0}$的平面为平衡面( ${l_0}$为流体粒子直径)，进入平衡面与壁面之间区域的流体粒子定义为近壁面流体粒子，该类粒子受到壁面斥力作用而重新回到平衡位置. 壁面斥力通过压力梯度力施加，表达式^[14]如下：

(9) $ {\left\langle {\nabla p} \right\rangle _i} = {\left\langle {\nabla p} \right\rangle _{i{\text{f}}}}+{\left\langle {\nabla p} \right\rangle _{i{\text{b}}}}, $

(10) $ {\left\langle {\nabla p} \right\rangle _{i{\text{b}}}} = \frac{{\rho \left| {{\bf{dr}}} \right|}}{{\Delta {t^2}}}\frac{{{{\boldsymbol{r}}_{i{\text{b}}}}}}{{\left| {{{\boldsymbol{r}}_{i{\text{b}}}}} \right|}}. $

式中： ${{\boldsymbol{r}}_{i{\text{b}}}}$、 ${\bf{dr}}$分别为从流体粒子i指向壁面的矢量和从流体粒子i指向平衡面的矢量，近壁面流体粒子 $ i $所受压力梯度力 $ {\left\langle {\nabla p} \right\rangle _i} $分为周围流体粒子作用部分 $ {\left\langle {\nabla p} \right\rangle _{i{\text{f}}}} $以及壁面边界作用部分 $ {\left\langle {\nabla p} \right\rangle _{i{\text{b}}}} $.

周围流体粒子作用部分通过式(5)计算，壁面边界作用部分通过引入距离函数(式(10))计算，具体推导过程可见文献[14]. GSW方法引入接触模型计算 ${\bf{dr}}$，从而避免直接使用式(10)造成的压力振荡问题^[12].

接触模型^[20]包含1个线性弹簧、1个阻尼以及1个摩擦滑块，如图2所示. 接触力分解为法向与切向，由以下公式^[14]计算：

(11) $ {{\boldsymbol{f}}_{i{\text{b}}}} = {\boldsymbol{f}}_{i{\text{b}}}^{\text{n}}+{\boldsymbol{f}}_{i{\text{b}}}^{\text{t}} , $

(12) $ {\boldsymbol{f}}_{i{\text{b}}}^{\text{n}} = - {k^{\text{n}}}{{{\boldsymbol{\delta}} }^{\text{n}}} - {\zeta ^{\text{n}}}{\boldsymbol{u}}_i^{\text{n}} ,$

(13) $ {{{\boldsymbol{\delta}} }^{\text{n}}} = 0.5l_i{\boldsymbol{e}}^{\rm{n}} - {{{\boldsymbol{r}}}_{{i}{\text{b}}}} - ({\boldsymbol{u}}_i^{\text{n}} - {\boldsymbol{u}}_{\text{b}}^{\text{n}})\Delta t, $

(14) $ {\boldsymbol{f}}_{i{\text{b}}}^{\text{t}} = \left\{ \begin{gathered} - {k^{\text{t}}}\left( {{\boldsymbol{u}}_i^{\text{t}} - {\boldsymbol{u}}_{\rm{b}}^{\text{t}}} \right)\Delta t - {\zeta ^{\text{t}}}{\boldsymbol{u}}_i^{\text{t}}, \;\;\; {\left| {{\boldsymbol{f}}_{i{\text{b}}}^{\text{t}}} \right| < \mu \left| {{\boldsymbol{f}}_{i{\text{b}}}^{\text{n}}} \right|} ; \\ - \left| {{\boldsymbol{f}}_{i{\text{b}}}^{\text{n}}} \right|{{ {\boldsymbol{u}}_i^{\text{t}}} / {\left| {{\boldsymbol{u}}_i^{\text{t}}} \right|, \;\;\; {\left| {{\boldsymbol{f}}_{i{\text{b}}}^{\text{t}}} \right| \geqslant \mu \left| {{\boldsymbol{f}}_{i{\text{b}}}^{\text{n}}} \right|} }}. \\ \end{gathered} \right. $

式中： $ {{{\boldsymbol{f}}}_{i{\text{b}}}} $表示边界对流体粒子产生的作用力；上标n与t分别代表法向与切向，下标i、b分别表示粒子i和壁面边界； $ {{\boldsymbol{\delta}} ^{\text{n}}} $为流体粒子i的法向移动矢量；l_i为i粒子的粒子直径；eⁿ为单位法向量. $ k $、 $ \zeta $分别为弹簧系数与阻尼系数，依据具体算例确定.

根据运动学定律，粒子向平衡面移动距离即

(15) $ \left| {{\bf{dr}}} \right| = |{{{\boldsymbol{f}}}_{i{\text{b}}}}|\Delta {t^2}/{m_i}. $

式中：m为粒子质量.

对于壁面附近液体粒子，其核近似域被壁面边界截断，粒子数密度偏小，导致各个算子计算不准. GSW基于经验公式构造壁面补偿函数，虽然在二维情况下补偿效果较好，但对于三维算例，较难得到合适的补偿函数. 本研究耦合无表面粒子判定技术（non-surface detection technique，NSD），引入虚拟粒子，补偿粒子数密度缺失的流体粒子，如图3所示.

虚拟粒子从数值层面补足粒子数密度，补偿作用与流体粒子受到的核截断程度相关，即与壁面边界的形状相关，与壁面粒子的分布无关，在壁面加密后无须再引入虚拟粒子来解决边界缺失问题. 新的粒子数密度定义为 $ {\langle n''\rangle }_{i} $，表达式如下：

(16) $ {\langle n''\rangle }_{i}=\mathrm{max}\;({\langle n\rangle }_{i},{\langle n'\rangle }_{i}), $

(17) $ {\left\langle {n'} \right\rangle _i} = {n^0}+\sum\limits_{i,j \in {\rm{fluid}}}^{N_{\rm{r}}} {[w(\left| {{{{\boldsymbol{r}}}_{{ij}}}} \right|) - w({l_0})]} . $

式中： $ {\left\langle n \right\rangle _i} $为按原始方法（式(7)）计算出的粒子数密度； ${\left\langle {n'} \right\rangle _i}$为补足后的粒子数密度；式(17)右侧第2项考虑了粒子挤压时的额外补偿值，N_r为i粒子搜索域内的流体粒子数. 对于壁面附近流体粒子，若按照原始方法计算出的粒子数密度小于补足后的粒子数密度，采用补足后粒子数密度. 压力泊松方程对应修改为

(18) $ \begin{split} &\frac{2d}{\lambda {n}^{0}}({\langle n''\rangle }_{i}^{*}{p}_{i}-{\displaystyle \sum _{j\ne i}^{M}w({{r}}_{{ij}}){p}_{j}})=\gamma \frac{\rho }{\Delta {t}^{2}}\frac{{\langle n''\rangle }_{i}^{*}-{n}^{0}}{{n}^{0}}- \\& \qquad (1 - \gamma )\frac{\rho }{{\Delta t}}\nabla \; \cdot \;{{\boldsymbol{u}}}_i^*. \end{split} $

式中：M为i粒子搜索范围内的其余液体粒子数， ${p}_{i} $为粒子i的压力， ${{\boldsymbol{u}}}_i^* $为流体粒子i的速度. 具体推导过程可见文献[15]. 黏性项采用同样的方式拆分计算.

在采用GSW-NSD模型时，壁面粒子不参与压力计算，只通过法向量传递壁面信息，粒子尺度（疏密）与斥力大小无关，可以抽象看作节点. 据此特性提出适用于任意三维曲面的壁面粒子建模方法.

如图4(a)所示为一子弹头形状的几何模型，在三维建模软件中生成三维模型，将该模型通过网格划分软件先划分相对较粗的表面网格(见图4(b))，提取网格节点坐标并对应生成壁面粒子(见图4(c))，再针对形状复杂、曲率较大或异面干涉较强的表面，进行局部加密或嵌套加密. 原则上，壁面粒子间距应小于液体粒子直径. 划分表面网格只为获取位于几何面的粒子坐标，适应性较强的非结构三角形网格即可满足需求.

壁面粒子的生成不依赖网格的拓扑结构，即使是扭曲网格或负体积网格，只要节点分布足够密集，仍可利用本模型进行计算，网格质量对计算结果无影响. 本研究使用ICEM划分面网格，由于对网格质量要求不高，其他可在几何曲面上生成节点的常用软件或程序也可满足要求. 得到的壁面模型由单层粒子组成，叶片内部空心，且尖角以及大曲率表面都可以得到较好表达(见图4(c)).

壁面粒子不参与压力计算，其间距对计算结果无影响，只体现壁面的离散程度. 3个壁面粒子即可精确表达1个平面，对于曲面离散，壁面粒子数量越多，离散精度越高. 本研究提出一种壁面粒子局部加密技术，即采用更多的壁面粒子来精确离散主曲率较大的曲面，如图5所示为离心泵蜗舌、倒圆角以及弯扭叶片的局部加密. 采用较少的粒子离散平面或主曲率较小的表面，从而在保证高离散精度的同时，减小计算量.

在本研究使用的光滑壁面(GSW)模型中，斥力的大小主要由流体粒子到壁面边界的距离决定，与壁面粒子的疏密无关；壁面粒子加密可以使流体粒子流经曲面时受到过渡更加光滑的斥力，从而减小非物理的压力振荡，提高计算稳定性并增强收敛性.

流体机械存在较多薄壁结构，如轮盘、轮盖和叶片等，在计算法向量时，薄壁结构端面粒子易受附近壁面（异面）粒子的干涉，导致法向量计算方向出错，如图6(a)所示.

常规基于网格法的加密技术由于拓扑结构的限制，较难实现相邻网格模块之间的网格（粒子）加密突跃. 本研究基于拉格朗日法无固定拓扑结构的特性，提出嵌套加密方法，可以将局部结构进行独立加密（见图6(b)）并与其他结构直接组合，实现该薄壁结构的高精度离散与法向量的准确计算. 嵌套加密的层数与具体几何形状以及法向量计算精度相关，本研究建议端面离散一般需6层以上.

为了评价离散后粒子模型的几何精度，定义法向光滑角度(normal smooth angle)为评价指标. 如图7所示，以二维曲线为例，对于壁面粒子 $ i $，沿着曲线切线方向在其左、右两侧分别搜索最近壁面粒子 $ {i_{\text{L}}} $、 $ {i_{\text{R}}} $，依次计算 $ i $与 $ {i_{\text{L}}} $、 $ {i_{\text{R}}} $的法向量夹角，平均化后作为 $ i $粒子法向光滑角度 $\theta _i^{{\text{NSA}}}$.

各个壁面粒子的法向光滑角度取标准差即得到该壁面模型的离散精度评估参数 ${\varepsilon ^{{\text{NSA}}}}$. 在三维情况下，沿曲面切平面搜索4个方向. 表达式如下：

(19) $ \left. \begin{array}{l} {\theta _i^{{\text{NSA}}} = \dfrac{1}{D}\displaystyle \sum\limits_{k = 1}^D {\arccos \; \dfrac{{{{{\boldsymbol{n}}}_{i}} \cdot {{{\boldsymbol{n}}}_{{k}}}}}{{\left| {{{{\boldsymbol{n}}}_{i}}} \right| \left| {{{{\boldsymbol{n}}}_{{k}}}} \right|}}} } ,\\ {{\varepsilon ^{{\text{NSA}}}} = \sqrt {\dfrac{1}{N}\displaystyle \sum\limits_{i = 1}^N {{{(\theta _i^{{\text{NSA}}} - {{\bar \theta }^{{\text{NSA}}}})}^2}} } } . \end{array} \right\} $

式中： $N$为壁面粒子总数； $D$为搜索方向数目，本研究三维情况取4； ${{{\boldsymbol{n}}}_{i}} $为i粒子的法向量； ${{\bar \theta }^{{\text{NSA}}}} $为所有壁面粒子的法向光滑角度的平均值. 离散复杂曲面的粒子越多，离散精度越高，离散精度评估参数 ${\varepsilon ^{{\text{NSA}}}}$越小. 在复杂曲面布置更密的壁面粒子不仅能更精确地表达壁面形状，也可以得到更加准确的法向量，使壁面附近流动更加稳定，减小壁面离散几何误差导致的流体粒子振荡. 法向光滑角度与计算误差的关系将会在后续研究中报道，本研究仅给出能满足计算需求的参考值.

对于如图5所示的离心泵模型，采用本研究方法加密处理后的壁面粒子数为75377，离散精度为 ${\varepsilon ^{{\text{NSA}}}} = 7.11 \times {10^{ - 2}}$. 在采用单一壁面粒子间距离散时，达到该精度需要的粒子数为117997，较本研究方法多36%.

为了验证上述模型处理复杂曲面的准确性与稳定性，本研究计算了2个带有典型复杂几何特征的静水算例，如蘑菇形（见图8(a)）和圆锥形（见图8(b)），分别考察本模型精确表达曲率变化较大的凹凸面和处理尖锐顶点的能力.

蘑菇和圆锥表面采用加密布置，壁面粒子间距为1 mm，水箱壁面采用稀疏布置，粒子间距为2 mm. 蘑菇算例与圆锥算例的总壁面粒子数依次为31331、29504，模型离散精度 ${\varepsilon ^{{\text{NSA}}}}$为 $5.71 \times {10^{ - 2}}$、 $5.55 \times {10^{ - 2}}$. 壁面粒子模型如图9所示.

X方向半截面的压力分布如图10所示. 图中，p为压力. 静水压力分布连续且与理论值一致，蘑菇算例的不规则凹凸面被液体粒子均匀贴合. 圆锥算例尖角并未导致液体粒子非物理振荡，2个算例的自由面在长时间计算过程中均保持水平稳定. 蘑菇与圆锥算例的底部压力误差分别为0.6%、1.3%.

为了验证本研究算法及离散模型的准确性，计算溃坝撞击挡板算例并与实验数据、其他高精度算法的结果进行对比. 模型几何参数、仿真参数以及压力测点位置与文献[21]一致. 水箱和挡块壁面的粒子间距分别为5、10 mm，壁面粒子数和液体粒子数分别为273323、43120. t=0.4 s时的实验图片和t=0.439 s时本研究结果的对比如图11所示，仿真结果仅展示了Y方向中截面视图的液体粒子.

测点处压力随时间波动如图12所示，相较于其他高精度算法如MPS-PW-IS和PMPS，本研究结果更加接近实验结果，且在t=3 s后压力振荡尖峰较少，结果相对更稳定.

为了验证复杂流动情况下，经本方法离散出的壁面模型仍具有较好的力学性能，选取离心泵内的诱导轮结构进行计算，几何模型如图13所示. 叶片厚度较小，与诱导轮直径比约为1∶47，且该厚度沿轴向逐渐变化(最小厚度0.969 mm，最大厚度1.520 mm). 诱导轮叶片呈螺旋形，叶根、叶顶及叶片端面处的法向量计算易受轮毂及叶片工作面壁面粒子的干涉.

利用本研究方法对诱导轮不同部件进行局部加密，轮毂粒子和水箱壁面粒子间距均为0.8 mm，螺旋叶片粒子间距为0.6 mm. 对于尺寸较小的叶片侧壁，采用间距为0.2 mm的壁面粒子进行嵌套加密，加密结果如图14所示. 诱导轮模型壁面粒子数为80584，圆柱形水箱壁面粒子数为27417，整体模型离散精度为 $8.10 \times {10^{ - 2}}$.

螺旋叶片沿Z轴正方向转动，持续向水箱底面输运流体并造成压力冲击，转速为200 r/min. 底部流体经由叶片与水箱间隙流回自由面，尺寸较小的间隙使得叶片后沿易产生多个叶尖漩涡，加剧流动的复杂程度. 流体机械内部流动多为湍流，本研究主要关注前处理工作并验证离散方法和壁面模型有效性，暂未考虑湍流模型，相关研究内容将在后续工作中报道. 流体运动黏度为5.0×10⁻⁴ m²/s，密度为1000 kg/m³，

为了对比计算结果，采用有限体积法对同参数条件下该模型进行计算，通过网格无关性验证的网格数为752万. 2种方法针对流体域的离散过程如图15所示.

传统网格方法往往须将流体域划分为旋转流域和静止流域，并建立动静交界面实现2个流域的信息传递，如图15（b）所示. 无网格粒子法利用带有物性的液体粒子整体离散流域，液体粒子依据NS方程流动到任意地方，再与当地周围的流体和固体粒子建立力学关系，无须将流体域分为动/静模块，也因此不涉及动静交界面.

初始流体域可以利用重力和粒子的自适应性对本方法建立的高精度壁面模型进行自动填充布置. 对于诱导轮算例，保持叶片静止，在叶轮上方布置足量的液体粒子并开始自由下落，液体粒子在重力作用下逐渐填充叶片间隙，当流体基本静止即平均粒子数密度稳定时，提取此时流体粒子坐标并作为初始分布.

在水箱底部、叶片表面设置压力监测点，坐标依次为(20.00,0,5.00)、(4.09,−10.78,9.02)、(4.06,−20.00,10.00) mm，2种方法计算出的压力随时间波动如图16、17所示，叶片表面监测点1处MPS和FVM法计算出的平均压力分别为399.46、409.30 Pa，误差为2.4%，叶片表面监测点2处MPS和FVM法计算出的平均压力分别为401.9、419.3 Pa，误差为4.1%. 流动进入稳定循环后的流场对比如图18所示.

从定量比较看，2种方法在监测点的压力波动频率、幅值一致；压力值误差较小，水箱底部平均误差为2.31%，峰值误差为4.90%，叶片表面压力误差依次为4.1%和2.3%. 从流动细节看，2种方法计算出的流场分布较接近，图示时刻5个叶尖漩涡都得到了较好的体现. 在本方法计算过程中液体粒子紧贴壁面，较好地表达了壁面形状，说明在复杂流动情况下，经本方法离散后的壁面模型仍可保持较好的力学性能.

本研究采用无网格MPS法，基于通用光滑壁面边界模型(GSW)与无表面粒子判定技术(NSD)，针对目前无网格粒子法较难处理流体机械中常见复杂型面的问题，建立了一套任意三维曲面的高精度通用离散方法，结论如下：

（1）提出基于MPS-GSW-NSD的分级离散方法. 先整体离散复杂壁面，借助面网格建立高精度壁面粒子模型，再向壁面模型内填充布置流体粒子；整体离散复杂流域，无须分块处理动/静流体域.

（2）提出离散过程的嵌套加密技术. 区别对待不同尺寸及曲率复杂型面的壁面粒子分布密度，避免了异面粒子干涉，有效提升了壁面法向量计算的准确性. 针对流体机械中的特征几何形状如弯扭叶片和薄壁结构，建立了一套有效的离散方法.

（3）提出单一型面内局部加密技术. 基于GSW模型的壁面粒子尺度对计算无影响的特点，在单一几何型面内粒子离散分布密度任意可调，实现了计算精度与计算效率的可调节.

（4）建立粒子离散精度评价方法. 基于法向光滑角度建立离散精度的评价方法，为评估和优化粒子离散精度提供了分析工具和优化依据.

此外，本研究通过计算包含不规则凹凸面和尖点的三维静水算例，验证了本方法处理复杂三维几何的能力. 通过计算溃坝撞击挡板算例，验证了离散后的壁面模型的准确性；通过计算诱导轮水中旋转算例，证明了在复杂流动情况下，本模型仍能保证正确的壁面作用以及压力计算的稳定性. 本研究仅给出某型号离心泵的泵壳、闭式叶轮和诱导轮等主要通流结构的离散评估方法，后续可将本方法扩展至密封等间隙流动结构，实现更加细致和全面的数值模拟.