本节描述了欧拉壁膜（EWF）模型使用的时间和空间差分方案，以及求解算法。

17.4.1 时间差分方案

17.4.2 空间差分方案

17.4.3 求解算法

17.4.4 耦合求解方法

17.4.1 时间差分方案

以下时间差分方案可用于欧拉壁膜（EWF）模型。

17.4.1.1 一阶显式方法

17.4.1.2 一阶隐式方法

17.4.1.3 二阶隐式方法

17.4.1.1 一阶显式方法

定义 $F$ 和 $G$ 如下，下标 $i$ 表示前一时间步的值：

$$\begin{array}{rl}& F_i\equiv {\left[-{\nabla }_s\cdot \left({\rho }_{l}h{\vec{V}}_l\right)+{\dot{m}}_s\right]}_i\\ & G_i\equiv {\left[-h{\nabla }_s{P}_L+{\rho }_{l}h{\vec{g}}_{\tau }+\frac{3}{2}{\vec{\tau }}_{fs}-\frac{3{\mu }_l}{h}{\vec{V}}_l+{\dot{\vec{q}}}_s+{\vec{\tau }}_{{\theta }_w}-{\nabla }_s\cdot \left[{\rho }_{l}h{\vec{V}}_l{\vec{V}}_l+{\bar{D}}_V\right]\right]}_i\end{array}\text{\hspace{1em}(17.39)}$$

我们拥有以下离散化的薄膜质量和动量方程，其中下标 $i+1$ 表示当前时间步的值，而 $\Delta t$ 则代表用于薄膜计算的时间步长。

$$\begin{array}{rl}\frac{{\left(\rho h\right)}_{i+1}-{\left(\rho h\right)}_i}{\Delta t}& =F_i\\ \frac{{\left({\rho }_{l}h{\vec{V}}_l\right)}_{i+1}-{\left({\rho }_{l}h{\vec{V}}_l\right)}_i}{\Delta t}& =-{\left(\frac{3{\mu }_l}{h}{\vec{V}}_l\right)}_{i+1}+G_i\end{array}\text{\hspace{1em}(17.40)}$$

随后，影片高度与速度的计算如下：

$$\begin{array}{rl}{h}_{i+1}& =\frac{{\left({\rho }_{l}h\right)}_i+\Delta t{F}_i}{{\left({\rho }_l\right)}_{i+1}}\\ {\left({\vec{V}}_l\right)}_{i+1}& =\frac{{\left({\rho }_{l}h{\vec{V}}_l\right)}_i+\Delta t{G}_i}{{\left({\rho }_{l}h\right)}_{i+1}+\Delta t\frac{3{\mu }_l}{h_{i+1}}}\end{array}\text{\hspace{1em}(17.41)}$$

上述方程组完成了显式差分方案，其中首先根据前一薄膜时间步长计算的$F_i$值来计算薄膜高度；然后利用最新的薄膜高度和前一薄膜时间步长计算的$G_i$值来计算薄膜速度。

17.4.1.2 一阶隐式方法

在显式方法中，$F$和$G$的评估是基于前一时间步长的薄膜高度和速度矢量进行的。为了提高显式方法的准确性，引入了一阶隐式方法，在该方法中，$F$和$G$值在薄膜时间步长内的迭代循环中进行更新。这种方法可以描述为一个预测-校正过程。首先，即在预测步骤中，使用显式方案来计算薄膜高度和速度矢量，

预测器：

$$\begin{array}{rl}{h}_{i+1}^0& =\frac{{\left({\rho }_{l}h\right)}_i+\Delta t{F}_i}{{\rho }_{li+1}^0}\\ {\left({\vec{V}}_l\right)}_{i+1}^0& =\frac{{\left({\rho }_{l}h{\vec{V}}_l\right)}_i+\Delta t{G}_i}{{\left({\rho }_{l}h\right)}_{i+1}^0+\Delta t\frac{3}{h_{i+1}^0}}\end{array}\text{\hspace{1em}(17.42)}$$

上标0表示迭代循环的第一步。

校正器：

使用最新的液膜高度和速度向量来更新$F$和$G$；随后重新计算液膜高度和速度。

$$\begin{array}{rl}{h}_{i+1}^{n+1}& =\frac{{\left({\rho }_{l}h\right)}_i+\Delta tF\left[h_{i+1}^n,{\left({\vec{V}}_l\right)}_{i+1}^n\right]}{{\rho }_{li+1}^{n+1}}\\ {\left({\vec{V}}_l\right)}_{i+1}^{n+1}& =\frac{{\left({\rho }_{l}h{\vec{V}}_l\right)}_i+\Delta tG\left[h_{i+1}^n,{\left({\vec{V}}_l\right)}_{i+1}^n\right]}{{\left({\rho }_{l}h\right)}_{i+1}^{n+1}+\Delta t\left(3{\mu }_l/h_{i+1}^{n+1}\right)}\end{array}\text{\hspace{1em}(17.43)}$$

上标 $n+1$ 和 $n$ 分别表示当前和前一次迭代。迭代过程在满足以下收敛准则时结束：

$$\begin{array}{rl}& \left|h_{i+1}^{n+1}-h_{i+1}^n\right|\le \epsilon \\ & \left|{\left(u_{\alpha }\right)}_{i+1}^{n+1}-{\left(u_{\alpha }\right)}_{i+1}^n\right|\le \epsilon \end{array}\text{\hspace{1em}(17.44)}$$

其中，$u_{\alpha }$ 表示速度矢量的每个分量。

17.4.1.3 二阶隐式方法

在上述讨论的显式和一阶隐式方法中，时间差分仅具有一阶精度。下面介绍一种二阶隐式方法。其迭代过程与一阶隐式方法非常相似，但在预测步骤中使用了两个时间步长（即 $i$ 和 $i-1$）进行时间差分。

预测步骤：

$$\begin{array}{rlrl}& h_{i+1}^0=& & \frac{{\left({\rho }_{l}h\right)}_i+\left(\Delta t/2\right)\left(3F_i-F_{i-1}\right)}{{\rho }_{li+1}^0}\\ & {\left({\vec{V}}_l\right)}_{i+1}^0=& & \frac{{\left({\rho }_{l}h{\vec{V}}_l\right)}_i+\left(\Delta t/2\right)\left(3G_i-G_{i-1}\right)}{{\left({\rho }_{l}h\right)}_{i+1}^0+\Delta t\left(3{\mu }_l/h_{i+1}^0\right)}\end{array}\text{\hspace{1em}(17.45)}$$

校对者：

$$\begin{array}{rl}& h_{i+1}^{n+1}=\frac{{\left({\rho }_{i}h\right)}_i+\left(\Delta t/2\right)\left\{F\left[h_{i+1}^n{\left({\vec{V}}_l\right)}_{i+1}^n\right]+F\left[h_{i+1}^{n-1}{\left({\vec{V}}_l\right)}_{i+1}^{n-1}\right]\right\}}{{\rho }_{li+1}^{n-1}}\\ & {\left({\vec{V}}_l\right)}_{i+1}^{n+1}=\frac{{\left({\rho }_{l}h{\vec{V}}_l\right)}_i+\left(\Delta t/2\right)\left\{G\left[h_{i+1}^n{\left({\vec{V}}_l\right)}_{i+1}^n\right]+G\left[h_{i+1}^{n-1}{\left({\vec{V}}_l\right)}_{i+1}^{n-1}\right]\right\}}{{\left({\rho }_{l}h\right)}_{i+1}^{n+1}+\Delta t\left(3{\mu }_l/h_{i+1}^{n+1}\right)}\end{array}\text{\hspace{1em}(17.46)}$$

迭代过程终止于以下收敛标准：

$$\begin{array}{rl}& \left|h_{i+1}^{n+1}-h_{i+1}^n\right|\le \epsilon \\ & \left|{\left(u_{\alpha }\right)}_{i+1}^{n+1}-{\left(u_{\alpha }\right)}_{i+1}^n\right|\le \epsilon \end{array}\text{\hspace{1em}(17.47)}$$

其中 $u_{\alpha }$ 表示速度矢量的各个分量。

薄膜能量和被动标量方程的离散表达式可以采用与上述类似的方式推导。

17.4.2 空间差分格式

计算方程 17.39（第 848 页）中定义的 $F$ 和 $G$ 需要对各个量在面中心处的空间梯度进行评估。应用格林-高斯定理，标量 $\phi$ 在面中心 0 处的梯度计算如下：

$$\nabla \phi =\frac{1}{A}\sum _e{\phi }_e{\overrightarrow{l}}_e\text{\hspace{1em}(17.48)}$$

在文本中，下标 $e$ 表示薄膜面每条边的中点，$A$ 表示面的面积，$\overrightarrow{l}$ 表示长度向量（其大小为边的长度，方向垂直于边），如图17.5所示：

图17.5：空间梯度

显然，梯度计算的核心在于如何获取 ${\phi }_e$。以下描述了当前实现中所采取的步骤。

(1) 初级梯度计算

初级梯度是通过格林-高斯定理简单估计边中心值来计算的，具体如下：

$${\phi }_e=\frac{1}{2}\left({\phi }_P+{\phi }_E\right)\text{\hspace{1em}(17.49)}$$

$$\nabla {\phi }_P=\frac{1}{A}\sum _e{\phi }_e{\overrightarrow{l}}_e\text{\hspace{1em}(17.50)}$$

(2) 边缘中心值的重构

随后，边缘中心值通过一阶迎风或二阶迎风方案进行计算，

一阶迎风：${\phi }_e^{\ast }={\phi }_{UP}$

二阶迎风：${\phi }_e^{\ast }={\phi }_{UP}+\overrightarrow{d}{r}_{UP}\cdot \nabla {\phi }_P$

其中，${\phi }_{UP}$ 是 $\phi$ 在迎风面中心的值；$\overrightarrow{d}{r}_{UP}$ 是从迎风面中心到边缘中心的距离向量。

(3) 最终梯度计算

最后，面中心的梯度计算如下，

$$\nabla \phi =\frac{1}{A}\sum _e{\phi }_e^{\ast }{\overrightarrow{I}}_e\text{\hspace{1em}(17.51)}$$

17.4.3. 解决方案算法

尽管薄膜运动始终随时间追踪，EWF模型既可以用于稳态流模拟，也可以用于瞬态流模拟。

17.4.3.1 稳态流

假设尽管流动影响其运动，薄膜本身并不影响流场。在这种情况下，流场已经收敛并保持不变（冻结流场）。薄膜发展所经过的时间通过薄膜时间步长 $\Delta t$ 增加。这个时间步长可以是

用户指定的常数，也可以使用自适应时间步进方法自动计算，该方法根据最大Courant数，由用户指定的增加和减少因子来决定，具体如下：

• 当计算出的最大Courant数超过用户指定的值时，薄膜时间步长会通过减少因子 $F_{dec}$ （默认：2）减少。

• 当计算出的最大Courant数小于用户指定值的一半时，薄膜时间步长会通过增加因子 $F_{inc}$ （默认：1.5）增加。

17.4.3.2 瞬态流

在每个流动时间步长 $\Delta t_{flow}$ 内，使用若干个薄膜子时间步长来推进薄膜时间，使其与流动的物理时间同步。薄膜子时间步长的确定方法如下，

$$\Delta t=\frac{\Delta t_{\text{flow }}}{N_{\text{film }}}\text{\hspace{1em}(17.52)}$$

$N_{\text{film}}$ 表示薄膜时间步数；$\Delta t_{\text{flow}}$ 表示流动时间步长。

通过调整薄膜时间步数 $N_{\text{film}}$ 并使用用户指定的增加和减少因子 $F_{\text{inc}}$ 和 $F_{\text{dec}}$，实现了自适应时间步进：

• 当计算出的最大库朗数超过用户指定的值时，薄膜时间步数 $N_{\text{film}}$ 将增加 $F_{\text{dec}}$ 倍，从而使薄膜时间步长减少 $F_{\text{dec}}$ 倍。

• 当计算出的最大库朗数小于用户指定值的一半时，薄膜时间步数 $N_{\text{film}}$ 将减少 $F_{\text{inc}}$ 倍，从而使薄膜时间步长增加 $F_{\text{inc}}$ 倍。

需要注意的是，薄膜时间步数可减少到的最小值为 2。

由于薄膜计算在每个流动时间步长的末尾开始，因此薄膜对流动的影响会滞后一个流动时间步长。

17.4.4 耦合求解方法

默认情况下，欧拉壁面薄膜模型中薄膜质量和动量方程是顺序求解的。在薄膜高度和速度紧密耦合的情况下，分离求解方法无法解析薄膜的特定特征，例如自由下落薄膜的波浪表面。

在耦合求解方法中，薄膜质量和动量方程并行求解，其中前一时间步的薄膜高度和速度解用于下一时间步中两个方程的求解。

为了提高求解稳定性，耦合求解方法中对用于薄膜曲率和表面张力力计算的变量应用了平滑处理。变量平滑通过在迭代操作中混合相邻面的变量值来实现：

$${\phi }_s=\frac{1}{W}\sum _{i=1}^n\left({\phi }_0\left(1-\alpha \right)+{\phi }_i\alpha \right)\text{\hspace{1em}(17.53)}$$

在给定面的薄膜变量$\phi$中，${\phi }_0$和${\phi }_s$分别表示原始值和平滑值；${\phi }_i$表示在第$i^{\text{th}}$相邻面（直接或间接）上的$\phi$值，如图17.6所示：面值平滑（第853页）；$\alpha$是一个平滑因子，范围从0到1（默认值为$0.5$）；$w$是由平滑算法确定的权重因子。

图17.6：面值平滑

加权平均操作会根据平滑级别输入参数（范围从1到4，默认值为2）指定的次数重复进行。