找回密码
 注册
查看: 18173|回复: 34

OpenFOAM的程序开发初步

  [复制链接]
发表于 2010-1-14 08:10:29 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?注册

x
一.OpenFOAM应用的类型:
使用OpenFOAM进行CAE模拟的,大致可分为三种类型:
1)直接利用OpenFOAM的标准的求解器进行模拟,把OpenFOAM替代商业软件来使用,OpenFOAM已基本具有这样的功能和人气,与Fuent,Star-CD等相比较,OpenFOAM显然具有更高的求解效率和灵活性。
2)用户自定义求解器,即利用OpenFOAM的基本类库,如finiteVolume,OpenFOAM库来按照自己的求解流程来编写针对某类应用的求解器。用户需要开发的求解器就是类似于在OpenFOAM的applications中所看到的标准求解器icoFOAM,simpleFOAM等。显然这一需求是非常大的,从OpenFOAM问世以来,已有很多用户定义了自己的求解器。这类需求的特点是,并不需要特别关心,离散和求解的最底层的知识,如时间项离散,空间项离散等,关注的重点是求解的步骤或者流程。在编程中,通常是顶层的求解流程的开发,在多数情况下可以不编译OpenFOAM的finiteVolume和OpenFOAM库。这种顶层的求解器的开发,是我们以前常常忽略的,或者是以前没有能力做到的。需要指出的是,商业软件中的所谓udf,user subroutine和这是不可相比的。
3)用户自己定义离散方法等。对于研究离散格式、代数求解器等人来说,更关注时间项ddt,扩散项Laplacian,对流项div是如何离散的,能否有更高效更高精度的离散方法,这需要修改finiteVolume库和OpenFOAM库中对应的代码。尤其是对流项,尽管OpenFOAM已经提供了基于NVD和TVD的模板和40多种有名的高阶高精度格式,但可以预见,这仍然是不够的,毕竟对流项的离散仍然是目前CFD的重点研究方向。
可以肯定的是,目前有很多人关注类型2的应用,毕竟将OpenFOAM当成Fluent或Star-CCM来使用,并不见得方便。但是将OpenFOAM作为类库来构建自己的求解器,这是其它软件无法实现的

评分

1

查看全部评分

 楼主| 发表于 2010-1-14 08:11:02 | 显示全部楼层

OpenFOAM的程序开发初步

二.OpenFOAM程序开发的基本知识
2.1OpenFOAM的基本术语
重要的环境变量:
$WM_PROJECT_USER_DIR  ―― OpenFOAM的用户目录
$FOAM_TUTORIALS          ------OpenFOAM的算例目录
$ FOAM _SRC                ------OpenFOAM库的源程序目录
$ FOAM_APP                ------ OpenFOAM的求解器目录
$ FOAM_APPBIN            ------- OpenFOAM的求解器执行文件目录
$ FOAM_RUN                ------用户的算例目录
重要的shell:
run    =  cd to $FOAM_RUN
src    =  cd to $FOAM_SRC
app    =  cd to $FOAM_APP
util    =  cd to $FOAM_APP/utilities
sol    =  cd to $FOAM_APP/solvers
tut    =  cd to $FOAM_TUTORIALS

求解器的基本文件结构
appName            包含求解器源代码的目录
+appName.C        求解器主程序
  +CreateFields.H  场变量的声明和初始化
  +Make/          编译指令
+files        编译需要的源程序文件和生成的目标文件
+options      编译选项,如链接库等
appName/appName.C是求解器的主程序
appName/createFields.H声明变量,并从文件中读入初值,如p,物性。
appName/Make/files 所有源程序的名称,一个文件一行,最后一行是目标代码的名称和存放位置,EXE=$(FOAM_USER_APPBIN)/appName
appName/Make/options设定查找头文件和库的路径,EXE_INCS,和需要链接的库EXE_LIBS

算例的基本文件结构
case/                        算例目录
+0/                        包含初始和边界条件
+constant/                  包含初次读入后,不随时间变化的数据
+polyMesh/            包含多面体网格数据
+transportProperties/      包含物性数据
  +system/                  包含计算控制和离散格式设定
+controlDict            包含计算控制,如时间步长等
+fvSchemes              包含离散格式设定
+fvSolutions            包含代数求解器或SIMPLE,PISO算法设定
具体而言
case/0                      每个需求解的变量需要一个文件设定其初始边界条件
case/constant/polyMesh        网格数据,如owner neighbour points faces boundary
case/system/transportProperties  物性数据
case/system/controlDict        设定起始终止时间,时间步长,输出控制
case/system/fvSchemes        设定程序用到的每个微分算子的离散格式
case/system/fvSolution      为每个变量选择代数方程求解器/收敛精度及PISO等算法设定
回复 支持 1 反对 0

使用道具 举报

 楼主| 发表于 2010-1-14 08:12:44 | 显示全部楼层
三.OpenFOAM程序开发的理论知识
作求解开发,必须能写出需要求解的控制方程及其定解条件,并且对于如何求解方程或方程组的步骤已经明确。
这些流体力学、传热学以及相关的理论是必需的,所谓连续介质力学中的数学模型,控制方程和定解条件就是表示它的语言。
在这里是不可能说清楚的,这要看个人的功底了。

四 .OpenFOAM程序开发的最简单的例子
下面采用OpenFOAM来开发一个用户自己的求解器。主要是利用OpenFOAM的标准求解器icoFoam,用户不需要写任何代码,只为为了熟悉OpenFOAM程序开发的环境和步骤。
步骤:
1)    将icoFoam目录拷贝到新的目录
可采用下面的Linux的命令实现:
到OpenFOAM的incompressible目录
cd applications/incompressible
cp –r icoFoam myicoFoam
以上只是复制目录icoFoam到新的位置,并且新目录名为myicoFoam
cd myicoFoam
进入新的目录,查看一下,可以看到里面的文件和icoFoam中是否一样

2)    原文件改名,并且删除依赖文件
将icoFoam.C改名myicoFoam.C
mv icoFoam.C myicoFoam.C
删除依赖文件
rm icoFoam.dep
3)    修改编译文件files和options
进入Make目录,打开files文件
    将
icoFoam.C      源程序文件名
EXE = $(FOAM_APPBIN)/icoFoam  可执行文件名
修改为
myicoFoam.C      源程序文件名
EXE = $(FOAM_APPBIN)/myicoFoam  可执行文件名

此例中options不需修改,可以打开看看
EXE_INC = \    头文件包含
    -I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \  链接库
-lfiniteVolume
4)删除原来的obj文件
    rm –rf linuxGccDPOpt
cd ..
5)编译
  wmake
6) 检验一下
  到tutorial目录,检验一下
  myicoFoam . cavity

六.OpenFOAM程序开发――例子一:在icoFoam中加入温度场求解
准备:
能量控制方程:
  dT/dt+div(den*U*T)=div(a gradT)

在壁面上给定值条件。
需要解决的问题:
a)如何创建标量场,T
b) 如何创建物性,a
c)如何定义温度方程,并求解
d) 如何在算例中设定T和a
e)如何设定T的离散格式
f)如何设定T的求解器的收敛标准等
步骤:
1)创建程序需要的新物性和新变量场
打开myicoFoam.C可以看到,程序开始运行时调用CreateFields.H,创建变量场。
打开CreateFields.H,可以看到程序首先从transportProperties文件中读入物性,
    Info<< "Reading transportProperties\n" << endl;

    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",              从字典文件transportProperties读入
runTime.constant(), //transportProperties文件位于目录runTime.constant()中
            mesh,            网格对象                  
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
);  创建了Iodictionary类型对象 transportProperties

    dimensionedScalar nu        //首先读入粘性系数
    (
        transportProperties.lookup("nu")
    );  创建有量纲标量nu,nu通过从字典transportProperties查找”nu”来赋值

可以加上新方程需要的物性
    dimensionedScalar DT        //首先读入热扩散率
    (
        transportProperties.lookup("DT")
    ); 创建有量纲标量DT,DT通过从字典transportProperties查找”DT”来赋值
此外还要从createFields中读入p,U场,我们要加入的新的变量场为温度场T,最快的加入温度场的方法是拷贝p场的代码,修改为
    Info<< "Reading field T\n" << endl;
    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
这样,创建了新的vol标量场T,从文件T中读入。
对于T的创建具体解释如下:
a)创建了标量场T
b)T通过读(IOobject::MUST_READ)在runTime.timeName()目录下名称为“T”的文件创建,在开始计算时,runTime.timeName()是contorlDict中设定的startTime值决定的。
c)T将自动写入(IOobject::AUTO_WRITE)计算结果到runTime.timeName()目录中,runTime.timeName()随迭代是变化的,写入控制由contorlDict中设定。
d)T是定义在mesh对象上的,这意味着T在内部cell上有值internalField,在边界上还需要边界条件,这与polyMesh/boundary中要一致。

2)在求解器中加入新的求解方程
    下一步回到myicoFoam.C加入新的微分方程,由于温度场依赖于速度场,可放在PISO循环后面。
#          include "continuityErrs.H"

            U -= rUA*fvc::grad(p);
      
            U.correctBoundaryConditions();
      
//          Add the temperature equation
   
            fvScalarMatrix Teqn 温度是标量方程  
          (
              fvm::ddt(T)
            + fvm::div(phi, T)      要用到界面流量
            - fvm::laplacian(DT, T)  扩散项
          );
              TEqn.solve();          求解

3)编译
wmake
4)在算例中加入新方程的初始和边界条件
4.1拷贝一个cavity算例到mycavity
4.2修改transportProperties字典文件,设定DT
cd constant
修改transportProperties文件,前面已提到DT要从该字典文件读入。设定DT=0.002m2/s
DT              DT [0 2 -1 0 0 0 0] 0.002;
4.3修改T文件,设定初始值和边界
cd 0 进入0目录
拷贝一个T文件
cp p T
修改T文件为
FoamFile
{
    version        2.0;
    format          ascii;

    class          volScalarField;
    object          T;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


dimensions      [0 0 0 1 0 0 0];

internalField  uniform 300;    初始内部点为300℃

    movingWall      
    {
type            fixedValue;
value uniform  350.;        边界为350℃
    }

    fixedWalls      
    {
type            fixedValue;
value uniform  300.;        边界为300℃

5)修改离散格式和代数求解器求解控制文件
A进入system目录
由于温度方程有非稳态项,对流项,扩散项,分别要在ddt,div,laplacian中设置
打开fvSchemes文件,添加
divSchemes
{
    default        none;
    div(phi,U)      Gauss upwind;
    div(phi,T)      Gauss upwind;
}

laplacianSchemes
{
    default        none;
    laplacian(nu,U) Gauss linear corrected;
    laplacian(DT,T) Gauss linear corrected;
    laplacian((1|A(U)),p) Gauss linear corrected;
}
在fvSolution中设置代数求解器选项
      T PBiCG
    {
        preconditioner  DILU;
        tolerance        1e-06;
        relTol          0;
    };
注意T方程形成的矩阵是非对称的,不要用PCG和DIC
6)运行
myicoFoam . mycavity
 楼主| 发表于 2010-1-14 08:13:10 | 显示全部楼层
七.OpenFOAM程序开发――求解器的详细分析1
进入icoFoam目录
可以看到
createFields.H  icoFoam.C  icoFoam.dep Make/
Make/为wmake编译所需的文件
IcoFoam.C为主程序文件,它包含createFields.H
编辑icoFoam.C
可以看到icoFoam.C首先引入的头文件为fvCFD.H。
所以你可以看到,在编译选项options中
EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude  //fvCFD.H的存放目录

EXE_LIBS = \
    -lfiniteVolume                  //需要链接的库

找到fvCFD.H,编辑,可以看出这些是主程序必须的类库
#ifndef fvCFD_H
#define fvCFD_H

#include "parRun.H"

#include "Time.H"              时间类
#include "fvMesh.H"            网格类
#include "fvc.H"                fvc类
#include "fvMatrices.H"          fvMatrix类
#include "fvm.H"                fvm类
#include "linear.H"
#include "calculatedFvPatchFields.H"
#include "fixedValueFvPatchFields.H"
#include "adjustPhi.H"
#include "findRefCell.H"
#include "mathematicalConstants.H"

#include "OSspecific.H"
#include "argList.H"

#ifndef namespaceFoam
#define namespaceFoam
    using namespace Foam;
#endif

#endif

再看看icoFoam的程序体,了解一下求解程序的结构

#include "fvCFD.H"        ――――――――――――――――(头文件)
                        通常位于main函数前,是程序所需的类的定义
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
                          ―――――――――――――――(包含文件)
#  include "setRootCase.H"  

#  include "createTime.H"
#  include "createMesh.H"
―――――――――――包含文件通常是程序片断,如创建时间、创建网格等

―――――――――――――――(求解器代码)――――――
#  include "createFields.H"         
需要根据应用,单独写的代码,如"createFields.H"和Main,以及Ueqn,pEqn等
    ―――――――――――――――――――――――――――――――――――――
#  include "initContinuityErrs.H"
                        
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

。。。。。
}
八.OpenFOAM程序开发――求解器的详细分析2
a.场变量的定义
引用前面的温度场
    Info<< "Reading field T\n" << endl;
    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
例如
volScalarField CO2
(
IOobject
(
"CO2",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
// Optional declaration, this can be done by accessing a file in "case/0/",量纲可在文件中读
// dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), value)
);
b.场变量的使用
场变量有定义在内部cell上的值,也有边界上的值
例如给组分限值
Example of a mass fraction limiter used in this project:
// Initialize the variable Y_i for use in a loop
scalarField& CO2Internal = CO2.internalField();  引用内部点

// Loop for all mesh points  遍历内部点
forAll(CO2, celli)
{
// Limits the mass fraction to a positive number
if (CO2Internal[celli] < 0.0)
{
CO2Internal[celli] = 0.0;
}
// Limits the mass fraction to max 1.0
if (CO2Internal[celli] > 1.0)
{
CO2Internal[celli] = 1.0;
}
}

c.定义输运方程
OpenFOAM 定义方程时要选择一种类型的fvMatrix,有fvScalarMatrix和fvVectorMatrix
离散格式在case/system/fvSchemes.中设定
// Define a ScalarMatrix as a object
fvScalarMatrix CO2Eqn 定义系数矩阵
(
fvm::div(phi, CO2)            对流项离散fvm::div
- fvm::laplacian(turbulence->nuEff(),CO2)  扩散项离散fvm::div
== S_CO2                            源项
);

// Apply underrelaxation to the equation
// Under relaxation factors defined in file: fvSolution
CO2Eqn.relax();  松弛
CO2Eqn.solve();  求解
 楼主| 发表于 2010-1-14 08:21:58 | 显示全部楼层
[img]c:\[1.jpg]

solvers.rar

137.9 KB, 下载次数: 1265

 楼主| 发表于 2010-1-14 08:23:38 | 显示全部楼层
tut files

tuts.rar

83.96 KB, 下载次数: 1093

tut file

发表于 2010-1-15 08:29:01 | 显示全部楼层
好贴
支持下,谢谢了
发表于 2010-1-16 00:30:14 | 显示全部楼层
不错 顶一个!
发表于 2010-3-7 15:25:30 | 显示全部楼层

回复 1# liuhuafeifei 的帖子

写得好,以后多发这种好贴呀!
发表于 2010-3-7 20:04:02 | 显示全部楼层
强。
学习。
发表于 2010-3-18 08:45:02 | 显示全部楼层
up
发表于 2010-4-4 08:40:47 | 显示全部楼层
发表于 2010-4-11 11:54:20 | 显示全部楼层
太好了
发表于 2010-5-6 10:56:36 | 显示全部楼层
谢谢楼主
好帖  强烈支持
发表于 2010-5-20 10:43:47 | 显示全部楼层
谢谢楼住,学习中,
您需要登录后才可以回帖 登录 | 注册

本版积分规则

快速回复 返回顶部 返回列表