星期六, 十二月 23, 2006

有限元法、有限差分法和有限体积法的区别

昨天在网上看到一篇关于有限差分,有限元和有限体积法区别的文章。现在贴在下面,也许对有兴趣的朋友有帮助(我对文章作了一些编辑)。

1. 有限差分法(Finite difference method)

有限差分方法(FDM)是计算机数值模拟最早采用的方法,至今仍被广泛运用。该方法将求解域划分为差分网格,用有限个网格节点代替连续的求解域。有限差分法以泰勒(Taylor)级数展开等方法,把控制方程中的导数用网格节点上的函数值的差商代替进行离散,从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法,数学概念直观,表达简单,是发展较早且比较成熟的数值方法。对于有限差分格式,从格式的精度来划分,有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑,可分为中心格式和逆风格式。考虑时间因子的影响,差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式,主要是上述几种形式的组合,不同的组合构成不同的差分格式。差分方法主要适用于结构网格,网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。

构造差分的方法有多种形式,目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式:一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等, 其中前两种格式为一阶计算精度,后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合,可以组合成不同的差分计算格式。

2. 有限元法(Finite element method)

有限元方法的基础是变分原理和加权余量法,其基本求解思想是把计算域划分为有限个互不重叠的单元,在每个单元内,选择一些合适的节点作为求解函数的插值点,将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式 ,借助于变分原理或加权余量法,将微分方程离散求解。采用不同的权函数和插值函数形式,便构成不同的有限元方法。有限元方法最早应用于结构力学,后来随着计算机的发展慢慢用于流体力学的数值模拟。

在有限元方法中,把计算域离散剖分为有限个互不重叠且相互连接的单元,在每个单元内选择基函数,用单元基函数的线形组合来逼近单元中的真解,整个计算域上总体的基函数可以看为由每个单元基函数组成的,则整个计算域内的解可以看作是由所有单元上的近似解构成。在河道数值模拟中,常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同,有限元方法也分为多种计算格式。从权函数的选择来说,有配置法、矩量法、最小二乘法和伽辽金法,从计算单元网格的形状来划分,有三角形网格、四边形网格和多边形网格,从插值函数的精度来划分,又分为线性插值函数和高次插值函数等。不同的组合 同样构成不同的有限元计算格式。对于权函数,伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数 ;最小二乘法是令权函数等于余量本身,而内积的极小值则为对代求系数的平方误差最小;在配置法中,先在计算域 内选取N个配置点 。令近似解在选定的N个配置点上严格满足微分方程,即在配置点上令方程余量为0。插值函数一般由不同次幂的多项式组成,但也有采用三角函数或指数函数组成的乘积表示,但最常用的多项式插值函数。有限元插值函数分为两大类,一类只要求插值多项式本身在插值点取已知值,称为拉格朗日(Lagrange)多项式插值;另一种不仅要求插值多项式本身,还要求它的导数值在插值点取已知值,称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标,有对称和不对称等。常采用的无因次坐标是一种局部坐标系,它的定义取决于单元的几何形状,一维看作长度比,二维看作面积比,三维看作体积比。在二维有限元中,三角形单元应用的最早,近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元,常采用的插值函数为有Lagrange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。

对于有限元方法,其基本思路和解题步骤可归纳为

(1)建立积分方程,根据变分原理或方程余量与权函数正交化原理,建立与微分方程初边值问题等价的积分表达式,这是有限元法的出发点。
(2)区域单元剖分,根据求解区域的形状及实际问题的物理特点,将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作,这部分工作量比较大,除了给计算单元和节点进行编号和确定相互之间的关系之外,还要表示节点的位置坐标,同时还需要列出自然边界和本质边界的节点序号和相应的边界值。
(3)确定单元基函数,根据单元中节点数目及对近似解精度的要求,选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的,由于各单元 具有规则的几何形状,在选取基函数时可遵循一定的法则。
(4)单元分析:将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近;再将近似函数代入积分方程,并对单元区域进行积分,可获得含有待定系数(即单元中各节点 的参数值)的代数方程组,称为单元有限元方程。
(5)总体合成:在得出单元有限元方程之后,将区域中所有单元有限元方程按一定法则进 行累加,形成总体有限元方程。
(6)边界条件的处理:一般边界条件有三种形式,分为本质边界条件(狄里克雷边界条件 )、自然边界条件(黎曼边界条件)、混合边界条件(柯西边界条件)。对于自然边界条件,一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件,需按一定法 则对总体有限元方程进行修正满足。
(7)解有限元方程:根据边界条件修正的总体有限元方程组,是含所有待定未知量的封闭方程组,采用适当的数值计算方法求解,可求得各节点的函数值。

3. 有限体积法(Finite volume method)

有限体积法(Finite Volume Method)又称为控制体积法。其基本思路是:将计算区域划分为一系列不重复的控制体积,并使每个网格点周围有一个控制体积;将待解的微分方程对每一个控制体积积分,便得出一组离散方程。其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分,必须假定值在网格点之间的变化规律,即假设值的分段的分布的分布剖面。从积分区域的选取方法看来,有限体积法属于加权剩余法中的子区域法;从未知解的近似方法看来,有限体积法属于采用局部近似的离散方法。简言之,子区域法属于有限体积发的基本方法。

有限体积法的基本思路易于理解,并能得出直接的物理解释。离散方程的物理意义,就是因变量在有限大小的控制体积中的守恒原理,如同微分方程表示因变量在无限小的控制体积中的守恒原理一样。 限体积法得出的离散方程,要求因变量的积分守恒对任意一组控制体积都得到满足,对整个计算区域,自然也得到满足。这是有限体积法吸引人的优点。有一些离散方法,例如有限差分法,仅当网格极其细密时,离散方程才满足积分守恒;而有限体积法即使在粗网格情况下,也显示出准确的积分守恒。就离散方法而言,有限体积法可视作有限单元法和有限差分法的中间物。有限单元法必须假定值在网格点之间的变化规律(既插值函数),并将其作为近似解。有限差分法只考虑网格点上的数值而不考虑值在网格点之间如何变化。有限体积法只寻求的结点值,这与有限差分法相类似;但有限体积法在寻求控制体积的积分时,必须假定值在网格点之间的分布,这又与有限单元法相类似。在有限体积法中,插值函数只用于计算控制体积的积分,得出离散方程之后,便可忘掉插值函数;如果需要的话,可以对微分方程中不同的项采取不同的插值函数。

星期四, 十二月 07, 2006

ANSYS 11.0终于即将发布

最近一两个月实在是太忙了,发文章很不正常,真是太内疚了。:(

这几天正在参加一个conference,晚上回来查看网页,无意中看到台湾虎门科技的网页(虎门科技是ANSYS在台湾的代理商)。居然得到了ANSYS即将发布的消息。台湾的发布会会在12月27号进行,也就是说ANSYS 11.0会在2-3周内发布。

星期一, 十二月 04, 2006

Fluent 6.3已经发布

昨天早上打开ANSYS的主页,发现上面的新闻稿说Fluent 6.3已经发布了。下面一个问题是ANSYS 11.0什么时候发布呢?现在ANSYS 11.0和10.0之间至少已经间隔16个月了。

Fluent6.3发布的相关新闻链接在这里。这里是关于Fluent 6.3新的特性的概述。

1. 网格、数值及并行计算:新的基于压力的耦合求解器,新的建模检查工具,新增加到多面体网格的支持,增加对64位windows的支持以及其它并行计算能力的提高。

F1赛车多面体网格

2. 传热、相变和辐射:面对面辐射模型应用范围和效率的提高,用户可以耦合离散坐标辐射模型和能量方程,支持多组分实际气体。

玻璃熔炉模拟

3.动网格:动网格可以应该在稳态问题中,动网格的定义更加直接和方便,滑移网格支持一边多个网格。

4.反应流:新的低速化学反应和微混合模型,更多的化学物质和反应可以使用无预混合或部分预混合模型模拟,更多模拟SOx、NOx生成的模型。

5.湍流和声学:更灵活的DES模型,LES模拟更多的时间统计结果,提高了RSM模型,用户可以自定义壁面函数。

6.多相流:瞬态多相流模拟精度得到了提高,欧拉、混合和VOF模型和壳传热模型兼容,扩展了欧拉模型的应用范围,在使用VOF模型时,可以使用新的界面跟踪模型,可用用UDF设定自由液面接触角,对离散相模型可以模拟多组分例子的气化,在使用动网格时,粒子可以从动网格表面射入。

7.后处理:一些新的报告功能,增强了和第三方后处理软件的协作,Fluent并行计算结果可以输出到tecplot。

8.自定义工具:扩展了用户定义标量的应用范围,更多的UDF链接功能,对扩展模块(add-on modules)的改进。

星期日, 十一月 12, 2006

ANSYS计算节点内力的宏

今天根据ansys.net上的一个宏稍微修改了一下,做了一个计算节点内力的宏,贴在下面。这个宏的原始作者是Shen-Yeh Chen。这里是原始描述和原始宏。

这个宏可以计算指定节点上的内力。通常情况下,ANSYS只计算被限定自由度的节点的反应力。如果模型里有两个材料,并且这两个材料用连续网格模拟(而不是接触),有些时候我们想知道这两个材料间的相互作用力,而在ANSYS里面并没有直接的命令。使用这个宏可以做到这一点。

这个宏的使用法方法是:0.当然是把这个宏命名并放在相应的目录下;1.选择要计算内力的节点;2.选择这些节点一边的单元,因为这些节点代表一个截面,而这个截面会有两边有单元。选择不同的边,内力的方向是不同,如果两边都选的话,内力就会平衡为零了。因此必须要先选择一边,但是不必要只选和节点相邻的一层,因为宏里面已经有命令去选择和节点相连的单元;3.运行宏;4.在工作目录下,会有一个x_nfor2f.ans文件,这个文件的最后一行,就是在x,y,z三个方向的内力分量了。

这个宏必须要在通用后处理器中运行。

! ========================================

! Shen-Yeh Chen

! 05/01/1998

! 12/10/2001

! =========================================

!

! This macro will transfer the element force on the selected

! node into cross-section force on each node, and write

! all the force on each node into the "F" command. Note that

! the macro go through each SELECTED node, sum up all SELECTED

! elements that are attached to the SELECTED node. Then the

! force on each nodes is added a negative sign and printed

! into the 'F' command.

! ---------------------------------------------------------

! Syntax :

! x_nfor2f

!

! Input :

! None

!

! Output :

! file : x_nfor2f.ans : nodal force for the section load

! x_nfor2f.cdb : nodal force & nodal coordinates/ID in

! CDB format

! nodal components : x_nfor2f : defined in x_nfor2f.cdb

! X_NFOR2F_FX : summation of FX

! X_NFOR2F_FY : summation of FY

! X_NFOR2F_FZ : summation of FZ

! Temporary components

! X_NFOR2F

! X_NFOR2F_E

!

! Temporary components

! x_nfor2f.db

!

! Temporary variables

! AR30,AR31,AR32,AR33

!

! External Programs called

! None

!

!

! ========================================

/POST1

*GET,ZNNODE,NODE,,NUM,MAX

CM,X_NFOR2F,NODE

CM,X_NFOR2F_E,ELEM

ESLN,R

!

AR30=0

*CFOPEN,x_nfor2f,ans

*VWRITE

('/PREP7')

X_FOR2F_FX=0.0

X_FOR2F_FY=0.0

X_FOR2F_FZ=0.0

!

! LOOP THROUGH ALL SELECTED NODES

:LINE01

AR30=NDNEXT(AR30)

*IF,AR30,EQ,0,THEN

*GO,:LINE02

*ENDIF

NSEL,R,NODE,,AR30

FSUM

!

!

! GET THE FORCE

*GET,AR31,FSUM,0,ITEM,FX

*GET,AR32,FSUM,0,ITEM,FY

*GET,AR33,FSUM,0,ITEM,FZ

!

AR31=-AR31

AR32=-AR32

AR33=-AR33

X_FOR2F_FX=X_FOR2F_FX+AR31

X_FOR2F_FY=X_FOR2F_FY+AR32

X_FOR2F_FZ=X_FOR2F_FZ+AR33

!

! WRITE TO FILE

*VWRITE,AR30,AR31

('F,',F6.0,',FX,',E12.5)

*VWRITE,AR30,AR32

('F,',F6.0,',FY,',E12.5)

*VWRITE,AR30,AR33

('F,',F6.0,',FZ,',E12.5)

!

! SELECT THE NODES

CMSEL,,X_NFOR2F

*GO,:LINE01

!

!

! SAVE THE STATE TO x_nfor2f.db

:LINE02

*VWRITE,X_FOR2F_FX,X_FOR2F_FY,X_FOR2F_FZ

('FX = ',E12.5,', FY = ',E12.5,', FZ = ',E12.5)

*CFCLOS


星期一, 十一月 06, 2006

STAR-CD和STAR-CCM+ :冠军的选择

最近CD-Adapco网站上登出了一则新闻。新闻的题目是:STAR-CD and STAR-CCM+: Still the choice of champions (STAR-CD和STAR-CCM+:冠军的选择)。这则新闻被很多网站转载。新闻说的是2006赛季F1冠军车队雷诺车队使用的是CD-Adapco提供的CFD软件。我每天上班的路上看的一份小报几乎每天都要介绍F1大赛,我对当中一些人得的名字还是比较熟悉的。

新闻稿中说,自2001年雷诺车队成立以来,CD-Adapco就是雷诺车队的CFD软件提供商。RS26引擎 和R26底盘 的设计都和CD-Adapco的STAR-CD或STARE-CCM+软件有关。CD-Adapco是F1大赛中最成功的CFD软件提供商,因为,在2006赛季有得分的车队中,有70%的车队使用STAR-CD或者STAR-CCM+作气动或者引擎流动模拟。看来STAR-CD在汽车行业的优势是当之无愧的。关于F1大赛的更多信息可以参见搜狐体育的F1频道 。

星期日, 十一月 05, 2006

多面体网格的优点

多面体网格(polyhedral mesh)是最近CFD求解器开发的热点,它已经被使用在CD-Adapco的STAR-CCM+中,STAR-CD V4也支持多面体网格,Fluent V6.3也将支持多面体网格。其实在CFX中,由于采用基于节点的控制体积(Vertex-based control volume),用户提供给求解器的网站是正常的四面体或六面体网格,而在求解器内部会生成多面体网格,并且声称已经使用这个技术超过15年。下面的文章翻译自CD-Adapco上面的一个新闻发布(链接),作者是Milovan Peric和Stephen Ferguson。

自CD-Adapco推出STAR-CCM+求解器和自动多面体网格生成器以来, 我们经常被问到一个问题:为什么我们会采用多面体网格,而不是四面体网格? 在这篇文章里,我们会就这个问题提供答案和实例。

四面体单元是最简单的体积单元;他们的表面是平面碎片 (plane segment), 因此面和体的中心点很容易定义。四面体网格也是相对容易自动生成的。它们几乎是自动网格生成技术中的标准应用,主流CFD求解器都支持四面体网格。但是,在另一方面,四面体网格不能被拉伸太多,因此为了保证在边界层、长的管道、或小的缝隙的精确度, 需要使用比结构化网格(六面体)多的多的四面体网格。棱形网格层可以部分解决这个问题。

四面体控制体只有四个相邻网格,在使用标准近似(线性形函数)来计算单元中心的梯度时,会比较困难。一个困难是相邻网格节点的空间位置,因为它们可能会几乎处在一个平面上,使得几乎不可能计算垂直于该平面的梯度。当单元邻近边界时会有另外一个可能碰到的问题:即使只有一个面是边界面,其它三个面的分布可能会很不好;尤其是在邻近线或者角落时,一个单元可能只有两个或者一个相邻单元,这不光会降低精确度,而且会引起严重的数值问题。

为了在四面体网格上实现精确求解和好的收敛性,需要特别的离散技术和大量的单元。这些解决方案并不是最优化的;第一种方法使得程序更加复杂、更难扩展和维护,而第二个方法增加了内存和计算时间需求。多面体网格提供了和四面体网格一样的自动网格生成功能,并且能解决这些缺点。

多面体网格的主要优点是它有多个相邻单元(通常10数量级) ,和四面体网格相比,使用多面体网格能更好的近似梯度(使用线性形函数以及最近单元的信息)。即使靠近线和角,一个多面体单元也可能有几个相邻单元,从而可以合理的预测梯度和当地流动分布。更多的相邻单元意味着更多的内存使用量和更多的计算操作,这可以从更高的精度中得到补偿。下面将要讨论这个问题。

多面体网格也不如四面体网格多拉伸敏感。智能化的网格生成和优化技术也提供了无限的可能性:通过加入新的节点、单元或面,单元可以自动结合、分裂、或修改。实际上,网格质量在将来会有显著提高,并提高求解器效率和精确度。而且,求解器中在使用四面体网格时需要的特别处理将不再必要。在使用多面体网格时,要模拟当地网格细化、滑移网格界面、周期性边界条件等,只需要生成特别的多面体网格,对求解器来说它们都是一样的。

多面体网格在处理回流时更有优势。比如测试显示在模拟正方体空腔驱动流时,要得到指定的精度所需的多面体单元数量要比笛卡尔六面体(通常认为对长方体计算域,六面体是最理想的)少得多。这可以解释为:对六面体单元,有三个最佳流动方向可以达到最大精度(和三对平行面垂直的方向) ;对有十二个面的多面体单元,有六个最佳流动方向,并且有更多的相邻单元,因此可以用较少的单元数得到较高的精确度。对于不同网格类型特点和一个测试结果可以参考Peric (2004) 发表的文章。

很多实际应用的比较分析显示,要达到和四面体网格同样的精确度,使用多面体网格只需要四分之一的单元数、一半的内存、十分之一到五分之一的计算时间。另外,多面体网格有更好的收敛性,通常不需要调整默认的求解器参数,如图1至3所示。

图1. 模拟发动机水套内流动所用的多面体网格


图2. 最细的四面体网格(左图)和最细的多面体网格(右图)所预测的水套壁面压力分布

图3. 不同四面体网格和多面体网格所预测的水套入口出口压力差

这个测试模拟的是发动机水套;监测的是进口和出口的压力差。计算分别在6个多面体网格(单元数21872至593888)和6个四面体网格(单元数39587至2322106)上进行。在所有情形中近壁面处均使用了棱形网格层。图1是一个多面体网格。图2是在最细的多面体网格和最细的四面体网格上计算所得的压力分布。因为网格很小,这两个计算结果非常相似。为了分析结果对网格的依赖性,图3显示了在不同网格上计算得到的压差。该图显示这两种网格类型都可以收敛得到网格无关解。所有的情形使用了同样的离散和求解方法(对流项采用二阶迎风差分格式)。

重要的是,在多面体网格上得到的解总是比有相近网格数量的四面体网格要精确。实际上,从有65513单元的多面体网格上得到的解比从有393273个单元(大约6倍多)的四面体网格上得到的解还要稍微精确。多面体网格所需的计算时间只是得到同等精度的四面体网格的十分之一。到目前为止所进行的比较得到了相似的结果,这验证了CD-Adapco发展这一技术的决定是正确的。随着CD-Adapco的全自动多面体网格生成器的面世,多面体网格将会成为业界标准。

参考文献:

1 M. Peric: Flow simulation using control volumes of arbitrary polyhedral shape, ERCOFTAC Bulletin, No. 62, September 2004.

一周没有发言了

原因是参加了一个关于CAE的conference。