一种界面锐化技术 Thinc 方法 -- CFD 文献阅读整理（8）

THINC 方法（Tangent of Hyperbola for Interface Capturing ，双曲正切界面捕捉法）是一种代数型 VOF（Volume of Fluid）方法 ，属于代数重构类界面捕捉算法 。其核心思想是利用双曲正切函数 来重构体积分数在单元内的连续分布，从而实现界面的平滑表示与可控耗散。

有限体积，间断有限元，单相流，多相流，都可以用到 THINC 方法。它的主要作用是让间断算的更加 sharp.

这里的 $\alpha (x)$就是 THINC 函数，可以使用如下的定义[1]

${\alpha }_i(x)=\frac{{\alpha }_{\max }}{2}\left\{1+\gamma \tanh \left[\beta \left(\frac{x-x_{i-1/2}}{x_{i+1/2}-x_{i-1/2}}-{\tilde{x}}_i\right)\right]\right\}+{\alpha }_{\min }$式中： ${\tilde{x}}_i$为双曲正切函数过渡层的中点，这个的值每次是不一样的，需要使用 ${\int }_{I_i}{\alpha }_i(x)/(\Delta x)dx={\bar{\alpha }}_i$算出来； $\beta$控制界面过渡层厚度，一般取 3.5，能保证界面控制在两个网格厚度； $\gamma$与界面方向相关：
$\gamma =\{\begin{array}{llc}1,& \text{if }{\alpha }_{i+1/2}<{\alpha }_{i-1/2}-1,& \text{otherwise}\end{array}$
取 ${\alpha }_{R,i-1/2}=\alpha (x_{i-1/2})$和 ${\alpha }_{L,i+1/2}=\alpha (x_{i+1/2})$可求出左右界面的体积分数：
$\begin{array}{rl}{\alpha }_{L,i+1/2}& =\frac{1}{2}{\alpha }_{\max }\left(1+\gamma \frac{\tanh \beta +C}{1+C\tanh \beta }\right)+{\alpha }_{\min }\\ {\alpha }_{R,i-1/2}& =\frac{1}{2}{\alpha }_{\max }(1+\gamma C)+{\alpha }_{\min }\end{array}$

其中：

$\begin{array}{rl}{\alpha }_{\min }& =\min ({\alpha }_{i-1},{\alpha }_{i+1})\\ {\alpha }_{\max }& =\max ({\alpha }_{i-1},{\alpha }_{i+1})-{\alpha }_{\min }\\ B& =\exp \left[\gamma \beta \left(2\frac{{\alpha }_i-{\alpha }_{\min }+\epsilon }{{\alpha }_{\max }+\epsilon }-1\right)\right]\\ C& =\left(\frac{B}{\cosh \beta }-1\right)/\tanh \beta \end{array}$

其中：$B,C$ 是临时中间变量， $\epsilon$是一个小量，取 $1.0\times 10^{-12}$。为保证体积分数单调性，需要满足：
$\text{if }\text{ sign}[({\alpha }_{i+1}-{\alpha }_i)({\alpha }_i-{\alpha }_{i-1})]=-1$
则：
${\alpha }_{L,i+1/2}={\alpha }_{R,i-1/2}={\alpha }_i$

实际应用的时候，还应当结合BVD(boundary variation diminishing) 方法来用，这样可以确保程序的稳定[2]

---

其实以上内容基本都是简单的摘抄，接下来的内容比较关键。 ${\tilde{x}}_i$的计算需要使用限制条件 ${\int }_{I_i}{\alpha }_i(x)/(\Delta x)dx={\bar{\alpha }}_i$，但是这个 ${\alpha }_i(x)$里面带一个 $\tanh$所以没那么好积分。我算了好半天才搞清楚。

最终的计算结果为

${\tilde{x}}_i=\frac{1}{\beta }\mathrm{arctanh}\left(\frac{1-\exp (\frac{2\beta \theta }{{\alpha }_{\max }}({\alpha }_i-{\alpha }_{\min })-\log (\cosh \beta )-\beta \theta )}{\tanh \beta }\right)$

详细的推导过程在这里：

其对应的latex代码如下：

\begin{equation}
    \sinh x = \frac{e^x-e^{-x}}{2},\quad \cosh x = \frac{e^x+e^{-x}}{2}
\end{equation}
\begin{equation}
     \cosh(x-y) = \cosh x \cosh y - \sinh x \sinh y
\end{equation}


\begin{equation}
     \begin{aligned}
               \alpha_i &= \frac{2\,\beta \,\alpha_{\min}+\beta\alpha_{\max}-\beta\alpha_{\max}\theta+\alpha_{\max}\theta\,\log \left(\tanh \left(\beta \,\tilde{x}_i\right)+1\right)-\alpha_{\max}\,\theta\log \left(\tanh \left(\beta \,{\left(\tilde{x}_i-1\right)}\right)+1\right)}{2\,\beta }\\ 
               &= \alpha_{\min} + \frac{\alpha_{\max}}{2}(1-\theta) + \frac{\alpha_{\max}\theta}{2\beta} \log \frac{\tanh(\beta \tilde{x}_i)+1}{\tanh(\beta (\tilde{x}_i-1))+1}\\ 
               &= \alpha_{\min} + \frac{\alpha_{\max}}{2}(1-\theta) + \frac{\alpha_{\max}\theta}{2\beta} \log \left( \frac{\frac{2e^{\beta \tilde{x}_i}}{e^{\beta \tilde{x}_i}+e^{-\beta \tilde{x}_i}}}{\frac{2e^{\beta (\tilde{x}_i-1)}}{e^{\beta (\tilde{x}_i-1)}+e^{-\beta (\tilde{x}_i-1)}}}\right)\\ 
               &= \alpha_{\min} + \frac{\alpha_{\max}}{2}(1-\theta) + \frac{\alpha_{\max}\theta}{2\beta} \log \left( e^{\beta}\cdot\frac{e^{\beta (\tilde{x}_i-1)} + e^{-\beta (\tilde{x}_i-1)}}{e^{\beta \tilde{x}_i} + e^{-\beta \tilde{x}_i}}\right)\\ 
               &= \alpha_{\min} + \frac{\alpha_{\max}}{2}(1-\theta) + \frac{\alpha_{\max}\theta}{2} + \frac{\alpha_{\max}\theta}{2\beta} \log \frac{\cosh(\beta(\tilde{x}_i-1))}{\cosh(\beta \tilde{x}_i)}\\
               &= \alpha_{\min} + \frac{\alpha_{\max}}{2}(1-\theta) + \frac{\alpha_{\max}\theta}{2} + \frac{\alpha_{\max}\theta}{2\beta} \log \frac{\cosh(\beta\tilde{x}_i)\cosh(\beta)-\sinh(\beta\tilde{x}_i)\sinh(\beta)}{\cosh(\beta \tilde{x}_i)}\\
               &= \alpha_{\min} + \frac{\alpha_{\max}}{2}(1-\theta) + \frac{\alpha_{\max}\theta}{2} + \frac{\alpha_{\max}\theta}{2\beta} \log \left( \cosh(\beta) - \sinh(\beta)\tanh(\beta\tilde{x}_i) \right)\\
               &= \alpha_{\min} + \frac{\alpha_{\max}}{2}(1-\theta) + \frac{\alpha_{\max}\theta}{2} + \frac{\alpha_{\max}\theta}{2\beta} \log(\cosh\beta) + \frac{\alpha_{\max}\theta}{2\beta} \log\big(1-\tanh\beta\,\tanh(\beta\tilde{x}_i)\big)\\
     \end{aligned}
\end{equation}


\begin{equation}
         A = \alpha_i - \alpha_{\min} - \frac{\alpha_{\max}}{2} - \frac{\alpha_{\max}\theta}{2\beta}\log(\cosh\beta)
\end{equation}
\begin{equation}
     \begin{aligned}
     \log\big(1-\tanh\beta\tanh(\beta\tilde{x}_i)\big) = \frac{2\beta}{\alpha_{\max}\theta} A
        &\quad\Longrightarrow\quad
        1 - \tanh\beta \tanh(\beta \tilde{x}_i) = \exp\Big(\frac{2\beta}{\alpha_{\max}\theta} A\Big)\\
        &\quad\Longrightarrow\quad \tanh(\beta \tilde{x}_i) = \frac{1 - \exp\Big(\frac{2\beta}{\alpha_{\max}\theta} A\Big)}{\tanh\beta}
     \end{aligned}
\end{equation}
\begin{equation}
     \begin{aligned}
          \tilde{x}_i &= \frac{1}{\beta} \operatorname{arctanh} \left( \frac{1 - \exp\Big(\frac{2\beta}{\alpha_{\max}\theta} A\Big)}{\tanh\beta} \right) \\
          &= \frac{1}{\beta} \operatorname{arctanh} \Bigg( \frac{1 - \exp\Big(\frac{2\beta}{\alpha_{\max}\theta}\Big(\alpha_i - \alpha_{\min} - \frac{\alpha_{\max}}{2} - \frac{\alpha_{\max}\theta}{2\beta} \log(\cosh\beta)\Big)\Big)}{\tanh\beta} \Bigg) \\
          &= \frac{1}{\beta} \operatorname{arctanh} \Bigg( \frac{1 - \exp\Big(\frac{2\beta}{\alpha_{\max}\theta}(\alpha_i - \alpha_{\min}) - \log(\cosh\beta) - \frac{\beta}{\theta}\Big)}{\tanh\beta} \Bigg)\\ 
          &= \frac{1}{\beta} \operatorname{arctanh} \Bigg( \frac{1 - \exp\Big(\frac{2\beta\theta}{\alpha_{\max}}(\alpha_i - \alpha_{\min}) - \log(\cosh\beta) - \beta\theta\Big)}{\tanh\beta} \Bigg)
     \end{aligned}
\end{equation}

为了验证我的计算结果是否正确，又编写了如下的matlab程序。（因为过于复杂，matlab不会处理这个符号运算，所以使用数值的方法计算。不放心的话可以多试几组数）

clear
syms x
uavem = 10;
uave = 8;
uavep = 6;
xa = -1 ;
xb = 1;
beta = 1.5;
Umin = min(uavem,uavep);
Umax = max(uavem,uavep)-Umin;
theta = sign(uavep-uavem);

% 计算的 \tilde{x}_i
x_tilde_i = 1/beta*atanh( (1-exp(2*beta/Umax*theta*(uave-Umin)-log(cosh(beta))-beta*theta)) / (tanh(beta))  );

% thinc 函数
T = Umin + Umax/2*(1+theta*tanh(beta*((x-xa)/(xb-xa)-x_tilde_i)));

% 积分结果 int_T 应该等于 uave
int_T = simplify((subs(int(T,x),x,xb)-subs(int(T,x),x,xa))/(xb-xa));
double(simplify(int_T))

参考

1. ^张伦, 牟斌, 蒋浩, 等. 气液两相流混合模型代数重构方法[J]. 空气动力学学报, 2023, 41(8): 117−123. doi: 10.7638/kqdlxxb-2022.0107
2. ^Deng, X., et al. (2019). “A fifth-order shock capturing scheme with two-stage boundary variation diminishing algorithm.” Journal of Computational Physics 386: 323–349.