在玩中理解科研：一个椭圆液滴在表面张力作用下回复变形，Level set中表面张力的作用

soliton

前面发了个帖子，向大家请教一个电磁场作用下金属液滴变形的问题，承蒙众多高手指点，对问题的认识增进了不少（http://forum.simwe.com/thread-948105-1-1.html）。在这里特别感谢hgniang,fachlaoer,abubob,zglgeo和版主TBE的持续指点！

昨天按照fachlaoer大侠指点的方向，尝试将Level Set和系数型PDE耦合起来。（通过PDE系数的的巧妙设置，将变化界面的一个“动态”边界条件引入，fachnaoer大侠真神来之笔啊）

模型似乎可以算，但和fachnaoer大侠一样，计算的收敛总是有问题。后来分析，如果是一个完完全全对称的液滴，在我们的那个物理过程中（只有电磁场产生的力和界面的表面张力）金属液内部受到的“力”会是完全对称的，可能都不会产生运动啊！一定是在初始情况下液滴有一个什么“扰动”，才会导致它动起来。（似乎有点道理啊）

后来把这个问题分解了一下。其实那个PDE方程，单独解起来，完全没有问题。似乎收敛的问题出在Level Set这边。

于是试啊试，出来了这个帖子的主题：什么因素都不考虑，只考虑表面张力对液滴形状的影响，初始状态是一个椭圆，由于各处的曲率不一样，比较大的地方表面张力比别的地方大，就会逐渐回复为圆形。

这应该是一个简单的物理过程，算起来应该没有什么问题啊。可是老不收敛！（这里，可能涉及到abubob在上一个帖子里说的，表面张力怎么在方程中起作用的问题。我在前面的计算中，也发现，在模型库中那个Rising Bubble的例子，界面张力sigma就非常小，只有0.02。如果调大后，就总是不收敛。）

这是我算的一个初步结果：确实椭圆的液滴在回复，流动状态似乎也是对的。计算中表面张力从sigma=0.0485开始的，而实际的界面张力大概是10倍，甚至更大。一个初步的感觉，似乎网格对计算过程的影响非常大，但也不是越细越好。

附件是做的这个简单模型，供有兴趣的人玩玩。所谓悬赏，不过是博大家一笑而已，主要是请教一下大家界面张力大了之后，计算要注意哪些问题？

约定：表面张力从0.02一直到1.5，以sigma最大，最先算出收敛结果，说明计算过程者为胜！

facklaoer

真是奇怪，我把过去做的气泡模型表面张力调到1.5，计算一点问题也没有。
不知问题从何而来。

soliton

这确实是我没有完全理解的地方，想请教一下大家。

我上面的这个例子，Sigma大之后确实就不收敛。程序的收敛性可能是个比较复杂的问题，但我总觉得这里有一些可以通过简单设置来解决的“诀窍”，我一直没有掌握。

我还做过一个例子，界面张力大之后，能够计算，但出来的结果明显和物理事实不符，一个气泡上浮过程中，却慢慢往左边漂过去（对称的计算域，气泡处于中央）。在界面张力小的时候，还能符合直觉。有高手跟我说，可能是网格不对称，导致气泡左右两边的界面张力不平衡在计算结果上表现了出来。所以，一直对Comsol中的界面张力那一项怎么处理的，理解不是很透彻。

于是就想到借这个简单的例子，请大家讨论，指点一下。

facklaoer

说实话，这种小问题在每个成功的模型做出前，至少会碰到10个8个的
如果每个都这么讨论的话，有什么意思么？
况且，这种东西别人即使做出来有对别人有什么意义么？
别人不可能每天都围着你的课题转的。
最好提出些具有普遍意义的问题，比如：某个命令，某种处理方法之类的
不好意思，实在忍不住想说。

soliton

说实话，这种小问题在每个成功的模型做出前，至少会碰到10个8个的
如果每个都这么讨论的话，有什么意思么？
况且，这种东西别人即使做出来有对别人有什么意义么？
别人不可能每天都围着你的课题转的。
最好提出些 ...
facklaoer 发表于 2010-9-12 10:25

先抱歉让你有这样的感觉。

恰恰是我觉得这里面有一些共性的问题存在，才想到以这样的方式请教大家。在两三个人的小范围里我们确实讨论过这个事情，但由于我们接触Comsol的时间不算长，并且Level set 确实也是个以前没有接触过理论，在认识上肯定有所欠缺。

我非常感谢你在上个帖子中的“神来之笔”，对你可能是举手之劳，但对我这种水平低的人，那就是切实地提高。但我真的没有让大家“每天围着我的课题转”的意思。可以看出，我上面的那个模型，已经尽可能地将问题简化到一个所有人都理解的简单问题上了。并且，一开始，我就说了，这个只是供大家“玩玩”的。

论坛的存在，主旨就是让大家自觉讨论，共同提高。我觉得，这个帖子所提出的游戏，基本也还契合这个宗旨吧？

再次谢谢你！包括问题本身，我觉得一切都可以心平气和地讨论。

水若无痕

1# soliton
1“金属液内部受到的“力”会是完全对称的，可能都不会产生运动啊！”我觉得这个计算必须考虑重力，所以即使张力和电磁力完全对称，还有重力作用；
2 level set 对表面张力的模拟据说容易出现不守恒现象，不知道cosmol里面的level set 如何；
3 另外对于vof算法来说，表面张力系数大小不会影响收敛性的，但是融化的金属表面张力系数非常大，所以是不是达到了和电磁力同样量级导致收敛性不好？

水若无痕

1# soliton
1“金属液内部受到的“力”会是完全对称的，可能都不会产生运动啊！”我觉得这个计算必须考虑重力，所以即使张力和电磁力完全对称，还有重力作用；
2 level set 对表面张力的模拟据说容易出现不守恒现象，不知道cosmol里面的level set 如何；
3 另外对于vof算法来说，表面张力系数大小不会影响收敛性的，但是融化的金属表面张力系数非常大，所以是不是达到了和电磁力同样量级导致收敛性不好？

soliton

1# soliton  
1“金属液内部受到的“力”会是完全对称的，可能都不会产生运动啊！”我觉得这个计算必须考虑重力，所以即使张力和电磁力完全对称，还有重力作用；
2 level set 对表面张力的模拟据说容易出现不守恒 ...
水若无痕 发表于 2010-9-12 21:46

水若无痕版主：

谢谢！

1、我上面的这个例子暂时还没有电磁力进来。因为考虑的是液滴在一个平面上的变形，变为二维模型了之后，重力是不起作用的。所以，没有重力似乎也说得过去。

2、 这一点，我也只是怀疑。目前似乎是设置大一点之后，要么收敛有问题，要不出来的结果不符合直觉。但Facklaoer大侠，说他的例子没有什么问题。这也是我写这个帖子，想要寻求的解答。

3、熔化的金属表面张力确实非常大，但它和电磁力之间的相互关系，谁大谁小，如何相互起作用，我还真没有度量过，直觉认为“相当”（这是我要好好考虑的）。同时，我猜，是因为电磁力和界面张力在界面上有个共同的特点，就是几乎是突变的。如何寻找一个有效的光滑函数，大概会是一个问题。不知道大家在这个上面有什么经验？

水若无痕

9# soliton
请问：
如果是2d模型的话 怎么会变形呢？
这个变形是由于表面张力和液体浸润角共同作用吧，就好比在桌面上的一个水滴或者叶子上的一个露珠，即使不考虑桌子法向的重力，也应该是3d 建模吧？

另外我不懂电磁学，请问电磁力方向和作用效果如何？
我想用fluent试试。

*午夜流星*

如果你的计算模型与你上传的模型设置是一样的话，表明你连LEVEL SET方法的用法都没有弄清楚。计算之前还是先看看算例以及说明。
所以谈论表面张力根本毫无意思。

soliton

10# 水若无痕

其实，我想要的效果是这样的：http://math.berkeley.edu/~sethian/2006/Applications/Geometry/（水平集开拓者Sethian做的一个二维复杂连通图案在表面张力作用下变形为圆形的一个视频： Moviecurvcollapse。）

我猜测，初始状态为横向的椭圆，在表面张力的作用下慢慢回复为圆形，但在内部流体的流动惯性作用下，还会变形为竖向的椭圆。如此循环往复，知道流体的内摩擦力将初始的表面能耗散掉为止。这是理想状态了。无奈，程序算不收敛啊。

这两天毫无头绪地调参数，确实有一点*午夜流星*版主所说的，对Level set 的用法似乎隔了一层，没有理解透彻。



按照你的建议，准备好好去看看说明书先。还望午版不吝指点啊！

水若无痕

11# soliton
水平集和level set 我一点不懂，但是这两个是不同的方法，你的问题到底想用level 还是水平集啊？

soliton

12# 水若无痕

你是想说VOF和Level set吧？我这个问题里面，想用Comsol 的Level set方法做。

Fluent里面用VOF可能更方便，可惜我不会啊。

水若无痕

13# soliton
恩 我觉得vof 做这种液滴没问题的。

soliton

确实是没有什么明确思路地试了几天。一个朋友建议：在外侧空气区域的外边界上，不宜采用“壁面”边界条件。可能有道理，这样把流场压缩在一个范围内，是容易出问题。 随后，将外侧四个边界设为outlet, pressure, no viscous stress, p0=0, 一直算到0.76s左右，变形几乎为圆形了。但最后还是不收敛。这个计算非常快，在我机器上大概不到一个小时。图片是目前得到的界面（红色线）和流场，但接下来就不收敛了。 但这是在Sigma=0.048的时候，如果稍微调整一下界面张力大小，还是不行。 10楼“午夜流星”版主似乎知道问题出在哪里，还请明示啊！不胜感激！

soliton

昨天又做了一些尝试，有一小点进展：
在界面张力sigma=0.048(比真实的小十倍)时，算到2.2秒，不收敛。确实是：横向椭圆——竖向椭圆——再回复。这本是个非常简单的物理问题，没想到模拟起来真是不容易啊。还请大家指点！

反复观察计算出来的结果，发现都是由于外围空气中的流速突然变大导致的（这一部分事实上并不关心），我这样理解：由于空气的密度和粘度都非常小，只要界面上有一点点力的作用，附近空气的速度就会突然变大。由于网格等的问题，可能有时候界面是凹进去的，这样在凹陷下去的空气部分会出现与“主流场”方向完全相反的流动，搞得程序无所适从而报“不收敛”的错。

所以参数循此思路做了调整：将周围空气的“流动性”减弱（怎么和目前宏观经济有关啊？）。计算能进行，但牺牲了一些真实性。

可惜的是，如果换到真实的表面张力系数，还是不行。这回大概是网格的问题了。

soliton

由于那个动画超过附件的限制，我把它放在我的空间中：

http://home.simwe.com/space-5330-do-album-picid-15994.html。请注意左下角在最后时刻速度突然变大，并结束计算的过程。这就是问题所在啊！

soliton

这个帖子居然已经过去一个月了！时间过得真快啊。

当初在玩这个椭圆液滴的问题时，我面临两种选择：水平集和动网格，选择水平集方法入手的原因在于：
1、 水平集处理表面张力这一项有优势，不用去求动态的曲率，由comsol直接去处理了。
2、 同时，动网格在这两个工具中我想对不熟一些。

水平集搞了好几个星期，由于计算的收敛性问题，网格必须划得比较细（因为界面处必须细），但还总是出现在界面某点速度突然变化的这种意外情况（似乎是有些文献中说的“虚假”流动），且计算非常非常慢，通常这个30mm半径的气泡，DOF要到50万左右，这几乎到了单机的极限了！

上周开始，尝试ALE来做这个问题。困难重重，最难的还是表面张力如何加的问题。咨询了论坛大侠hgning，尝试把边界“化整为零”，计算每一段弧的曲率，并施加表面张力。他这个办法相当好，对于二维的问题，简单明了。但在我这里却没有成功，我分析了一下，大概是把我这个椭圆周长分为100等分还是不够细，能算，但结果不理想。

后来，我生搬了一套边界弱约束的理论，尝试在边界上施加一个弱约束（资料见两天前的帖子）：结果出乎意料的好！计算的动画我放在这里：http://home.simwe.com/space-5330-do-album-picid-33341.html

这个过程可以一直计算下去，直到液滴的振荡衰减到完全静止。且初期振荡的频率和理论值几乎完全一样。这个结果不能说错！

soliton

接下来要请教各位高手了，尤其是对弱解及边界弱约束有体会的大侠！

问题是：对于ALE下面，一个运动的界面处边界条件如何施加？

问了个老外，他提供了两个资料，具体位置在这里：http://forum.simwe.com/thread-953717-1-1.html 经过一番在边界上的表面张力乘以“试函数”的积分推导后，那个PDF文档说，表面张力象下面这张图这样写到weak里面去！

我的问题是：这个weak有10行，每行都代表什么意思啊？为什么要写在5/6行？为什么？到那里去获得关于这方面的东西。看了一下说明书，没找到。

soliton

从单纯的玩的角度来说，这个液滴振荡的问题是解决了。我试过的两种方案中，动网格+边界上表面张力的弱约束明显占优，计算速度比level set快20倍，且不容易有收敛性问题！在处理这个问题上，level set在界面附近“过渡”，并通过狄拉克函数处理表面张力（事实上当成了体积力）的方法，有些“和稀泥”的味道，不够好。

但ALE中，也不是没有问题。比如19#提出的那个边界弱约束的写法，到底是怎么实现的？我觉得，理解了这个问题，有助于解决一系列和“自由界面”表面张力相关的问题。

论坛里面除了几位斑竹，高手众多。不知道有没有对这个问题稍微感兴趣的同志？一起努力了！