中国前100家数字经济企业的运营效率测量：基于 DEA 和核密度估计的方法 张琳燕、 张玉萌、 杨坤、 张健

《Mathematics》：The Operational Efficiency Measurement of China’s Top 100 Digital Economy Firms: An Approach Based on DEA and Kernel Density Estimation Linyan Zhang, Yumeng Zhang, Kun Yang and Jian Zhang

　　**亮点与主要发现**

1. 开发了一种高斯-立方混合核函数，用于无网格配置方法求解四阶流函数方程。
2. 幽灵点方法提高了瞬态不可压缩流模拟中的数值稳定性和准确性。
3. 皮卡迭代有效地线性化了纳维-斯托克斯流函数公式中的非线性雅可比算子。
4. 三个基准测试证实了所提出方法的高精度和收敛性。

**主要发现的意义**

无网格方法消除了网格生成的需求，简化了复杂二维粘性流的计算建模。增强的稳定性和精度使得对非稳态流动现象进行长时间模拟成为可能。该方法为工程流体动力学应用提供了一种高效的传统基于网格方法的替代方案。

**摘要**

本文提出了一种无网格配置技术，用于纳维-斯托克斯方程的四阶瞬态流函数公式。该技术采用混合核函数，并通过幽灵点方法和皮卡迭代进行改进。这种流函数方法固有的未知数减少简化了解决过程。引入涡度和流函数使得耦合的、时间依赖的纳维-斯托克斯系统能够被重新表述为单变量的四阶偏微分方程。使用高斯-立方向后代换方法和时间差分方法来求解相应的方程，其中非线性部分通常通过皮卡迭代方法转化为线性方程。通过模拟三个流动案例，证明了该方案的可行性。

**1. 引言**

近几十年来，由于数学建模的推动，对不可压缩粘性流体流动的理解取得了显著进展。不可压缩流体流动的数学模型通常基于动量和质量守恒方程。这些方程描述了速度场、压力场及其演化过程，为理解和预测流体行为提供了基本框架。不可压缩流体流动的建模在航空[1]、机械工程[2]、化学过程[3]和生物工程[4]等各种工程应用中至关重要。通过数学建模，工程师可以优化设备设计、提高分离和反应效率，并预测各种与流动相关的现象。不可压缩粘性流体的多个特性（如流动稳定性、涡流、层流和湍流现象）可以通过数学分析得到研究。计算的进步使得包括有限元方法（FEM）[5]、有限差分方法（FDM）[6]和有限体积方法（FVM）[7]在内的数值方法成为模拟不可压缩流动的流行方法。这些数值算法能够有效解决复杂的流动问题，并促进对实际流动现象的理解。纳维-斯托克斯方程中固有的强非线性和速度-压力耦合带来了相当大的计算挑战。为了提高求解精度，四阶流函数公式作为一种简化替代方案而受到重视，用于模拟二维不可压缩粘性流动[8]。在速度-涡度公式中系统地用流函数表示速度变量，得到了一个四阶纯流函数控制方程。这种方法通过消除压力项，本质上也避开了由压力引起的计算难题。在二维不可压缩粘性流动中，所得的四阶方程、涡度-流函数公式和原始变量纳维-斯托克斯方程代表了可互换的数学框架[9,10]。因此，通过求解这些公式中的任何一个得到的数值结果都可以用来表征不可压缩流体流动现象。

本研究建立了一种高斯-立方向后代换方法（GC-BSM），结合了全隐式时间方案，用于求解二维不可压缩粘性流动中的自包含流函数公式。该方法具有三个新颖特性：
- 强形式的无网格实现，无需网格生成；
- 代数方程构造，绕过了数值积分过程；
- 与高阶流函数算子的天然兼容性。

本文的结构如下：第2节从数学上推导出纯流函数方程。第3节介绍了核心的高斯-立方向后代换算法。第4节通过三个测试案例验证了其准确性，第5节提供了结论性意见。

**2. 数学公式**

在流体动力学中，控制二维不可压缩粘性流动的数学模型由连续性方程和动量方程决定。边界假设为Lipschitz连续，以推导流函数公式：
(1)
在二维坐标中，
(2)
众所周知的纳维-斯托克斯方程明确表述如下：
(3)
(4)
其中u和v分别是x和y方向上的无量纲速度。计算域由边界界定。P表示压力，Re表示雷诺数。为了在动量方程中分离速度u和压力p，引入流函数ψ和涡度函数ω来简化系统。这些变量来自以下定义：
(5)
将这些关系代入方程(2)-(4)后，原始变量系统被严格转化为流函数-涡度公式：
(6)
(7)
由于在推导过程中消除了压力梯度项，计算的稳定性得到了提高。结合这两个方程，可以得到包含时间项的四阶纯流函数公式：
(8)
其中J表示雅可比算子，由以下公式定义：
(9)
所推导出的偏微分方程的四阶性质要求施加Neumann边界条件和Dirichlet边界条件，以确保数学系统的适定性。这些条件在下面给出了数学定义：
(10)
这里，dU/dx和dP/dy分别规定了Dirichlet和Neumann边界约束。边界上的单位外法向量由n表示。

本文研究了由方程(8)中的时间依赖公式和方程(10)中的边界约束所决定的二维瞬态不可压缩流动动力学。

**3. 解决算法**

3.1. 混合核函数

径向基函数（RBF）插值相比克里金方法具有网格独立的优势，因为它消除了变差图的需求[11]，特别适用于稀疏数据集。其更快的收敛性有助于在无网格偏微分方程（PDE）求解器中构建形状函数。然而，传统的RBF在大规模应用中会产生病态矩阵，从而指数级放大输入误差。Mishra的混合核通过高斯基函数和无形状参数的立方基函数的协同耦合[12]解决了这个问题，实现了通过谱互补性达到条件稳定性。Mishra的数值分析表明，这种混合核有效降低了插值矩阵的条件数，同时提高了算法稳定性。混合核是高斯核和立方径向基函数的线性组合：
(11)
这里，d表示点u和v之间的欧几里得距离，形状参数α控制高斯核的衰减特性。对于混合核，权重系数β调节了高斯基函数和立方基函数的相对贡献，从而调整它们的耦合强度。

3.2. 幽灵点方法

Chen等人[12,13,14]提出的幽灵点方法显著增强了径向基函数配置方法（RBFCMs），解决了传统Kansa方法[15]的关键局限性。传统方法将虚拟节点分布在边界和计算域内。相比之下，幽灵点方法将RBF中心策略性地均匀放置在扩展的幽灵域内，同时将边界点和配置点限制在边界上。图1展示了Kansa方法的标准节点分布模式，而图2展示了幽灵点配置，其中Δ（一个矩形域）完全包含了计算域Δ，使得Δ?Δ。数值验证[12]证实，这种方法不仅确保了解的稳定性，而且与经典实现相比，相对误差降低了40-60%。在此基础上，本研究采用了带有随机优化节点选择的幽灵点框架来进一步提高稳定性。

3.3. G-C BSM

核心的BSM方法将边界近似解纳入数值解中，以确定校正函数的展开系数。这些校正基系数由规定的边界条件决定。虽然近似函数本身满足边界约束，但校正函数确保数值解遵循均匀边界条件。首先，我们将区域分为三个部分：(1) 计算域Δ，(2) 幽灵点方法区域Γ，以及(3) 边界区域Β。分配给这三个区域的点数分别为NΔ、NΓ和NB，其中NΔ、NΓ和NB分别为Δ、Γ和Δ的元素数量。

为了说明，四阶纯流函数公式（方程(8)）和边界条件（方程(10)表示为简写形式：
(13)
其中Δ×?表示在边界上定义的微分算子，δ表示在内部域中作用的微分算子。Δ×δ表示方程(8)的左侧。我们考虑Dirichlet边界条件：
(14)
和Neumann边界条件：
(15)
流函数的数值近似表示为：
(16)
其中δ表示从皮卡迭代中得到的校正函数。这个函数通过数值解确保了控制方程的满足。边界近似函数和校正函数都是之前定义的基函数的线性组合：
(17)
(18)
这里，展开系数C可以通过将ψ代入边界条件方程(14)和(15)得到：
(19)
同时，为了保证数值解与规定的边界约束一致，与校正函数相关的混合核函数必须满足均匀边界条件：
(20)
方程(11)中的参数β通过以下分析过程确定：
(21)
(22)
确定参数β和展开系数后，通过方程(11)、(17)和(18)计算边界近似和混合核函数。然后将方程(16)中定义的流函数代入控制方程(13)，得到：
(23)
由于雅可比算子的非线性，采用皮卡迭代方案来线性化系统。

3.4. 皮卡迭代

在获得初始系数后，迭代计算确定随时间步变化的所有后续系数。使用时间步长Δt，在BSM框架内构建连续时间级的数值解。每个解是由边界近似和校正函数的叠加组成。如前所述，时间t处的流函数表示为：
(24)
这里，ψt表示时间t处的近似流函数，δt表示同一时间步的相应校正函数。这种关系与之前建立的框架一致：
(25)
(26)
其中，ψt的系数可以使用边界条件得到。而δt的系数将通过本节中的方法求解。为了书写方便，将其缩写为：
(27)
混合核函数的参数β应根据上述方法重新计算，δt可以使用皮卡迭代求解。首先，控制方程在时间上隐式求导，结果如下：
(28)
由于控制方程的非线性，采用皮卡迭代方案系统地线性化非线性算子：
(29)
这里，υt表示时间t处的离散数值解。将构成关系(24)和(27)代入控制方程(29)，非线性系统被重新表述为线性化方程：
(30)
这里，δt是通过在前一节描述的方法在时间t级施加边界条件得到的。在时间t级，迭代过程如下实现：
(1) 初始迭代：对于初始迭代，采用前一个时间步的数值解作为初始猜测。
(2) 后续迭代：在第二次迭代中，将上一次迭代的更新解作为新的初始值，重新计算系数。这个过程重复进行，直到第(n + 1)次迭代满足收敛标准：
(31)
在达到迭代收敛后，将解作为时间t级的最终数值解。

**4. 数值示例**

本节通过三个基准案例验证了所提出的方法。其中两个案例的解析解是已知的。对于缺乏显式解的问题，数值结果与已建立的文献进行了对比。所有模拟都在MATLAB R2024a环境下，使用Inter(R) Core (TM) i5-13600KF @5.1 GHz处理器和32 GB DDR5 6800 MHz内存上执行。在所有情况下，高斯核形状参数固定为某个值， ghost点域的长度设置为某个值，其中表示正方形域的边长。通过均方根（RMS）误差来量化算法在不同节点密度下的鲁棒性。为了便于计算比较，进一步定义了具有解析解的问题的最大误差。为了量化收敛速率，计算了在不规则节点上最大绝对误差范数的收敛阶数。示例1：研究了单位正方形域[0,1] × [0,1]内的二维行波流动。参考文献[16]提供了对比数值解（Re = 1）。在0.001的固定时间步长下，10次迭代后的数值结果（T = 0.01）显示速度分量的计算解与精确解之间有极好的一致性，如图3的云图所示。图3：T = 0.01时的对比云图。在（对应于1000次迭代）时，图4展示了速度分量的数值解与精确解之间的逐点绝对误差散点图。每个数据点代表绝对误差的空间分布。示例2：考虑了[0, 1] × [0, 1]正方形内的泰勒-格林涡旋。参考文献[16,17]给出了数值解。使用0.001的时间步长，模拟进行了500次迭代，总模拟时间为0.5。图6比较了速度分量的数值解与精确解。图7和图8展示了在不同形状参数值下，数值解与精确解之间速度分量的最大误差。形状参数显著影响插值矩阵的条件数和基于RBF的近似的准确性。不适当的形状参数值可能导致数值不稳定或过度平滑。在这项工作中，研究了数值解对形状参数的敏感性，并确定了一个适当的值以确保在模拟不可压缩粘性流动时的稳定性和高保真度。数值解在关键区域（如行波流动的波峰和波谷位置）与精确解高度一致，验证了捕获波动动力学的高保真能力。水平速度分量的波前形态和振幅分布与理论解非常吻合，证明了数值方法的可靠性。此外，如表1所示，计算过程中没有出现发散或非物理振荡，突出了算法的鲁棒稳定性。表1：示例1中不同点密度下的最大误差、RMS误差和计算时间的比较。表2研究了不同ghost点域对数值结果的影响。研究发现，虽然扩大ghost点域仅对精度有适度改善，但增加节点数会产生更显著的效果，这突出了节点密度比域扩展更为重要。示例3：考虑了单位正方形[0, 1] × [0, 1]内的方形腔体流动。边界条件在参考文献[18,19]中给出（Re = 100）。计算域定义为二维正方形区域。顶部边界施加了均匀的水平速度边界条件，垂直速度分量保持为零。所有其他边界（左、右和底部）均规定为无滑移条件，其中两个速度分量均满足相应条件。这种配置产生了一组完整的边界条件（也在图9中展示）。时间步长设置为0.01，计算时间至25，此时流动达到完全发展状态。Re = 100时的速度剖面（图10）与参考文献[18,19]的结果非常吻合。图11还绘制了沿和的速度剖面，将我们的结果与GFDM[20]、时间分裂方法[18]以及Ghia的[21]公认基准解进行了对比。表5：计算效率比较（T = 5）。结论：本研究提出了一种新颖的无网格径向基函数（RBF）基于点的配置框架，用于求解二维瞬态不可压缩粘性流动的四阶流函数公式。通过三个系统设计的基准案例（行波流动、泰勒-格林涡旋和盖子驱动的方形腔体流动），所提出的方法展示了以下优越特性：（1）所提出的无网格配置方法便于对具有不规则节点分布的瞬态二维不可压缩粘性流动进行计算建模。（2）ghost点方法增强了数值解的稳定性和准确性。当与隐式时间离散化和Picard线性化结合使用时，该方法在长时间模拟中保持了高精度。（3）所提出的方法在工程预测性方面表现出色，能够准确重建包括涡旋结构的时间发展和标志性速度分布在内的典型流动特征，与-established基准数据的偏差小于5%。后续研究将重点关注对所提出方法的局部改进，以提高计算效率。同时，引入了一个部分迎风方案，使改进后的方法能够应用于更复杂的流场。