一个非常有趣却又奇怪的问题——matlab画的图倒入comsol作为计算域的问题

soliton

这个问题出现的物理背景很简单：一个椭圆（或者任意形状）的液滴，在表面张力的作用下会慢慢回复为完美的圆形。

针对这个问题，先试了水平集，后用了动网格(http://forum.simwe.com/thread-948772-1-3.html)，为了施加表面张力，不惜弄上了weak form!( http://forum.simwe.com/thread-955092-1-2.html)。

对于椭圆的情况，目前看来是比较完美地解决了（计算结果见http://home.simwe.com/space-5330-do-album-picid-33341.html）。

最近想把初始形状弄得复杂一点。比如说有几个角的.用上述算椭圆完全相同模型和边界设置，但却一步也算不下去。通常的错误是两种：（1） Failed to find consistent initial value……(2) Error: 6199, Attempt to evaluate real square root of negative number.,- Function: sqrt. 同样的模型，按道理说不应该这样呀。


思来想去，因为只是更换了计算域的图形，所以我怀疑这里出了问题。

对于这些复杂的形状，我是按照论坛上一个帖子的指点，先在matlab中写函数，plot线，然后通过comsol中的import→Geometry objects，再把线变成实体。我用matlab画图的方法做了“椭圆”的情况，同样也不能算。所以几乎可以肯定，就是图形的事情！

唯一的区别就是：由于采用matlab做那条边界上的线，比如说写了下面的函数：

a=0.03;b=0.002;mode=3;

sita=0:2*pi/100:2*pi;

for i=1:length(sita)

x(i)=(a+b*cos(mode*sita(i)))*cos(sita(i));

y(i)=(a+b*cos(mode*sita(i)))*sin(sita(i));

end

c=line1(x,y);

事实上，这个曲线C是含有100等分的“线段”组成的，所以在模型里有100个边界条件要设置！（而原来可以成功计算的“椭圆”是直接在comsol里画出来的）那么，在求表面张力的时候，两个线段相交的那个点上的曲率就可能是不连续的？导致采用弱解形式设置的表面张力不能计算（有冲突）。区别就在这里吗？

我也不是很确定，请教一下大家：从matlab里画的图导入Comsol就只能有这么多条“边界”吗。还是有什么办法，可以把一些比较特殊的，必须用函数写出的“特殊”形状“融合”起来啊？

soliton

再重复一下我的问题：
如果一个比较特殊的图形，要采用matlab的形式写函数，再往comsol里import，有没有可能避免最后的图形在微观上是“一条一条的线段”，而是一个处处可导的曲线？
（这个问题，对于处理与表面张力有关的问题来说，曲线处处连续可导，似乎还是比较重要的）

谢谢指点！

soliton

又想了想：到底comsol建模的时候，在draw中，选择画一个“圆”，它内部是怎样定义这个圆的呢？

而如果我们由外部“作图工具”画好后，import进comsol，它内部又是怎样处理这个“图形”的呢。

以我的经验，同样的一个椭圆，在comsol直接画就是一个数学上完美的“椭圆”，用matlab画好导入，就是一条条线段组成的“多边形”。为什么呢，为什么？

有没有其他的函数，可以代替这个line1?请高手救命！

soliton

翻到一个好早的一个帖子：http://forum.simwe.com/viewthrea ... amp;highlight=curve

但还是没有回答上面的问题。看来要研究一下各种画曲线函数的用法。有人对rational bezier curves 有心得吗？好像有点复杂。

soliton

顺便问一下，下面图中的length2,3,4，还有weights是什么意思啊？

为什么我用
a=0.03;b=0.002;mode=3;
sita=0:2*pi/100:2*pi;

for i=1:length(sita)
x(i)=(a+b*cos(mode*sita(i)))*cos(sita(i));
y(i)=(a+b*cos(mode*sita(i)))*sin(sita(i));
end

c=curve2(x,y);

得不到一条曲线呢？

看起来似乎是这个意思啊？！

soliton

杯具了，这个问题居然是一个计算机图形学的问题！有难度了。

facklaoer

朋友，我觉得，你的课题一开始就走入了误区
请重新审视COMSOL和有限元方法，有限元方法要比你做的这些美好的多
你需要仔细考虑控制方程
然后通过电磁场和相场法才能解决你的问题
不要徒劳无益的做无用功

soliton

朋友，我觉得，你的课题一开始就走入了误区
请重新审视COMSOL和有限元方法，有限元方法要比你做的这些美好的多
你需要仔细考虑控制方程
然后通过电磁场和相场法才能解决你的问题
不要徒劳无益的做无用功
facklaoer 发表于 2010-10-29 19:43

感谢你不懈地指点。我完全理解你的出发点。

关于控制方程的部分，尤其是电磁场的部分，我有一个比较“讨巧”的推导，在二维的情况下，那个简单的方程是个非常好的近似！我也试过二维平面感应电流下磁场强度的方程，两个过程的解是完全一样的（诚然，电磁场上，你认为应该仔细考量，那是十分自然的。）

如果说有误区，我认为误区在两点：
1、本质上，这是一个参数稳定性的问题。采用这种“数值试验”的方法，原则上可以告诉我们“是什么”，而还是不能告诉“为什么”。这是有限元模拟的局限。因为那个我不明就里的“初始扰动”是什么，完全不能穷举。但能够知道是什么，也聊胜于无。

2、从三维到二维的简化。尽管初步的分析也表明，可以做二维的近似。但毕竟，真实的世界是三维的。这个我更担心。

从水平集到ALE的初步尝试表面，在椭圆振荡上解的结果是一样的（无论有无磁场）。只不过水平集的算法本身收敛性较差一些，还有相泄露的问题，比较起来不如ALE好。但这只是数值方法的问题，无关物理本身。你所提到的相场法，我粗略看了一下。是否真的比水平集好，我不确定。

通过这个玩票，坦率说，学了很多：ALE，ALE下通过弱解形式施加表面张力，ALE下的REMESH，微分几何的一点点。论坛上好多人都给了宝贵的指点。即使课题貌似没有进展，也不能说是做了无用功。有限元方法的美妙，我只能看出万分之一。但，这大概就是我们上论坛来交流的原因吧。

再次感谢你的持续关注和指教！其实，我欠你一个贴图，希望有一天真的能贴上来。

followtime

我觉得楼主的问题很简单，关键是你没搞清点阵图和矢量图，

你用matlab通过设定点后给出的园只是一些点连接后的集合，而comsol内部存储的是一些点的数据关系，例如你到处comsol的脚本后发现，它的园有两种，中心和半径与脚点两种形式，所以comsol存储的不是那些点的集合，而是部分点的关系，

所以如果你想要把自己matlab底下的图形导入comsol，就必须以点的关系的形式给出，也就是按照comsol脚本中区域的定义的各种组合生成。

由存储的某些点的关系的生成图形的方式通常叫做矢量图方式。它与点阵图的最大区别在于点与点之间是数学的关系，所以放大或者缩小精度自动可以调整，但是点阵图存储的就是点的大小，颜色和点之间的距离和角度，放大后往往精度损失，缩小后也会差生很多问题。

另外，如果按照楼主所示把圆剖分100份导入comsol后生成网格，那么无论在怎么加密，边界的逼近程度就是固定的了。如果按照comsol存储的点的关系生成网格，随着加密，边界可以无限的逼近精确的圆。

soliton

谢谢9#指点，隐约是觉察到了你说的这个意思。
一开始是想自己写程序来解决的，但似乎搞不定。转而想，comsol一定有它自己内部的函数干这个，要不圆、椭圆什么的它怎么处理的？
看了一下说明书，从reference guide里面找了个命令geomspline，就是通过多个点的样条曲线。似乎是解决这个问题了，据说，点和点之间的连接是C2的了。现在在试它的效果。

piggycan

10# soliton
朋友你好！我本来想说用geomspline命令的，看来你已经发现了。同样的，变成实体的工作也可以在命令行模式解决。geomspline命令含有多个参量，是插值样条多项式的次数问题，对本质影响不大。

我曾经做过类似的工作，不过是半导体仿真的。这个领域还将出现一个问题就是当把模型导入到AUTOCAD制作掩模图形后你的边界会变成很多分段，也是一个难题!

piggycan

10# soliton
朋友你好！我本来想说用geomspline命令的，看来你已经发现了。同样的，变成实体的工作也可以在命令行模式解决。geomspline命令含有多个参量，是插值样条多项式的次数问题，对本质影响不大。

我曾经做过类似的工作，不过是半导体仿真的。这个领域还将出现一个问题就是当把模型导入到AUTOCAD制作掩模图形后你的边界会变成很多分段，也是一个难题!

soliton

楼上朋友：谢谢指点。

其实我并不完全清楚geomspline命令的算法。我想你那个问题，“模型导入到AUTOCAD制作掩模图形后你的边界会变成很多分段”，如果能够找到算法，做一个AUTOCAD的“接口”，不知能不能解决。当然，谈何容易？！

piggycan

13# soliton 在今年年会上听说未来将出现与AUTOCAD的无缝接口！现在我只能用一些妥协的方案。

关于geomspline的详细说明，可以在软件自带的参考文档里找到。或者更简单的是，在matlab中使用help命令，一目了然。