【读论文(paper reading)】Hydrophobic and Hydrophilic Solid-Fluid Interaction

亲-疏水表面流固耦合(SIG2022)

阅读时间：2022年2月上旬

Hydrophobic: 梳水
Hydrophilic: 亲水

从左到右的几个场景：

• 纸因为浸润自然卷曲包住液滴，
• 较重的软体被强表面张力撑起来，
• 液体从疏水叶子上滚落，
• 不同浸润性流体的接触角

论文地址：

https://jinyuan-liu.github.io/projects/2022hydro/

Video:
https://o9ixctz0o7.feishu.cn/file/boxcnUXKeQSGw4yQXxAJbDhYGFh

https://youtu.be/X_nPvKjl5k0

Goal是什么

模拟新现象：软体和流体的流固耦合，重点在表面张力。包括以下几种小尺度下的特殊现象：

• 弹性毛细管现象（就是液滴滴到小纸片上让纸片卷起来，封面图1）
• 莲花效应（液滴在叶子上滚落, 封面图3）
• 不同材料浸润角度（封面图4）

由于有双向流固耦合，所以第一个场景那个纸能够受力自动卷起来。这在之前是做不到的。

模拟的一个核心是将contact line（接触线）的动态演化模拟出来。

背景知识：关键词和术语解释

1. 表面张力

(ref:wiki)

定义：表面张力（英语：Surface Tension）在物理上的狭义的定义是指液体试图获得最小表面位能的倾向；广义地说，所有两种不同物态的物质之间界面上的张力被称为表面张力。表面张力的单位是$\frac{N}{m}$或$\frac{J}{m^2}$亦即，单位长度的力或单位面积的能。

直观大小：水的表面张力约为72.75毫牛每米。水分子由于是极性分子，有范德瓦尔力，因此明显高于同密度的其他物质。例如乙醇只有22.55毫牛每米。因此乙醇可以做表面活性剂（大部分化合物都可以）。且随着温度升高，表面张力会降低。

通俗解释：

• 物理上说，分为力的角度和能量角度。

  • 力的角度：在液体内部，每个分子都在每个方向都受到邻近分子的吸引力（也包括排斥力），因此，液体内部分子受到的分子力合力为零。然而，在液体与气体的分界面上的液体分子在各个方向受到的引力是不均衡的（见图解），造成表面层中的分子受到指向液体内部的吸引力，并且有一些分子被“拉”到液体内部。
    
  • 从能量的角度：液体内部分子周围有大量分子，因此，内部分子的分子势能较低。然而，表面层中的分子周围的分子明显小于液体内部分子的，所以，表面层的分子有较高的分子势能。为了达到低能量的稳定状态，表面层中的分子有向液体内部移动的趋势，从而导致表面层中的分子数量减少，宏观表现为液体表面积减小。

• 数学上说，就是找到最小面积。或者最小能量。可以用吉布斯自由能G对面积A的偏导数得到表面张力。（T,p代表特定的压力和温度下）
  

Bridson书中中给出浅显解释表面张力的原理：

相比于空气分子，表面水分子会被内部水分子拽回去

倾向于暴露给其他物质的面积最小化

最简单的表面张力模型就是0.073*接触面积

2. small-scale elastocapillary phenomena：小尺度的弹性毛细管现象

定义：引自Wiki: Elasto-capillarity is the ability of capillary force to deform an elastic material. 弹性毛细管现象就是毛细作用力使一个弹性的材料发生变形的现象。

这种力在宏观上 体现不出来，但是微观上可以其主导。一种非常经典的实验就是毛细折纸(Capillary origami)。 就是第一个场景所展示的。

https://en.wikipedia.org/wiki/File:PDMS_folding.ogv

这是弹性壁与液滴之间双向耦合才能做到。

杨-拉普拉斯公式可以描述毛细压力：

3. 接触线 contact line

接触线。就是气液固的三相交汇处。如下图所示。

绿色的是三相接触点。在二维是点，三维是曲线。因此它就是接触线。实际上他应该是一个小的控制体。如图中虚线区域所示。它的位置由体积应力力和表面应力共同控制。

在这个小的控制体上，我们可以列出质量和动量守恒。

该公式中，由于空气的密度很小，因此rho_a可以认为是0。同理sigma_a也可以忽略。

其中

• sigma是体积应力
  • sigma_s是固体弹性力
  • sigma_l是流体压力
  • sigma_a忽略
• Gamma代表的是表面应力。
• 下标a表示air, s表示solid，l表示liquid。

浸润与接触角与亲疏水

接触角小于90°，就是浸润，或者说亲水
接触角大于90°，就是不浸润，或者说疏水

疏水有个著名的现象：荷叶的自洁现象。

和刚性壁之间的接触角，可以采用Young-Dupre定理计算。

或者写成

SOTA是什么以及当前问题

Ruan 2021: 模拟了三向耦合的表面张力

Fei 2017: 模拟了浸润的头发丝

• 以往的表面张力驱动的，大都是刚性壁。因此是无法卷曲起来的。
• 以往的工作，接触线是静态的，无法动态演化。因此无法模拟弹性毛细管。

作者特别与阮良旺的工作进行了对比

https://dl.acm.org/doi/10.1145/3450626.3459862

Liangwang Ruan, Jinyuan Liu, Bo Zhu, Shinjiro Sueda, Bin Wang, and Baoquan Chen.
2021. Solid-! uid interaction with surface-tension-dominant contact. ACM Transactions on Graphics (TOG) 40, 4 (2021), 1–12.

作者claim:

• ours可以模拟弹性壁，而Ruan只能模拟刚性壁
• ours可以模拟亲水，而Ruan只能模拟疏水
  

本文的Solution（实现）

离散化

提示：grid指的是背景网格。而mesh指的是固体的网格，包括surface mesh和tetrahedron mesh(volumetric mesh)。

流体

流体使用MAC，然后采用advection-projection的方式求解。advection采用semi-Lagrangian。忽略粘度。

使用level set来追踪液体表面。表面张力使用半隐式积分。使交界面速度的sharp change diffuse掉。在narrow band的速度上求解Poisson方程

其中delta是在narrow band epsilon上定义的函数。本文中epsilon取3倍dx

计算其他交界面的时候，有两种表面张力：

1. Y_la（如图蓝色菱形）
2. Y_ls（如图黄色三角）

（Y_sa被忽略）

在计算接触线的时候，则三种表面张力都要考虑。

固体

固体使用volumetric mesh来离散。3D中就是四面体。该mesh比背景网格密1.5倍。

采用Bridson2002提出的半隐式积分。也就是将位移驱动的力（弹性力）和速度驱动的力（阻尼）分开处理。

实际使用的是简化版本

Ds代表隐式积分力的系数矩阵，Is是显式积分的力的系数矩阵

力的微分参照Sifakis2012 Course

接触线

接触线几何

接触线的几何离散包括三个要点：

1. 余维1 level set的生成
2. level set转换mesh
3. mass lump

核心要点是先在固体mesh上生成level set，然后转换为拉格朗日mesh。

然后使用这个余维1的level set来计算余维2的单纯形网格。在三维上就是线段。

在生成完拉格朗日表示之后，我们计算每个mesh vertex上的质量。该质量是通过周围grid中的固体和液体质量lump出来的。

以上这三步具体来说如下：

1 生成solid surface mesh上的level set（也就是接触线）

首先从固体表面mesh上生成level set: phi_c。这是液体和固体之间的交界线。

这是根据余维0的液体level-set和余维1的固体表面网格的交界来计算的。如下图所示

在生成这个level set的时候，采用了一种自创的enlarge and shrink算法。

这套算法要处理的就是如何准确捕捉contact line的level set。

输入：

• 液体的level set
• 固体的表面网格

输出：

• 固体表面上的接触线

先尝试扩大液体的level set 一个小量，
然后找到扩大了的液体level set和固体表面mesh的交界，也就是生成固体的表面上的leve set。
然后仅在固体表面的level set缩小epsilon。

epsilon, 取dx。
epsilon取值越大，生成的contour越 smooth。epsilon越小，生成的contact line越精确但contour 越bumpy。

2 重新初始化（优化一下level set)

使用fast marching算法重新初始化 solid mesh的level set。

fast marching的意思就是维护一个优先队列（可以是最小堆），从一直的grid point出发，计算他们离surface的距离。这个最小堆中存储的是未知的grid point，其键是他们的当前距离。首先以已知grid point的neighbors初始化这个堆。然后从这个堆中取出最小值，设为已知点。更新其未知邻居的距离和最近点信息。不断重复直到该优先队列为空。对于n个未知网格其复杂度为O(nlgn)。

具体细节可以看Zhu2015

然后求解了Eikonal 方程（如下）来求出曲面上的测地距离。

实际求解中会转化为这个方程。其中tau是虚拟的时间。也就是在虚拟时间上行进（也就是迭代），使得grad phi的绝对值趋于1。观察可发现当|grad phi| =1的时候，phi是不会变化的。

3 拉格朗日表示（就是接触线的粒子化）

使用marching triangle算法（Lorensen and Cline 1987）来获得拉格朗日表示。算法基于以下原理：当一条边的两个入射顶点的符号相反的时候，这个点就是接触线上的点。如图所示

然后进行一下拓扑修补（ [Brochu and
Bridson 2009]）就是把过长的线段细分下，过短的线段给合并下。如图所示。保证线段长度在0.5-1.5倍原始网格的平均边长内。

4 mass lumping

由于每个contact point都是有控制体的（列控制方程时要计算）。因此也会有质量。contact point的体积是通过以其为中心的一个球计算的。在球内对不同相的cut-cell进行加和得到其体积。

接触线物理

根据动量方程

体积应力：

• 空气体积应力为0
• 流体体积应力为压力
• 固体体积应力为弹性力

表面应力：Y*t
Y是个系数。
怎么找t?

怎么找切向

表面应力的方向是沿着切向的（下图的红箭头）

输入：交界面上sample points
输出：交界面的切向

做SVD分解

有了表面力的切向t_ij，然后表面张力系数Y_ij是参数给定的。界面力就是Y_ij t_ij

流固耦合

无接触线的流固耦合

流固
Two-way coupling of fluids to rigid and deformable solids and shells

通过SDF来判断grid的类型（固液气）。

grid指的是背景网格，mesh指的是固体上的动网格。两种网格同时存在（dual cell）。

使用插值矩阵Ws来在grid和mesh之间插值。具体来讲，是从固体mesh插值到流体的grid。被插值的量是速度和力。比如速度插值公式为：

Ws中每行存的是交叉的面积。 交叉面积是通过先在三角形内撒点，然后计算在dual cell内的比例来计算的。最后还要在每行normalize一下。

(个人猜想：Ws应该是个正交阵，所以W^T=W的逆。求逆相当于从grid插值回mesh，只要转置就行了。normalize其实为了保证是单位正交阵。）

在液体的pressure projection，增加solid的速度项。

其中

是压力对应的冲量（单位体积的）

动量方程

组合起来

其中V是dual cell的体积。Ws是从固体mesh到grid的插值矩阵。vs是积分完所有显式力（也就是弹性力）的中间速度，u是做project前的速度（也就是非压力梯度力作用后的速度）。Ms是固体的质量对角阵。M是液体质量对角阵，可以通过level set求出。

Ms~是把液体质量也叠加到固体上的质量矩阵。

这就是流固双向耦合。

下面讲的是有接触线的流固双向耦合。

有接触线的流固耦合

首先，选取的控制体区域R=3 dx

然后同样设定一个插值矩阵Wc。Wc用来在contact line和grid face之间插值。

然后传递动量。

其中下表c代表contact line

Wc也是插值矩阵。是从contact line插值到grid的。

是在接触线附近的控制体内的质量。

vc*是积分完界面力的接触线的中间速度。

结果

结果如何评价

很棒！

Amazing!

Limitations:

• 目前只有软体模型，而没有面或者线的模型。即使是叶子也是有厚度的。

• 并且使用的是显式积分。因此表中都限制了CFL数<1

• 性能没有做优化。用的Jacobi preconditioned CG，非matrix free。显存限制了模拟规模。并且Eq11和Eq15明明是非正定的，但仍然使用CG…不一定收敛