当前位置:首页 >> 计算机软件及应用 >>

OpenFOAM 中的 div 与 snGrad 操作符










Home Archives Rss About Hom e Arc hives

Rss

About

?

Giskard's CFD Learning Tricks
CFD and Scientific Computing

2015-05-17 ?

O P E N F O AM

OpenFOAM 中的 div 与 snGrad 操作符
OpenFOAM 的方便之处之一是利用 C + +的类模板和函数重载等技术定义了很多各种离散 操作符,如 d i v, l a p l a c i a n, g r a d 等等。利用这些操作符,很容易就能对偏微分方程 进行离散,并构建起线性方程组。但是,这些操作符真正执行的运算,却需要结合有限体 积方法的本质来理解一番才能真正掌握。下面尝试着对 OpenFOAM 中的 d i v和 s n G r a d 操作符进行一点解读。

1. div 操作符的本质
d i v操作符表面看,是计算散度的,实际上,在OpenFAOM中,div 操作符的作用是加 和 ,比如说 ? ? (U U ) ,在OpenFOAM中表示为 f v m : : d i v ( p h i , U ),这段代码真正执 行的是∑f
U f ?f

运算,即将每个网格包含的面上的流率与速度乘积,然后再加起来。再

比如, t w o P h a s e E u l e r F o a m的 U a E q n方程有一项 是f v m : : s p ( f v c : : d i v ( p h i a ) , U a ),其对应的公式是 Ua (? ? Ua ) 。为什么是这样 呢?以下试图对背后的原理做一点解释: 单相流的动量守恒方程的微分形式如下:
?U ?t
T

+ ? ? (U U ) + ? ? dev(?νeff (?U + (?U )

)) = ??p + Q









为了使用有限体积方法,需要将动量方程写成积分形式,即
?U ?t
T


V

[

+ ? ? (U U ) + ? ? dev(?νeff (?U + (?U )

))] dV = ∫
V

[??p + Q] dV

这里只分析∫V

? ? (U U )dV

这一项,利用高斯定理,可以将这一个体积分,转化成对包

围该体积微元的表面的面积分:∮?V (U U ) ? dS ,其中 dS 是面积微元矢量。又因为实际操 作中,一个体积微元总是由有限的几个面组成的,所以


?V

(U U ) ? dS = ∑ ∫
f Sf

(U U ) ? dSf = ∑ (U U )f ? Sf
f

,其中
1 (U U )f = mag(Sf )


Sf

(U U ) ? dSf

。 这里做一个近似:
(U U )f ≈ (U f U f )

于是得到
∑ (U U )f ? Sf ≈ ∑ (U f U f ) ? Sf
f f

运用张量计算规则[uv ? w]

= u(v ? w)

,得到

∑ (U f U f ) ? Sf = ∑ U f (U f ? Sf ) = ∑ U f ?f
f f f

综上,OpenFOAM中的代码 f v m : : d i v ( U U )对应的公式其实是∫V 据推导,在有限体积方法中 ∫V 格的所有面上的面心上得到
? ? (U U )dV = ∑ (U f ?f )
f

? ? (U U )dV

,而根

,所以, f v m : : d i v ( U U )
U

实质上进行的运算是,把网格当作体积微元,将网格中心的速度
Uf

插值到包围该网

,并计算每个面上的速度通量

?f

,然后返回每一个

面上的速度与通量乘积的加和 。Uf 的计算需要用到本网格与邻近网格的速度值,离散格 式的作用就体现在如果利用本网格与邻近网格的速度得到面上的速度。









再来看上面提到的另一项Ua (? ? Ua ) ,根据上面类似的推导

V

[U a (? ? U a )] dV = U a ∫
V

(? ? U a )dV

注意,这里之所以能这样变换,是因为在一个体积微元dV内,Ua是常数。根据高斯定理

V

(? ? U a )dV = ∮
?V

U a ? dS = ∑ (U a )f ? Sf = ∑ (?a )f
f f

于是得到


V

[U a (? ? U a )] dV = U a ∑ (?a )f
f

这一项在OpenFOAM中的表达是 f v m : : s p ( f v c : : d i v ( p h i a ) , U a ),含义就很明显了, 这一项相当于是一个系数与需要求解的量Ua 的乘积,所以被当作隐式的源项来处里。注 意我先前以为把(? ? Ua ) 当作显式处理是人为简化的结果,其实不然,这是自然而然的结 果。而在这里也可以看出,∑f (?a )f 对应的代码是 f v c : : d i v ( p h i a ),也印证了上面 的观点,即 d i v操作符本质上是在作加和运算 。

源项

2. grad 与 snGrad
在t w o P h a s e E u l e r F o a m中,
?α α

在两个不同的地方用了两种不同的表示。
?α α

? ? {[νeff , a

][U a ]}

对应的代码是: f v m : : d i v ( p h i R a , U a ),其中

1 2

p h i R a = f v c : : i n t e r p o l a t e ( n u E f f a ) * m e s h . m a g S f ( ) * f v c : : s n G r a d ( a l p h a ) / f v c : : i n t e r p o l a t e ( a l p h a+s c a l a r ( 0 . 0 0 1 ) ) ;

而另一项











?α [ α ] ? Rca

对应的代码是 f v c : : g r a d ( a l p h a ) / f v c : : a v e r a g e ( a l p h a+s c a l a r ( 0 . 0 0 1 ) )& R c a 下面是我对这个的理解 : 对于
?α α

? ? {[νeff , a

][U a ]}

其处理方法跟? ? (Ua Ua ) 是一样的,因为
[ νeff , a
?α α

?α α

也是一个矢量。 p h i R a相当于是

]

这个矢量的界面通量,类比于 p h i a。 p h i a的定义

是l i n e a r I n t e r p o l a t e ( U a )&m e s h . S f ( ),而 p h i R a理论上应该也可以定义成类 似于 l i n e a r I n t e r p o l a t e ( g r a d a l p h a )&m e s h . S f ( ),前提是要先定义一个 v o l V e c t o r F i e l dg r a d a l p h a = n u E f f a * f v c : : g r a d ( a l p h a ) / a l p h a。但是考虑 到界面通量本质上是先将一个 v o l V e c t o r F i e l d插值到面上,并乘以面积矢量,对于
?α α ?α α

,完全可以直接求出每个面上的

?α α

,然后乘以面积矢量,而不需要先建立体中心的

再插值到面上。 s n G r a d就是这样一个用来求面上的梯度量的函数,它求解面上的梯

度时,采用如下公式:
?N ? ?P |d|

(??)f =

即用相邻网格的值减去面所属网格的值,再除以两个网格中心矢量的模。注意这个处理对 于正交网格是精确的,对于非正交网格,只是一个近似处理。而且要注意,这样求得的面 上的梯度值已经是一个标量了。再回到 p h i R a的定义,既然 f v c : : s n G r a d ( a l p h a )已 经是标量了,那只要再乘以面积矢量的模 m e s h . m a g S f ( ),便是界面上的通量 了。 f v c : : i n t e r p o l a t e ( n u E f f a )和 f v c : : i n t e r p o l a t e ( a l p h a+ s c a l a r ( 0 . 0 0 1 ) )则分别是将volField插值到面。 再来看
?α [ α ] ? Rca









这一项是被当作显式的源项来处理,所以,我们需要得到的是一个volVcetorField。所以这 里将?α 处理成 volVectorField,再与作为 volTensorField 的 R c a 进行点乘,得到的就是 volVcetorField。因此,这里的
?α α

对应的代码是

f v c : : g r a d ( a l p h a ) / f v c : : a v e r a g e ( a l p h a+s c a l a r ( 0 . 0 0 1 ) )。注意这里的 f v c : : a v e r a g e ( ) 函数的定义如下:

1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 2 0 2 1 2 2 2 3 2 4 2 5 2 6 2 7 2 8 2 9

4 3t e m p l a t e < c l a s sT y p e > 4 4t m p < G e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h >> 4 5a v e r a g e 4 6( 4 7 4 8) 4 9{ 5 0 5 1 5 2 5 3 5 4 5 5 5 6 5 7 5 8 5 9 6 0 6 1 6 2 6 3 6 4 6 5 6 6 6 7 6 8 6 9 7 0 7 1 a v . i n t e r n a l F i e l d ( )= ) ; ) ) , m e s h , s s f . d i m e n s i o n s ( ) t m p < G e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h >>t a v e r a g e ( n e wG e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h > ( I O o b j e c t ( " a v e r a g e ( " + s s f . n a m e ( ) + ' ) ' , s s f . i n s t a n c e ( ) , m e s h , I O o b j e c t : : N O \ _ R E A D , I O o b j e c t : : N O \ _ W R I T E c o n s tf v M e s h &m e s h=s s f . m e s h ( ) ; c o n s tG e o m e t r i c F i e l d < T y p e ,f v s P a t c h F i e l d ,s u r f a c e M e s h > &s s f

G e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h > &a v=t a v e r a g e ( )









3 0 3 1 3 2 3 3 3 4 3 5 3 6 3 7 3 8 3 9 4 0 4 1 4 2 4 3 4 4 4 5 4 6 4 7 4 8 4 9 5 0 5 1 5 2 5 3 5 4 5 5 5 6 5 7 5 8 5 9 6 0 6 1 6 2 6 3 6 4 6 5

7 2 7 3 7 4 7 5 7 6 7 7 7 8 7 9 8 0 8 1 8 2 8 3 8 4 8 5 8 6 8 7} 8 8 8 9

( s u r f a c e S u m ( m e s h . m a g S f ( ) * s s f ) / s u r f a c e S u m ( m e s h . m a g S f ( ) ) ) ( ) . i n t e r n a l F i e l d ( ) ; t y p e n a m eG e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h > : : G e o m e t r i c B o u n d a r y F i e l d &b a v=a v . b o u n d a r y F i e l d ( ) ; f o r A l l ( b a v ,p a t c h i ) { b a v [ p a t c h i ]=s s f . b o u n d a r y F i e l d ( ) [ p a t c h i ] ; } a v . c o r r e c t B o u n d a r y C o n d i t i o n s ( ) ; r e t u r nt a v e r a g e ;

9 0t e m p l a t e < c l a s sT y p e > 9 1t m p < G e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h >> 9 2a v e r a g e 9 3( 9 4 9 5) 9 6{ 9 7 9 8 9 9 1 0 0 1 0 1 1 0 2 1 0 3} 1 0 4 1 0 5 1 0 6t e m p l a t e < c l a s sT y p e > 1 0 7t m p < G e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h >> ) ; t s s f . c l e a r ( ) ; r e t u r nt a v e r a g e ; t m p < G e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h >>t a v e r a g e ( f v c : : a v e r a g e ( t s s f ( ) )

c o n s tt m p < G e o m e t r i c F i e l d < T y p e ,f v s P a t c h F i e l d ,s u r f a c e M e s h >>









6 6 6 7 6 8 6 9 7 0 7 1 7 2 7 3 7 4 7 5 7 6 7 7 7 8 7 9 8 0 8 1 8 2 8 3 8 4 8 5 8 6 8 7 8 8

1 0 8a v e r a g e 1 0 9( 1 1 0 1 1 1) 1 1 2{ 1 1 3 1 1 4} 1 1 5 1 1 6 1 1 7t e m p l a t e < c l a s sT y p e > 1 1 8t m p < G e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h >> 1 1 9a v e r a g e 1 2 0( 1 2 1 1 2 2) 1 2 3{ 1 2 4 1 2 5 1 2 6 1 2 7 1 2 8 1 2 9 1 3 0} ) ; t v t f . c l e a r ( ) ; r e t u r nt a v e r a g e ; t m p < G e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h >>t a v e r a g e ( f v c : : a v e r a g e ( t v t f ( ) ) r e t u r nf v c : : a v e r a g e ( l i n e a r I n t e r p o l a t e ( v t f ) ) ; c o n s tG e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h > &v t f

c o n s tt m p < G e o m e t r i c F i e l d < T y p e ,f v P a t c h F i e l d ,v o l M e s h >> &t v t

可以看出, f v c : : a v e r a g e ( )函数的返回类型是 t m p < G e o m e t r i c F i e l d < T y p e , f v P a t c h F i e l d ,v o l M e s h > >,对应 a l p h a,返回类型就 是t m p < G e o m e t r i c F i e l d < s c a l a r ,f v P a t c h F i e l d ,v o l M e s h > >, 即t m p < v o l S c a l a r F i e l d >。 a v e r a g e函数的核心定义见代码71-74行,可见 a v e r a g e采用的是面积加权平均,即
?? ?



f

?f ? |Sf | |Sf |

? = ∑
f

。如果给 a v e r a g e的参数类型是surfaceFiled,那么直接计算平均值后返回volField,如果 给的参数的volField,那就先将volField插值到面,计算平均值后再返回volField,见代码









106-114行。

参考资料
1. http://www.openfoam.org/docs/cpp/ 2. OpenFOAM: A little User-Manual, Gerhard Holzinger, https://github.com/ParticulateFlow/OSCCAR-doc

#Code Explained #OpenFOAM

? Share

N E WE R
twoPhaseEulerFoam 全解读之一

OLDER
利用functionObjects 对指定区域内进行后处理

0条评论 还没有评论,沙发等你来抢 社交帐号登录: 微信 微博 QQ 人人 更多?

最新 最早 最热

说点什么吧…

发布
多说

C AT E G O R I E S
C++ (1)


相关文章:
更多相关标签: