自定义初始场

前言

以一维爆轰波为例，大家可以点击下载测试算例1查看。计算域左端设置一小段高温高压区域用于初始起爆，其他区域设置新鲜气体，在区域 $0.1m \leq x \leq 0.2m $ 设置温度场如下：

$$
T = 400 + 100 \sin \left[2\pi\frac{(x-0.1)}{0.1}\right]
\tag{1.1}
$$

这里仅设置了温度的初始场为自定义函数，相信大家可以举一反三。

查看做的修改

  1. 首先打开0.orig/T文件，第一处修改是format设为ascii，如果是binary记得修改：

    1
    2
    3
    4
    5
    6
    7
    8
    FoamFile
    {
        version     2.0;
        format      ascii; // not "binary"
        class       volScalarField;
        location    "0";
        object      T;
    }

    第二处修改是本次的重点 codeStream 编程，详见算例文件。
    下面仅贴出主体部分代码，大家可根据需要修改：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    code
    #{
        const IOdictionary& d = static_cast<const IOdictionary&>(dict);
        const fvMesh& mesh = refCast<const fvMesh>(d.db());
        scalarField value(mesh.nCells(), 400.); // 400K by default
        scalar Pi = acos(-1.);
        forAll(value, i)
        {
            // Access cell centers coordinates (x,y,z)
            // only use x in this 1D simulation
            const scalar x = mesh.C()[i][0]; 
            // const scalar y = mesh.C()[i][1];
            // const scalar z = mesh.C()[i][2];    

            // you can define value as a function of (x,y,z)        
            if (x >= 0.1 && x <= 0.2)
            {
                value[i] = 400. + 100*sin(Pi*(x-0.1)*20);                    
            }
        }
        os.writeEntry("", value);
    #};

  2. 然后打开system/setFieldsDict文件，其实里面就注释掉了一行，避免覆盖之前做的修改： 

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    defaultFieldValues
    (
        // volScalarFieldValue T $init.default.T
        volScalarFieldValue p $init.default.p
        volScalarFieldValue H2 $init.default.H2
        volScalarFieldValue O2 $init.default.O2
        volScalarFieldValue H2O $init.default.H2O
        volScalarFieldValue N2 $init.default.N2
        volScalarFieldValue Ydefault 0
    );

算例运行

在算例目录下直接运行./Allrun即可。

补充说明

  1. 上面T为标量场，如果你想对矢量场（比如U）进行修改，可以参考下面的写法：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    code
    #{
        const IOdictionary& d = static_cast<const IOdictionary&>(dict);
        const fvMesh& mesh = refCast<const fvMesh>(d.db());
        vectorField value(mesh.nCells(), vector(0, 0, 0)); 
        forAll(value, i)
        {
            const scalar x = mesh.C()[i][0]; 
            const scalar y = mesh.C()[i][1];
            const scalar z = mesh.C()[i][2];    
            
            // no sense, just for example   
            value[i] = vector(sin(x), sin(y), sin(z));      
        }
        os.writeEntry("", value);
    #};

  2. 本帖基于 OpenFOAM-v1812 版本，不同版本的 codeStream 在写法上可能略有出入。