一个自由液面的瞬态问题：Level set还是动网格？

soliton

这两天遇到一个非常有趣的问题，似乎和前一段时间论坛上hgniang大侠的问题有些接近，写出来请教一下大家：

一个半径为R金属液放在一个平板上，初始为圆形。

这时候，在垂直于这个金属液的平面上施加一个低频交变的磁场，这个金属液圆盘会在电磁场的作用下流动、变形，其边界由表面张力和电磁力、流动联合控制。

可以推导磁场的变化方程及电磁力的表达式如图中所示.

我初步的思路：

计算域中包括内部的金属液和外部的空气，其界面就是要追踪的那个界面。流动就用N-S方程，界面用Level set表征，然后耦合进去上面的Poisson方程。N-S方程的体积力中加进F(T，t）写的函数。

但我有一个问题：就是在Level set中，怎么把每时每刻的那个边界条件T=0写进去？因为这个边界是不断变化的，并且在Comsol中，这个界面是内部的边界，且和解出来的函数Phi有关。

或者，计算域中只包含金属液，采用动网格？不够动网格的我不太熟，有没有高手推荐几个例子学习学习？

soliton

这个问题里面有些东西，似乎是个共性问题，比如hgniang大侠的帖子，一个气泡在周期性压力作用下的变形http://forum.simwe.com/thread-946336-1-2.html。

上面图中的照片，是有人做实验的结果（水银液滴，直径30mm）。这个过程会出来如此对称的结果，还有三次，六次的对称的，真的让人很惊讶！当然，还出来一些特别复杂的花样。他后面的方程，几乎就这么多了。如果能用Comsol求出来，那就太好了。

如果用Level set来求解这个问题，我目前遇到的最大困难，就是那个金属液和空气界面的那个边界条件，如何施加？在Level set 中，每次计算完，界面有个重新初始化的问题，如果能在重新初始化过程中加进去这个边界条件的约束，就好了。关键是，如何加呢？

请高手指教！

hgniang

这两天遇到一个非常有趣的问题，似乎和前一段时间论坛上hgniang大侠的问题有些接近，写出来请教一下大家：

一个半径为R金属液放在一个平板上，初始为圆形。

这时候，在垂直于这个金属液的平面上施加一个低频交 ...
soliton 发表于 2010-9-6 22:03

我不是很清楚你的边界条件是如何变化的？是已知的有规律的变化，还是在求解过程中未知的变化？后者的话我也不知道了，我也在摸索。

至于动网格应该比较简单吧，耦合上ale模块，确定边界的运动情况就好了。

我的气泡模型没有用到LEVEL SET，我只是用了Imcompressible N-S。

关注一下这个问题，希望高手进来解答

soliton

我不是很清楚你的边界条件是如何变化的？是已知的有规律的变化，还是在求解过程中未知的变化？后者的话我也不知道了，我也在摸索。

至于动网格应该比较简单吧，耦合上ale模块，确定边界的运动情况就好了。
...
hgniang 发表于 2010-9-7 09:40

谢谢指点！我的边界是在求解过程时刻变化的，和各参数有关。这两天在考虑动网格这个事情。

关于水平集法在两相界面上，边界条件施加的问题，我觉得是一个大家都需要解决的问题：我以前一直都把这个边界当做是内部边界，是系统默认的内部边界。现在看来，系统内部的边界是怎么回事，没有搞清楚。

期待高手来指点一下。

facklaoer

非常好玩的现象
虽然没有看懂楼主的问题，但是忍不住想做一下
两相流好久没做了

facklaoer

请问，T是什么？

soliton

谢谢大侠指点！T是计算域中的一个标量场。在液滴和周围空气接触的界面上满足T=0，在域中满足上面图中写的Poisson方程。

假设采用Level set 方法，大致的思路是：求解域包含中间的液滴和空气（外面）。

1、 先求解关于T的标准PDE方程，在边界上保持T=0。得到T的分布。

2、 根据Fx=-1/2*sigma*omiga*B^2*sin(2*omiga*t)*Tx,   Fy=-1/2*sigma*omiga*B^2*sin(2*omiga*t)*Ty,求解流动方程，由水平集法得到新的瞬态界面。（其中电导率sigma在金属和空气中差别非常大，于是空气中基本没有力的作用。）

3、 这时候，Level set 重新初始化这个界面。此时界面上的仍然必须维持T=0。在新的条件下，再求T,得到新的力。循环计算。

关键的是，在现在Comsol的框架下，第三步该怎么做啊？并且它的初始边界条件设好之后，后面就很难再对这个边界进行操作了？！

soliton

关于本问题的一些参数：

频率omiga=1.6Hz,
电导率sigma=1e6[S/m]
Bz=0.2T,
液滴直径30mm
金属的密度13600Kg/m^3,
粘度1e-7m^2/s,
和空气的界面张力0.485N/m.

facklaoer

我做过一些高电场下的气泡变形，因为气泡在界面处的场强及介电常数等发生剧变，而体积力是与其梯度有关的，所以这种体积力自然地加到了交界面上，引起了两相形态的变化。
还是没能理解T的含义，是不是有人用这些控制方程模拟成功过？如果不是，可以考虑重新研究物理过程与控制方程。

另：如果楼主想在交界面上始终保持T=0，可通过在整个域中设定T*（phi==0.5）=0来保证吧。

soliton

楼上高手啊！

关于标量场T，是在推导其中电磁场、电磁力的过程中引入的一个新的标量。我这个是同课题组另外一个人推导后给我的，我正在消化中。应该还没有人对这个问题模拟成功过。

“想在交界面上始终保持T=0，可通过在整个域中设定T*（phi=0.5）=0来保证吧。”
能不能给点明确的指点：如果是在comsol水平集的框架下，具体在什么地方实现这一设定？

感激不尽！

facklaoer

非常抱歉，具体的实施方法我也没有，只是一个思路。
但是，我觉得用这几个控制方程就能模拟出你想要的结果，很难。
特别是电磁力的控制方程，应该比这个复杂挺多。
个人意见，仅供参考。
祝好运，一起思考，一起提高。

soliton

非常抱歉，具体的实施方法我也没有，只是一个思路。
但是，我觉得用这几个控制方程就能模拟出你想要的结果，很难。
特别是电磁力的控制方程，应该比这个复杂挺多。
个人意见，仅供参考。
祝好运，一起思考，一起 ...
facklaoer 发表于 2010-9-8 20:28

还是非常感谢！最近在论坛上承蒙各路大侠指点，学到了非常多的东西。

我觉得，上面的这个问题，几乎满足一个好的、有趣的问题的所有条件：

1、 过程简单，现象却又异常复杂；
2、 初一思考，以为能做。仔细一想，有难度。

原始的文献在：Journal of Fluid Mechanism, 2005: 285-301. 谢谢！

abubob

soliton 你好，我的MSN这几天无法登陆，不知道我发个你的邮件是否收到。关于这个界面力如何施加的问题。
    大家好，关于这个界面上的力如何施加的情况，我和soliton讨论过模型，论坛上也有相关帖子（8月发的）。这个月我没有关注模型本身，主要关注的是这个level set 方法。看这个comsol水平集函数如何捕捉界面。我发现的comsol这个水平集的有些不一样，不知道理解是否正确：（1）phi函数应该充当了传统水平集中Heaviside函数；（2）表面张力是个表面力，其表达式里采用了Delta函数进行定义界面上的力，但有个量纲？传统中应该没这个量纲；（如果采用这个delta函数定义界面上的力，我试过，力的体积密度就会发生变化，耦合方程就会出错）；（3）水平集函数数值解的右边表达式，未能明确其采用的数值解思路，看不懂其用有限元求解的套路。

soliton

abubob: 欢迎讨论，很抱歉，似乎没收到你的邮件。

关于你提出的这几个问题，我的理解是这样的，不一定对：

1）phi函数本质上是计算域中每一点到界面的距离函数，每一个时间步长计算时，在每一个节点上都得到一个Phi值，然后规定Phi=0.5这些节点组成新的“界面”。重新初始化的过程中，针对这个界面，每个节点的phi值重新再定义一遍。
而Heaviside函数只是用来定义“界面”的交界面上的物性参数而已，它的自变量是Phi。这两个有本质不同。
在这个过程中，Leve set里面设置的过渡层厚度、以及重新初始化速度起了作用。（但关于重新初始化速度到底怎么起作用，我还是没想明白）

2）表面张力的这一项，在N-S方程中写进去的，里面包括了界面张力系数、曲率、Delta函数和法向的表达式。应该是一个非常复杂的计算过程。如果要自己写函数来算，不知道你是怎么实现的。从物理上讲，它应该是个应力的单位。不知道你说传统中没有量纲，是什么意思。

3）水平集函数数值解的右边表达式，我的理解就是一个光滑函数。直接定义为零是一个太强的“约束”，矩阵计算过程不容易收敛。但关于它内部的有限元求解思路，我就不甚了了了。

欢迎讨论，共同提高！

abubob

solito:欢迎讨论。同时欢迎大家一起讨论。 我下面对第一个问题，我把我的疑惑说出来。以下这段话是帮助中对phi函数的描述： the level set function, . In COMSOL Multiphysics, is a smooth step function that equals zero in a domain and one in the other. Across the interface, there is a smooth transition from zero to one. The interface is defined by the 0.5 isocontour, or the level set, 从它的表述，其定义与Heaviside函数基本一致，无论H函数的表达形式是什么。在对过渡层的定义时，level-set for two phase flow中，也是采用phi作为H函数的。因此，我对这里的phi函数理解的是定义为两个作用：H函数与传统phi的函数的作用？
     传统phi=0应该是界面，这里的0.5是怎么通过对phi函数的定义得出来的呢？这个我还没想清楚。
     欢迎大家讨论。

zqlgeo

1)COMSOL中phi函数本质上已经不是到界面的距离函数，因为水平集函数并不一定就非要是是距离函数，只不过距离函数是符合水平集函数特点的比较直观的一个函数,所以初期大家都采用距离函数.一切符合水平集特点的函数都可以用来作为水平集函数的。如现在comsol中采用平滑后的Heaviside函数
2）界面张力系数、曲率第一个是物理特性，第一个采用了平均曲率计算方法，利用phi梯度积分平均化，Delta的作用只是为了保证界面总的表面张力等于张力系数。因为对界面上对delta作积分是等于1，也就是说把表面张力作加杈分配在表面张力一定范围。
3）函数右力包函了两部分，一部分用于跟踪界面，就是传统的水平集推进，而一项用于重新初始化，将两项合并，则使推进过程中伴随初始化。使原来要分开计算的进行了方程的叠加，因为两相步骤的左项是完全相同的。这也是作者的高明处。

facklaoer

看来大家有点跑题了，我再来探讨下此问题解决的思路。
首先，使用的是pde系数型和相场法两相流 ，因为我觉得相场法比较稳定，水平集法也可这样参考。
设常数b1=1，b2=-1，标量表达式b=b1+(b2-b1)*Vf2_chns；很明显，b在一相中为1，另一相中为-1，交界面处自然有过渡的0。
在pde系数型求解域设定中，设c=b，a=abs(abs(b)-1)，f=b。这样，当在两相中时，控制方程为ΔT=-1，而在交界面处会出现T=0，计算时能进行一些，不过总是以不收敛收场。我想是由于b在交界面附近的平滑度较差所致。不过，这个思路比较接近楼主所要解决的问题了。
这种解法是个思路，但离问题的解决有多远尚未可知。任何问题都应努力寻求各种解决的思路，其过程是艰难而且复杂的。楼主所列的方程，或许动网格法更适合些。
还想说两句：由于时间有限，这个题目也就能按照楼主的方程解决到此了，毕竟这不是我的课题，不能分出更多的精力。但是，对于楼主要解决这种问题的思路还是持保留态度，因为我做过类似的问题（只涉及电场），当时解决的思路要比这个复杂，但我觉得更有效，这在9楼也说了些。
祝楼主好运，如果能够解决，别忘了晒晒图。

soliton

非常感谢楼上各位的讨论，学到了很多！

尤其感谢facklaoer 大侠多次极具建设性的指点。我会领会学习。再次感谢！

abubob

16# zqlgeo
    感谢zqlgeo 大侠的讨论，长了不少知识。我把它的这个算法基础，就是参考里2005年的那篇文章下载下来了，我再看看。
    关于如何在这个界面上施加力的情况，就是如何定义delta函数的问题，我正在找comsol 是否有逻辑表达式，这样的东西。我沿着这条路走走看，有结果了告诉大家。欢迎大家继续讨论。

soliton

看来大家有点跑题了，我再来探讨下此问题解决的思路。
首先，使用的是pde系数型和相场法两相流 ，因为我觉得相场法比较稳定，水平集法也可这样参考。
设常数b1=1，b2=-1，标量表达式b=b1+(b2-b1)*Vf2_chns；很明显 ...
facklaoer 发表于 2010-9-9 22:03

最近，液滴在纯粹表面张力下回复-变形的例子，尽管还不完全真实，但有解决的苗头。

开始试一下facklaoer大侠这个想法。

初步的结果显示有一点点问题：

在纯粹的两项中，phi=0和phi=1都没有问题。问题就是界面上的邻域（过渡层）上，b有值，但既不等于1也不等于-1，此时PDE方程中系数c/a/f都有值。

就这一层没有满足上面界面上T=0的条件，解出来数值差了好远。

试了直接赋值，效果也不好。还在思考中。