发一个二维Euler方程DG程序 曲边四边形单元

glandetian

自己写了一个二维Euler方程的间断有限元程序

上次发了一个三角形单元的程序 因为不是曲边单元 所以在圆柱后面容易形成涡
现在把程序改为曲边四边形单元了 没有涡出现

单元是8节点四边形单元 节点编号顺序是 1 5 2 6 3 7 4 8 也就是四个角点依次是1 2 3 4  然后是边的中点编号 5 6 7 8.  

时间推进采用 Runge-Kutta 方法  

数值通量采用全局Lax-Friedrichs通量

仍然不能捕捉激波 因为没有做重构或者加人工粘性 等这个做出来了 再发一次。

程序没有进行优化 比如说内存的消耗没有优化  比如直边单元的边界积分仍然采用了曲边的积分方法 增加了计算量  比如面积分、线积分都是采用的是Gauss-Legendre-Lobatto积分 积分精度会比一般的Gauss-Legendre积分精度低一阶  等等问题。  二维的 纯属交流性质 就没有考虑这些问题  ^_^  

如果物面全部是直边 那么只要改变一个参数N 就可以获得不同的计算精度  且具有谱精度 因为单元的节点是Gauss-Legendre-Lobatto积分点。 其实就是谱元法  （物面是曲边的情况我不清楚是不是也可以通过提高基函数的阶数 也就是增加N 来提高计算精度）

欢迎交流！

zqb138

有没有算过亚音速内流问题啊，对出入流的数值边界处理碰到过没有？

Gerald135

顶了！

glandetian

加入口、出口边界条件就可以用

zqb138

你这个应该是超音速入流问题，我说的是亚音速入流问题，一般亚音速问题在出入流边界处由于特征关系只能给出部分物理边界条件，还有一部分变量是属于数值边界条件不用给出，但是间断有限元在边界处计算数值通量要用到边界两侧的所有物理量状态，那缺少的数值边界条件就需要采取某些办法得到才行的。这种情况的计算问题你用DG算过没有呢？

glandetian

这个我还不会   

我的这个程序 是计算亚音速圆柱绕流的 来流Ma=0.38
在物面用的镜像对称边界条件
在远离物面的地方 用的来流条件

[ 本帖最后由 glandetian 于 2012-4-13 08:15 编辑 ]

zqb138

哦，你的出入流是远场

[ 本帖最后由 zqb138 于 2012-4-13 07:07 编辑 ]

slaon

请问LZ，能否列举编写的参考文献啊？
新手，不太懂。但是很感兴趣。

glandetian

说一点自己不成熟的建议吧。
       首先要懂基本的有限元知识，特别是要明白加权余量法是怎么回事。我的理解是，有限元就是用的加权余量法，只不过形函数和权函数都是些紧支撑的函数。最好能自己编一个小的有限元程序，比如说求解一个二维的Laplace方程 （势流问题） ， 可以看章本照老师的书 《流体力学中的有限元方法》。会编这个程序了，就懂了有限元的基本概念了。
    然后，就可以学习时间相关的偏微分方程的求解。主要是两方面，一是在时间方向上，通常采用差分法，或者用得最多的Runge-Kutta法。最好编一个简单的用Runge-Kutta法求常微分方程初值问题的程序。这个很简单，但是很必要。可以自己想一个简单的常微分方程，有解析解的，编程计算，并与解析解对比一下计算结果。编这个程序时，可以参考一般的程序集的书，一般都有这个程序的。比如《常用算法程序集（第四版）》（徐士良）。二是空间上，就是有限元离散，针对不同的问题，会有不同的有限元方法，如CBS ，Taylor-Galerkin等方法，可以看O.C.Zienkiewicz的《有限元方法 第3卷 流体动力学（第5版）》清华大学出版社。
    上面两个程序自己编了后，会有些成就感，呵呵。。。然后可以学用间断有限元方法解Euler方程了。
    学习间断有限元时，要理解它与常规有限元的不同的地方，也就是单元间是间断的，所以不需要做质量矩阵的集成。另外，要理解什么是数值通量。关于数值通量的意思，可以看讲差分法的计算流体方面的书，比如《计算流体力学教程》 张德良 高等教育出版社。具体编程时可以看一些论文，比如《流体力学方程的间断有限元方法》（蔚喜军）。这篇文章用的捕捉激波的方法已经比较落伍了，可以看英文的用WENO和Hermite WENO限制器捕捉激波的论文。这方面的论文自己搜索吧。主要是CHI-WANG SHU舒其望和Jianxian Qiu丘建贤的论文，还有Hong Luo发在JCP上的论文《A Hermite WENO-based limiter for discontinuous Galerkin method on unstructured grids》。（这个我现在正在看，呵呵）
    对于边界条件，物面处好办，用反射边界条件就好。关于反射边界条件，可以看文献《High-order accurate implementation of solid wall boundary conditions in curved geometries》的第13个式子。远场边界条件，分入口和出口，比较麻烦，我也不懂，呵呵。。。前几天发了个帖子，也没人回答。书上说，对于亚音速情况，采用一维流动的特征关系来处理，对于超音速情况，入口用来流参数，出口用外推的方法。所谓一维流动的特征关系，可以看傅德薰老师的书，《计算空气动力学》《流体力学数值模拟》。我也没完全看懂。。。刚开始编程时，比如圆柱绕流，圆柱直径为1，远场用直径为20的圆，在这个圆上用远场条件，也可以算得比较好的结果。我的程序就是这么办的。
    说了这么多，也不知道对你有没有用。总之是希望这个论坛的气氛越来越好，有人提问，就有人回答。
    可能有说得不对的地方，请指正！

[ 本帖最后由 glandetian 于 2012-4-24 20:09 编辑 ]

slaon

哈哈，谢谢LZ。受益匪浅。

以前捣腾过fluent，也不知道CFD的核心是什么。现在课题要做可压缩流动的FSI，选择FEM计算CFD。看了一些文献，才发现激波捕捉是个大难题。
其实我一直想用连续FEM。我搜索了一段时间，发现基于FEM的Flux Corrected Transport方法可以在连续FEM中求解间断问题。

谢谢你的启发！

glandetian

不客气     CBS+人工粘性也可以捕捉激波 CBS也是连续的

[ 本帖最后由 glandetian 于 2012-4-24 20:07 编辑 ]

zqb138

楼主所说的那个边界条件，其实对于外流来说，入流一般给定所有物理量，出流按自由边界处理即可，而对于内流来说，才会去考虑什么特征关系什么的。

zqb138

还有楼主的Hermite WENO限制器看的怎么样了？我没有看过，现在计算的超音速流总是出现负密度和负压问题，而限制器如果没有保正性就会出问题，所以用舒其望的TVB限制器有时还要再加些特殊处理来保正性，不知道这个WENO限制器是否同时有保正性？

zqb138

总觉得这个限制器的计算量也越来越大了，重构过程走向越来越复杂。

glandetian

自由边界是什么意思啊？如果计算域是矩形，左右上下四条边界 物体在里面。来流从左边进入，给定全部来流条件，那上下和右边界的条件怎么提？能详细说说不？
现在正在看WENO的限制器。Hermit WENO的还没怎么看。到现在也没有看明白WENO为什么要搞个加权。在舒其望的论文里，加权后确实可以使得非“问题”单元恢复高精度，可是，很麻烦。而另外一个作者的论文，我看不出加权后有什么好处。还不如直接选一个斜率最小的拉倒。