数值风云

OpenFOAM造波-U文件内边界问题

在使用codedFixedValue生成一段波浪时,alpha.water文件内已写出,U文件内对应的边界条件不知道该怎么给。
参考了wiki上的一个教程案例(这里只贴出了波浪生成边界的代码,其他边界应该没有影响)
alpha.water

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

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    left
    {
        type            codedFixedValue;
        value           uniform 0;
        name            wave_alpha;
        code            
        #{
            const fvPatch& boundaryPatch = patch(); 
            const vectorField& Cf = boundaryPatch.Cf(); 
            scalarField& field = *this; 

            field = patchInternalField(); 	//take value from initialization

            const scalar l = 5.;
            const scalar A = 0.1;
            const scalar gs = 9.81;
            const vector gv(0.,9.81,0);

            const scalar t = this->db().time().timeOutputValue();

            scalar k = 2.0*3.14159/l;
            scalar w = sqrt(k*mag(gv));

            forAll(Cf, faceI)
            {
                const scalar y = Cf[faceI].y();

                if (y <= A*cos(-w*t)+0.5*k*A*A*cos(2*(-w*t)) )
                {
                    field[faceI] = 1.;
                }
            }
        #};
    }
    ...
// ************************************************************************* //

其中公式 y <= A*cos(-w*t)+0.5*k*A*A*cos(2*(-w*t)) 应该是代表二阶斯托克斯波,我对这一部分做了修改,将公式改为两个已知的数组,波高时间序列,具体如下

修改后的alpha.water文件

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

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;
boundaryField
{
    left
    {
        type            codedFixedValue;
        value           uniform 0;
        name            wave_alpha;
        code            #{
            const fvPatch& boundaryPatch = patch(); 
            const vectorField& Cf = boundaryPatch.Cf(); 
            scalarField& field = *this; 

            field = patchInternalField(); 	//take value from initialization
            const scalar t = this->db().time().timeOutputValue();
            //时间数组
            float time[] = {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,4,4.1,4.2,4.3,4.4,4.5,4.6,4.7,4.8,4.9,5,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6,6.1,6.2,6.3,6.4,6.5,6.6,6.7,6.8,6.9,7,7.1,7.2,7.3,7.4,7.5,7.6,7.7,7.8,7.9,8,8.1,8.2,8.3,8.4,8.5,8.6,8.7,8.8,8.9,9,9.1,9.2,9.3,9.4,9.5,9.6,9.7,9.8,9.9,10,10.1,10.2,10.3,10.4,10.5,10.6,10.7,10.8,10.9,11,11.1,11.2,11.3,11.4,11.5,11.6,11.7,11.8,11.9,12,12.1,12.2,12.3,12.4,12.5,12.6,12.7,12.8,12.9,13,13.1,13.2,13.3,13.4,13.5,13.6,13.7,13.8,13.9,14,14.1,14.2,14.3,14.4,14.5,14.6,14.7,14.8,14.9,15,15.1,15.2,15.3,15.4,15.5,15.6,15.7,15.8,15.9,16,16.1,16.2,16.3,16.4,16.5,16.6,16.7,16.8,16.9,17,17.1,17.2,17.3,17.4,17.5,17.6,17.7,17.8,17.9,18,18.1,18.2,18.3,18.4,18.5,18.6,18.7,18.8,18.9,19,19.1,19.2,19.3,19.4,19.5,19.6,19.7,19.8,19.9,20};
            //波高数组
            float n[] = {0,0.0998334,0.198669,0.29552,0.389418,0.479426,0.564642,0.644218,0.717356,0.783327,0.841471,0.891207,0.932039,0.963558,0.98545,0.997495,0.999574,0.991665,0.973848,0.9463,0.909297,0.863209,0.808496,0.745705,0.675463,0.598472,0.515501,0.42738,0.334988,0.239249,0.14112,0.0415807,-0.0583741,-0.157746,-0.255541,-0.350783,-0.44252,-0.529836,-0.611858,-0.687766,-0.756802,-0.818277,-0.871576,-0.916166,-0.951602,-0.97753,-0.993691,-0.999923,-0.996165,-0.982453,-0.958924,-0.925815,-0.883455,-0.832267,-0.772765,-0.70554,-0.631267,-0.550686,-0.464602,-0.373877,-0.279415,-0.182163,-0.0830894,0.0168139,0.116549,0.21512,0.311541,0.40485,0.494113,0.57844,0.656987,0.728969,0.793668,0.850437,0.898708,0.938,0.96792,0.988168,0.998543,0.998941,0.989358,0.96989,0.940731,0.902172,0.854599,0.798487,0.734397,0.662969,0.584917,0.501021,0.412118,0.319098,0.22289,0.124454,0.0247754,-0.0751511,-0.174327,-0.271761,-0.366479,-0.457536,-0.544021,-0.625071,-0.699875,-0.767686,-0.827826,-0.879696,-0.922775,-0.956635,-0.980936,-0.995436,-0.99999,-0.994553,-0.979178,-0.954019,-0.919329,-0.875452,-0.822829,-0.761984,-0.693525,-0.618137,-0.536573,-0.449647,-0.358229,-0.263232,-0.165604,-0.0663219,0.033623,0.133232,0.23151,0.327474,0.420167,0.508661,0.592074,0.66957,0.740376,0.803784,0.859162,0.905955,0.943696,0.972008,0.990607,0.999309,0.998027,0.986772,0.965658,0.934895,0.894791,0.845747,0.788252,0.722881,0.650288,0.571197,0.486399,0.396741,0.303118,0.206467,0.107754,0.00796318,-0.0919069,-0.190859,-0.287903,-0.382071,-0.472422,-0.558052,-0.638107,-0.711785,-0.778352,-0.837142,-0.887567,-0.929124,-0.961397,-0.984065,-0.9969,-0.999774,-0.992659,-0.975626,-0.948844,-0.912582,-0.867202,-0.813157,-0.750987,-0.681314,-0.604833,-0.522309,-0.434566,-0.342481,-0.246974,-0.148999,-0.0495356,0.0504227,0.149877,0.247834,0.343315,0.435365,0.523066,0.60554,0.681964,0.751573,0.813674,0.867644,0.912945};
            float eta;
            for (int i = 0;i <= 200;i++)
            {
                if (t >= time[i] && t < time[i+1])
                {
                    //插值得到新的波高
                    eta = (n[i] + (t-time[i]) * (n[i+1]-n[i]) / (time[i+1]-time[i]));            
                }
            }
           
                forAll(Cf, faceI)
                {
                    const scalar y = Cf[faceI].y();
                    if (y <= eta )
                    {
                        field[faceI] = 1.;
                    }
                }
        #};
    }
// ************************************************************************* //

修改后能够正常编译,setFields能够正常设置,问题出在了速度文件上
以下是原算例U文件

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

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
left
    {
        type            codedFixedValue;
        value           uniform (0 0 0);
        redirectType    wave_velocity;

        code
        #{
            const fvPatch& boundaryPatch = patch(); 
            const vectorField& Cf = boundaryPatch.Cf(); 
            vectorField& field = *this; 

            field = patchInternalField(); 	//take value from initialization

            const scalar l = 5.;
            const scalar A = 0.1;
            const scalar gs = 9.81;
            const vector gv(0.,9.81,0);

            const scalar t = this->db().time().timeOutputValue();

            scalar k = 2.0*3.14159/l;
            scalar w = sqrt(k*mag(gv));

            forAll(Cf, faceI)
            {
                const scalar y = Cf[faceI].y();

                if (y <= A*cos(-w*t)+0.5*k*A*A*cos(2*(-w*t)) )
                {
                    field[faceI] = vector(A*w*exp(k*y)*cos(-w*t), A*w*exp(k*y)*sin(-w*t), 0);
                }
            }
        #};
    }

按照我的理解,应该同样将公式改为对应的波速,然后插值
但是我这里的数据只有时间和波高,没有速度,想请教大佬这种情况该采用哪种边界条件呢?

我想到一些思路:

  1. 类似造随机波的思路,把时间序列做FFT,分解为很多的线性正弦波,然后把线性波的速度解析解合并叠加,产生上述时间序列的入口速度,从而可解决该问题。

  2. 不使用VOF而使用动网格进行波浪数值模拟,这样只需在波面上给定近似压力边界条件p=\rho g \eta ,入口边界速度直接解出而不用指定。只是个想法,我自己没有试过

  3. 尝试现有波浪openfoam包:wave2foam,IHfoam、、、看看有没有直接输入时间序列的。

1 个赞

非常感谢!我试着做一下