CFD中文网

lizhisongsjtu

网格用blockMesh生成，纯流体，没有结构物。
对沿线分布数据的时间分辨率要求很高，如果用paraview中导出线上的数据，时间上的分辨率不够。
也许，用topoSet选取局部网格，大概也可以

lizhisongsjtu

这是其中的两个位置的垂向分布。在Mesh4中，第二个位置的数据就出现了缺损。目前搞不清楚怎么回事。这块数据缺失，很影响对计算结果的解读。

lizhisongsjtu

对某溃坝算力，用sampling沿垂向取密度，获取自由面的位置。

在做网格独立性验证时，突然发现某些位置的数据，存在数据异常，甚至数据缺失。

目前已经采取的改进措施：（1）改变并行方法（2）略微修改x坐标

这两个方法都没有效果。

哪位遇到过类似的问题，如何解决？

Lines
{
    type            sets;
    libs            (sampling);
    enabled         true;
    writeControl    adjustableRunTime;
    writeInterval   0.01;

    interpolationScheme cellPoint;
    setFormat       raw;
    fixedLocations  false;

    fields
    (
        rho
    );

    sets
    (
        line01
        {
            type    uniform;
            axis    xyz;
            start   (2.002 0 -0.24);
            end     (2.002 0  0.04);
            nPoints 1000;
        }
        line02
        {
            type    uniform;
            axis    xyz;
            start   (2.502 0 -0.24);
            end     (2.502 0  0.04);
            nPoints 1000;
        }
        line03
        {
            type    uniform;
            axis    xyz;
            start   (3.002 0 -0.24);
            end     (3.002 0  0.04);
            nPoints 1000;
        }
        line04
        {
            type    uniform;
            axis    xyz;
            start   (3.502 0 -0.24);
            end     (3.502 0  0.04);
            nPoints 1000;
        }

        line05
        {
            type    uniform;
            axis    xyz;
            start   (4.0 0 -0.24);
            end     (4.0 0  0.04);
            nPoints 1000;
        }

        line06
        {
            type    uniform;
            axis    xyz;
            start   (4.5 0 -0.24);
            end     (4.5 0  0.04);
            nPoints 1000;
        }
        line07
        {
            type    uniform;
            axis    xyz;
            start   (4.65 0 -0.24);
            end     (4.65 0  0.04);
            nPoints 1000;
        }
    );
}

结果
line04_rho.xy

lizhisongsjtu

可能是版本不同的问题，在OF2206中，添加缓存没有成功。

用另一种方法，算是曲线完成任务了。

在H文件中，添加 volScalarField NN;，在C文件初始化的列表中，加入

 NN
(
    IOobject
    (
        "NN",
        this->runTime_.timeName(),
        this->mesh_,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
     ),
    this->mesh_
),

然后在恰当的位置，完成赋值 NN = N2 即可。

不过有一个缺陷，需要再初始文件夹中设定 NN 的初始文件。

(先解决问题，麻烦就麻烦点吧)

lizhisongsjtu

对SST湍流模型作修正，引入了浮力项N2和流线弯曲修正项Fr1。希望在输出k、omega的时候，同步输出N2和Fr1.请问要怎么写？

尝试了N2.write();，这个语句把每个dt的结果全都输出了。

请问要怎么改写？

    // Calculate the Brunt-Vaisala frequency
    volScalarField N2 = gField_&fvc::grad(rho1)/rho1;
    N2.write(); // this exports data every time setp.
    //curvature correction
    const scalar cr1_ = 1.0;
    const scalar cr2_ = 2.0;
    const scalar cr3_ = 1.0;
    volScalarField Fr1    
    (
        max
    	(
            min
            (
                (scalar(1.0) + cr1_)*2.0*rStar /
                (scalar(1)+rStar) *
                (scalar(1.0) - cr3_*atan(cr2_*rTilda))-
                cr1_
                ,
                scalar(1.25)
            )
            ,
			scalar(0.0)
        )
    );

    // Turbulent kinetic energy equation
    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(alpha, rho1, k_)
      + fvm::div(rhoPhi, k_)
      - fvm::laplacian(alpha*rho1*DkEff(F1), k_)
     ==
        alpha()*rho1()*Pk(G) * Fr1() // curvature correction
      - fvm::SuSp((2.0/3.0)*alpha()*rho1()*divU, k_) 
      - fvm::Sp(alpha()*rho1()*epsilonByk(F1, F23), k_)
      - fvm::Sp(nut*alphaBS_*N2*rho1/max(k_,this->kMin_),k_) // Buoyancy production term

      + rho1*kSource()
      + fvOptions(alpha, rho1, k_)
    );

lizhisongsjtu

解决方案：

在of8及同系列版本中，单相流和多相流均使用codeAddSup定义源项。
在of2206及同系列版本中，单相流使用codeAddSup定义源项，多相流使用codeAddSupRho定义源项。

lizhisongsjtu

@李东岳我想确认一下，7楼的代码真的在2206上运行成功了吗？谢谢

附：of2206、of8和of12的日志。of2206失败。of8直接可以用7楼的代码，of12需放入fvModels中并略作修改。

log.interFoam.2206

log.interFoam.8

log.foamRun.12

lizhisongsjtu

感谢回复解惑。

1楼、2楼和4楼，即代码和算例，是在of2206上测试的，失败。

下面的一段代码，在of2206上仍然失败。在of8上暂时可以运行，源项效果有待检验。此上！

codedSource
{
    type            vectorCodedSource;
    name            sourceTime;

    selectionMode   cellSet;
    cellSet         myUSourceCellsZoneSet;
    
    fields          
    (
        U
    );
    
    codeAddSup
    #{
        const vectorField& C = mesh_.C();
        const scalarField& V = mesh_.V();
        vectorField& Usource = eqn.source();
        forAll(C, i)
        {
            // Info << "cellI = " << i << nl;
            const scalar y = C[i][0];
            const scalar vol = V[i];
            const scalar fx = vol * 1.0;

            Usource[i] -= vector(1.0, 0, 0);
        }
        Info << "***codeAddSup***" << nl;
    #};

    codeConstrain
    #{
    #};

    codeCorrect
    #{
    #};
}

lizhisongsjtu

上传一下case文件

interFoam.sourceU.zip

lizhisongsjtu

addSup的构造函数

//- Explicit/implicit matrix contributions
        virtual void addSup
        (
            fvMatrix<Type>& eqn,
            const label fieldi
        );

        //- Explicit/implicit matrix contributions to compressible equation
        virtual void addSup
        (
            const volScalarField& rho,
            fvMatrix<Type>& eqn,
            const label fieldi
        );

lizhisongsjtu

用interFoam计算两相流，在计算域内选择一块区域作为源区域，施加动量源项。
在 system/fvOptions 中添加了 vectorCodedSource类型的源项，并参考了网站链接的方法设置了源项。

根据报错信息，初步确定是添加源项的量纲出了问题。请高手指点，如何解决这个问题

最核心的一句报错信息是

--> FOAM FATAL ERROR: (openfoam-2206)
Not implemented

    From virtual void Foam::fv::codedMomentumSourceFvOptionvectorSource::addSup(const volScalarField&, Foam::fvMatrix<Foam::Vector<double> >&, Foam::label)

附fvOptions代码

codedSource
{
    type            vectorCodedSource;
    active          true;
    name            codedMomentumSource;

    // selectionMode   all;
    selectionMode   cellSet;
    cellSet         myUSourceCellsSet;  // 可以用 topoSet 定义的 cellSet
    fields          
    (
        U
    );

    codeAddSup
    #{  
        // https://www.cfd-china.com/topic/3255/%E5%88%86%E5%8C%BA%E7%BD%91%E6%A0%BC%E7%9A%84%E4%B8%80%E5%B0%8F%E6%AE%B5%E4%BB%A3%E7%A0%81?_=1731583465684
        Info << "++++From codedSource in fvOptions" << nl;

        const label cellZoneID = mesh().cellZones().findZoneID("myUSourceCellsSet"); // find ID for the cellZone "BANANA"
        
		const cellZone& zone = mesh().cellZones()[cellZoneID];
		const cellZoneMesh& zoneMesh = zone.zoneMesh();
		const labelList& cellsZone = zoneMesh[cellZoneID]; //list of all cellZone cell ID's  
        
        const Time& time = mesh().time();
        const scalarField& V = mesh_.V(); // Cell volume
            
        Foam::Field<Foam::Vector<double> >& USource = eqn.source(); // get source from fvMatrix
		//scalarField & Udiag = eqn.diag(); // get diagonal of fvMatrix
                
        const scalarField& rho = mesh().lookupObject<scalarField>("rho"); // get density field
        
        const vector g (0, 0, -9.81);      // gravitational vector


        forAll(cellsZone,i)
        {
            const label celli = cellsZone[i];
        
            // USource[celli] -= rho[celli]*g/V[celli] * sin(time.value());
            USource[celli] = sin(time.value());
        }
    #};

    //  chatGPT 的解决方案 + CFD-online解决方案
    // codeAddSup
    // #{
    //     Info << "++++From codedSource in fvOptions" << nl;

    //     const scalarField& V = mesh_.V();
    //     vectorField& USource = eqn.source();
    //     const scalarField& cellx = mesh_.C().component(0) ; 
    //     const scalarField& celly = mesh_.C().component(1) ;
        
    //     forAll(USource, cellI)
    //     {
    //         Info << "cellI = " << cellI << nl;
    //         // scalar v1 = 512*( 2*pow(cellx[i],6) - 6*pow(cellx[i],5) + 7*pow(cellx[i],4)
    //         //                  - 4*pow(cellx[i],3) + pow(cellx[i],2))*pow(celly[i],7) ;
    //         // // .... and the other 9 scalars
            
    //         // // Please ADD to the source, don't overwrite it.
    //         // // Overwriting it means you don't care about source terms in 
    //         // // DESCRETISATION OPERATORS (laplacian, div ...)
    //         USource[cellI][0] += 1.0 * USource[cellI][1];
    //     }
    // #};

    codeInclude
    #{
        #include "fvCFD.H"
    #};

    codeCorrect
    #{
        // Use Correct to make adjustments to the U values after they get calculated
        Info << "++++codeCorrect" << nl;
    #};
    codeConstrain
    #{
        // Use Constrain to modify the U equation before it gets solved
        Info << "++++codeConstrain" << nl;
    #};
}

lizhisongsjtu

来填坑！

@zzcfd @李东岳

https://doi.org/10.1063/5.0202407

请多指教。

lizhisongsjtu

@队长别开枪高手。膜拜。

lizhisongsjtu

waves2Foam是海岸工程研究中非常常用的造波模块。安装waves2Foam经常遇到版本匹配和编译错误。这里提供一个新的解决方案，用户可一键安装解决。

解决方案可从本人github中下载，地址是 https://github.com/zhisongli/publicData，文件名是 waves2Foam2156-of2206-github.zip。

关于openfoam和waves2Foam匹配方案，推荐openfoam2206和waves2Foam2156。

编译错误多来自于Fortran程序中，实参和形参的类型不匹配，也有来自程序语句过长。
此处提供一个新的ThirdParty/settings，用户可一键安装。

一键安装脚本

# install necessary packages
sudo apt-get install libgsl-dev
sudo apt install subversion

# download OceanWave3D-Fortran90 and move it to ThirdParty
git clone https://github.com/boTerpPaulsen/OceanWave3D-Fortran90.git 
cp -r OceanWave3D-Fortran90 ThirdParty/OceanWave3D-Fortran90BK

# download waves2Foam
svn co http://svn.code.sf.net/p/openfoam-extend/svn/trunk/Breeder_1.6/other/waves2Foam

# go to waves2Foam and copy my Thirdparty.
cd waves2Foam
rm -rf ThirdParty
cp -f ../ThirdParty .
cd ..

# move waves2Foam to $WM_PROJECT_USER_DIR/applications/utilities/
mkdir -p $WM_PROJECT_USER_DIR/applications/utilities/
cp -r waves2Foam $WM_PROJECT_USER_DIR/applications/utilities/

# compile
cd $WM_PROJECT_USER_DIR/applications/utilities/waves2Foam

./Allwmake

lizhisongsjtu

原型：森林、草原发生火灾，关心货源周围的温度和空气流动。

基于buoyantPimpleFoam求解器设置了一个算例。设定底部中心温度，比如1300K，计算得到的速度场并不对称。
温度和速度分布

上传计算文件，期待有热传经验的道友解答疑惑。

lizhisongsjtu

解惑重任承担不了。说说个人看法吧

个人认为，旋转流动和剪切流动中的涡结构不一样。根据我的计算结果，Liutex在识别旋转流动中的涡，涡结构对等直面参数取值不敏感；对剪切流动中的涡，任何涡识别方法的结果对等涡面取值敏感。#58楼的图片可以佐证这个说法。

关于编辑、审稿人对Liutex的不中立的态度，咨询过业内大佬。可能和“学派”有关，当然也有新方法本身的原因。新事物的出现，总需要一定的时间才能被接受。

引用郭德纲老师的一句话“文无第一，武无第二”，引用邓爷爷的话“黑猫白猫，抓住老鼠就是好猫”。
至于哪个方法更好，能解决问题就好。

以上抛砖，期待引玉讨论。

lizhisongsjtu

@camellia openfoam自带动边界造波，推板运动算法来自Goring，目标波形来自Boussinesq。对0.2以上的孤立波，输入波形和输出波形并不对应。

openfoam自带的dirichelet孤立波造波方法，Boussinesq方法没问题，Grimshaw和McCowan造波程序有严重的bug。

我修正了Grimshaw和McCowan程序的bug，植入了稳定造波0.7的造波方法。等几个月可以开源出来给大家用。

lizhisongsjtu

这只是某个算例佐证了liutex算法的优势。个人感觉(主流)学界对liutex的态度有非常大的差异，有些大佬对liutex法的态度偏消极(审稿意见)。

lizhisongsjtu

用两种方法识别流场中的涡，第一种是lambda法，openfoam自带的无量纲方法，第二种是liutex方法。

123看出，改变阈值后，中心部位的漩涡变小；456看出中心漩涡并没有随阈值改变而明显改变。

根据S和Omega模的大小，把涡区分为剪切占优和旋转占优两种。123+456显示，阈值对剪切占优的涡结构的显示特征影响挺大，对旋转占优的涡结构影响并没那么明显。如此，凸显了liutex方法在旋转涡涡结构识别中的优势。

lizhisongsjtu

总结一下:

求任意矩阵的逆矩阵，OpenFOAM中有三种方法:

(1) LU-scalarSquareMatrix，实数、方阵；

(2) QR-RectangularMatrix<Type> or SquareMatrix<Type> ，实数/复数、方阵/非方阵；

(3) SVD-scalarRectangularMatrix，实数、方阵/非方阵。

CFD中文网

lizhisongsjtu

帖子