CFD中文网

16c

看了一篇论文，说是在openfoam上搞了个dem扩展，原文是这个:
OpenHFDIB-DEM.pdf

然后对应的代码库在github是这个：
https://github.com/techMathGroup/openHFDIB-DEM

表面上看起来好像挺有用的，但是毕竟没用过也不大知道这东西的弊端，有没有大佬尝试过这东西，评价一下给小弟开开眼界。

16c

@贾光普后来这个工作有什么进展么

讲道理我现在在做类似于这个玩意的东西，只不过我做的是气液+颗粒三相的东西，现在卡在筛选颗粒这块好久了。目前已经确定筛网用多孔介质做了，nasa他们有很久很久之前的阻力系数测试的数据，但是这个颗粒筛选我在想要不要用liggghts重新写一个带筛选功能的反射边界。

16c

然后又去pimpleFoam跑了一下，代码改成这样的：


momentumSource
{
     type            vectorCodedSource;
     active          yes;
     name            sourceTime;
        vectorCodedSourceCoeffs
        {
            selectionMode   all;
          //  cellZone        pZone;
            fields          (U);
            codeInclude
            #{
                
            #};
            codeCorrect
            #{
        //        Pout<< "**codeCorrect**" << endl;
            #};
            codeAddSup
            #{
       //         Pout<< "**codeAddSup**" << endl;
         
               //    const vectorField& C = mesh_.C();
                  const scalarField& V = mesh_.V();
                  vectorField& Usource = eqn.source();
                  const vectorField& U = mesh().lookupObject<volVectorField>("U");
                 // const scalarField& Rho = mesh().lookupObject<volScalarField>("rho");
                  const scalarField& nu = mesh().lookupObject<volScalarField>("nu");
                //  const scalarField& magU = mag(U);
                scalar A = 1e2;
                scalar B = 1e2;
              //  vector C(0,1e4,0);
               
                   
                   
                    forAll(V,i)
                    {
                        const scalar x = mesh_.C()[i][0];
                        const scalar y = mesh_.C()[i][1];
        
                        if(x <  0.5 && x > 0 && y <  0.5 && y >  0.45)

                       {  
                        Usource[i] +=  (1e-5 * A + mag(U[i])* B * 0.5 ) * U[i]* V[i];
                       
                          // Usource =  (A * U[i] + B * mag(U[i]) * U[i]) * V[i];       
                          //  Usource[i] += - C * V[i]; 
                       }
                    
                   
                    
                    }   
                // Info << "***codeAddSup***" << nl;     
                
                                                  
                               
            #};
            codeSetValue
            #{
               // Pout<< "**codeSetValue**" << endl;
            #};
            // Dummy entry. Make dependent on above to trigger recompilation
            code
            #{
                $codeInclude
                $codeCorrect
                $codeAddSup
                $codeSetValue
            #};
        }
        sourceTimeCoeffs
        {
            $vectorCodedSourceCoeffs;
        }
}

跑完的速度场是这样：

就完全处于A和B只能在1e2这个数量级，但凡再大一点就直接浮点溢出了。。。完全不知道是咋回事。。

16c

interFoam跑完了的图是这样的

16c

给源项公式改了一下，但是发现一个非常难受的事情，本身这个东西是仿照darcy-forchheimer写的，但是openfoam原生的模型的D和F可以提到1e10这么高的数量级，但是我这个写的感觉只能到1e2这样子，1e3开始就会浮点溢出，但是1e2这个水平只是能看出来这个特定区域确实有阻力存在，并不能实现极大多孔介质阻力并形成类固体区域的状态，所以就很好奇这个是什么情况，新改的代码是这样：

momentumSource
{
     type            vectorCodedSource;
     active          yes;
     name            sourceTime;
        vectorCodedSourceCoeffs
        {
            selectionMode   all;
          //  cellZone        pZone;
            fields          (U);
            codeInclude
            #{
                
            #};
            codeCorrect
            #{
        //        Pout<< "**codeCorrect**" << endl;
            #};
            codeAddSup
            #{
       //         Pout<< "**codeAddSup**" << endl;
         
               //    const vectorField& C = mesh_.C();
                  const scalarField& V = mesh_.V();
                  vectorField& Usource = eqn.source();
                  const vectorField& U = mesh().lookupObject<volVectorField>("U");
                  const scalarField& Rho = mesh().lookupObject<volScalarField>("rho");
                  const scalarField& nu = mesh().lookupObject<volScalarField>("nu");
                //  const scalarField& magU = mag(U);
                scalar A = 1e3;
                scalar B = 1e3;
              //  vector C(0,1e4,0);
               
                   
                   
                    forAll(V,i)
                    {
                        const scalar x = mesh_.C()[i][0];
                        const scalar y = mesh_.C()[i][1];
        
                        if(x <  0.5 && x > 0 && y <  0.5 && y >  0.45)

                       {  
                        Usource[i] +=  (nu[i] * A + mag(U[i])* B * 0.5 ) * Rho[i] * U[i]* V[i];
                       
                          // Usource =  (A * U[i] + B * mag(U[i]) * U[i]) * V[i];       
                          //  Usource[i] += - C * V[i]; 
                       }
                    
                   
                    
                    }   
                // Info << "***codeAddSup***" << nl;     
                
                                                  
                               
            #};
            codeSetValue
            #{
               // Pout<< "**codeSetValue**" << endl;
            #};
            // Dummy entry. Make dependent on above to trigger recompilation
            code
            #{
                $codeInclude
                $codeCorrect
                $codeAddSup
                $codeSetValue
            #};
        }
        sourceTimeCoeffs
        {
            $vectorCodedSourceCoeffs;
        }
}

16c

最近在用codedsource写个阻力源项，但是写完之后感觉效果都不太好，昨天用paraview看了一眼矢量图，发现用多孔介质模型的区域是没有速度矢量的，但是用codedsource的话从t=0开始就存在速度矢量，然后导致介质满天飞。所以openfoam是怎么实现让多孔介质模型明明作为一个动量源项但是不沾任何速度矢量的。。。

16c

@李东岳
这个东西我也很奇怪，就是我看fluid mechanics 101的多孔介质视频给出的多孔介质公式是

这个形式的，然后看动量方程里面各项基本都是力的形式，我之前就想当然在后面的公式里面多乘了个体积V,后来看of帮助文档里面写的是

我就又改成上面那个形式，但是感觉好像写完了之后无论哪个形式也不收敛。。。

16c

我在fvoptions里面想设置一个区域性的阻力源项，但是因为需要实现的结构比较特殊，所以没法用toposet去进行捕捉获取，然后知乎看见一个大佬用codedsource写的一个特定区域的动量源项合计自己试试去写一个这种，为了尝试设置的阻力源项可不可用我就toposet了一个区域合计试试，但是无论我的源项是纯数字还是仿darcy-forchheimer公式类型的他用interFoam跑都不收敛，找了好几天论坛视频啥的也找不到原因，vscode报错也只是给代码报错，偶尔代码看起来没错误但是跑起来就是不收敛。现在目前写的是如下这样的，因为回头要耦合dem，所以of版本是5.x。

momentumSource
{
     type            vectorCodedSource;
     active          yes;
     name            sourceTime;
        vectorCodedSourceCoeffs
        {
            selectionMode  cellZone;         // all;
            cellZone        pZone;
            fields          (U);
            codeInclude
            #{
                
            #};
            codeCorrect
            #{
        //        Pout<< "**codeCorrect**" << endl;
            #};
            codeAddSup
            #{
       //         Pout<< "**codeAddSup**" << endl;
         
               //    const vectorField& C = mesh_.C();
                const scalarField& V = mesh_.V();
                  vectorField& Usource = eqn.source();
                  const vectorField& U = mesh().lookupObject<volVectorField>("U");
                  const scalarField& Rho = mesh().lookupObject<volScalarField>("rho");
                  const scalarField& nu = mesh().lookupObject<volScalarField>("nu");
                  const scalarField& magU = mag(U);
                    scalar A = 10;
                    scalar B = 10;
                    forAll(V,i)
                    {
                    Usource -= (A * nu * Rho + B * Rho * magU / 2) * U ;                 
                    }  
                // Info << "***codeAddSup***" << nl;                                         
                               
            #};
            codeSetValue
            #{
               // Pout<< "**codeSetValue**" << endl;
            #};
            // Dummy entry. Make dependent on above to trigger recompilation
            code
            #{
                $codeInclude
                $codeCorrect
                $codeAddSup
                $codeSetValue
            #};
        }
        sourceTimeCoeffs
        {
            $vectorCodedSourceCoeffs;
        }
}

请各位大佬各位巨佬帮着指点一下哪里出现了问题。

16c

@李东岳我因为可能要做cfdem所以用的5.x跑的

5.x改了好多次网格尺寸和dt都在浮点溢出，后来看网上教程给piso改成pimple模式倒是不溢出了，但是这个结果就很离谱了

16c

@李东岳 test2.rar
改好了大佬您在帮着看一眼

16c

在openfoam里用interFoam算了一个很简单的cas，大概就是一个平板前面有个inlet,往里面注水，想看这个水流的情况，但是出现了非常离谱的状况，我先放上结果图
我设置的速度方向是(1 0 0),因此我觉得他不应该出现向上喷射的这个状态。
同样的我又在fluent里面简单跑了一下这个cas，边界条件设置都是相同的，结果是这样的：

今天着重的了解了一下这个p_rgh文件，大概可以理解为引入的一个方便求解的中间参量，等于原来的p-rhogh。
两者设置里面最唯一出现的差别可能就是压力，fluent的设置更加傻瓜式一点，直接对边界区域的相对压力给定值就好了，这个p_rgh原则上感觉不需要去关心具体是多少，但是总感觉这个东西好像跟不同的坐标有一定的关系，比方说坐标（ 1 1 0）的点和坐标（ 1 1 1 ）的点在初场设置时想让他们的静压或者说P文件中的压力为0的时候，需要设置在的p_rgh中的值应该是不一样的。
感觉出现这个问题也可能是我openfoam这边设置的问题。所以也得请求各位大佬帮着看一下，openfoam的设置文件我放在下面了。
test2.zip

因为文件稍微有点大我就给网格删了，不过都是均匀网格所以画起来也不麻烦，直接blockMesh就好，不过需要topoSet+createPatch一下来创造一个inlet来着

16c

我自己理解这个这个p_rgh应该是其他求解器0文件下的p-rg·h，但是interFoam似乎只能读取p_rgh而不能读取p文件，但是在对边界和整个初场只已知p的时候p_rgh怎么设置，尤其是inlet和outlet这种在某些时候需要给定一些数值的时候。

16c

最近在跑几个非常简单的case，发现网格画的比较精细的时候反而经常库朗数爆炸，网格画的相对粗糙的时候反而跑的杠杠的。是与电脑配置有关还是说求解器有个适配的敏感度啥的，我是一点头绪也没有。
而且报错之后浮点溢出的东西在我看起来就像是一坨乱码，unset FOAM_SIGFPE之后也没有具体的错误原因，这个报错怎么看有没有啥关键字或者资料啥的，还得请各位大佬提点提点。

16c

@学流体的小明在 toposet之后出现空集。。。中说：

感觉是topoSet选择网格，提取各种面的时候有些问题了，我没用过这方面的功能，不太了解。
你如果只是想把边界面提取出来的话，在blockMesh当中就可以设置吧。

@李东岳在 toposet之后出现空集。。。中说：

豁牙子注意一下toposet点位置的精度。

感谢两位大佬提点，豁牙子问题解决了...

16c

@学流体的小明在 toposet之后出现空集。。。中说：

你blockMesh文件convertToMeters 0.001;
实际就是把下面的坐标乘以0.001画网格，画出来的计算域是(0,0,0)到(2,0.5,0.05)。
而你toposet直接从(0,0,10)开始画box，当然找不到了。

大哥按你说的改了topoSetDric , box的坐标和网格尺寸我也确定了是能对应上的，但是createPatch并ParaFoam之后出现了inlet这个面变成了豁牙子。。。

改之后的topoSet是这样的：

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5                                     |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      topoSetDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

actions
(
   {
        name    porous;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box     (0 0 10e-3) (2000e-3 500e-3 12e-3);
        }
    }

   
   
   {
        name    porousZone;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        sourceInfo
        {
            set porous;
        }
    }
        

   
    {
        name    inletA;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box (0 0 12e-3) (4e-3 500e-3 22e-3);
        }
    }

    {
        name    inletB;
        type    faceSet;
        action  new;
        source  cellToFace;
        sourceInfo
        {
            set     inletA;
            option  all;
        }
    } 

    {
        name    inletB;
        type    faceSet;
        action  subset;
        source  patchToFace;
        sourceInfo
        {
            name walls;
        }
    }

    {
        name    outletA;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box (1.996 0 12e-3) (2 500e-3 50e-3);
        }
    }

    {
        name    outletB;
        type    faceSet;
        action  new;
        source  cellToFace;
        sourceInfo
        {
            set     outletA;
            option  all;
        }
    } 

    {
        name    outletB;
        type    faceSet;
        action  subset;
        source  patchToFace;
        sourceInfo
        {
            name walls;
        }
    }
        
);

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

createPatch是这样的：

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5                                     |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      createPatchDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Do a synchronisation of coupled points after creation of any patches.
// Note: this does not work with points that are on multiple coupled patches
//       with transformations (i.e. cyclics).
pointSync false;
// writeCyclicMatch  false;
// Patches to create.
patches
(
    {
        name inlet;
        patchInfo
        {
            type            patch;
        }
        constructFrom set;
        set inletB;
    }

    {
        name outlet;
        patchInfo
        {
            type            patch;
        }
        constructFrom set;
        set outletB;
    }
);

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

这个豁牙子问题我之前跟着阿B上田东的视频做案例的时候也出现了不知道是什么玩意而且case跑完之后他明显是个实际存在的数值上的豁牙子，不是paraview显示问题

16c

toposet之后polymesh下set文件夹的文件内全部显示是空集，终端显示识别了你的范围是XXX但是这个范围内狗der没有，我是真不到哪里出问题了。
toposetdric代码是这样的：

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5                                     |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      topoSetDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

actions
(
   {
        name    porous;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box     (0 0 10) (2000 500 15);
        }
    }

   
   /* 
   {
        name    porousZone;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        sourceInfo
        {
            set porous;
        }
    }
        */

   
    {
        name    inletA;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box (0 0 15) (0.1 500 25);
        }
    }

    {
        name    inletB;
        type    faceSet;
        action  new;
        source  cellToFace;
        sourceInfo
        {
            set     inletA;
            option  all;
        }
    } 

    {
        name    inletB;
        type    faceSet;
        action  subset;
        source  patchToFace;
        sourceInfo
        {
            name walls;
        }
    }

    {
        name    outletA;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box (1999.9 0 15) (2000 500 50);
        }
    }

    {
        name    outletB;
        type    faceSet;
        action  new;
        source  cellToFace;
        sourceInfo
        {
            set     outletA;
            option  all;
        }
    } 

    {
        name    outletB;
        type    faceSet;
        action  subset;
        source  patchToFace;
        sourceInfo
        {
            name walls;
        }
    }
        
);

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

终端是这样说的：

openfoam@openfoam-virtual-machine:~/OpenFOAM/openfoam-5.x/run/1$ topoSet
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5.x                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 5.x-7f7d351b741b
Exec   : topoSet
Date   : Nov 25 2024
Time   : 09:45:57
Host   : "openfoam-virtual-machine"
PID    : 181025
I/O    : uncollated
Case   : /home/openfoam/OpenFOAM/openfoam-5.x/run/1
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Reading topoSetDict

Time = 0
    mesh not changed.
Created cellZoneSet porous
    Applying source boxToCell
    Adding cells with center within boxes 1((0 0 10) (2000 500 15))
    cellZoneSet porous now size 0
Created cellSet inletA
    Applying source boxToCell
    Adding cells with center within boxes 1((0 0 15) (0.1 500 25))
    cellSet inletA now size 0
Created faceSet inletB
    Applying source cellToFace
    Adding faces according to cellSet inletA ...
    faceSet inletB now size 0
Read set faceSet inletB with size 0
    Applying source patchToFace
    Adding all faces of patch walls ...
    Found matching patch walls with 31250 faces.
    faceSet inletB now size 0
Created cellSet outletA
    Applying source boxToCell
    Adding cells with center within boxes 1((1999.9 0 15) (2000 500 50))
    cellSet outletA now size 0
Created faceSet outletB
    Applying source cellToFace
    Adding faces according to cellSet outletA ...
    faceSet outletB now size 0
Read set faceSet outletB with size 0
    Applying source patchToFace
    Adding all faces of patch walls ...
    Found matching patch walls with 31250 faces.
    faceSet outletB now size 0
End

blockMesh是这样的：

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5                                     |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
convertToMeters 0.001;
vertices
(
    (0 0 0)                 // 0
    (2000 0 0)               // 1
    (2000 500 0)            // 2
    (0 500 0)              // 3

   /* 
    (0 0 15)                // 4
    (2000 0 15)              // 5
    (2000 500 15)           // 6
    (0 500 15)             // 7
    */

    (0 0 50)                // 4
    (2000 0 50)              // 5
    (2000 500 50)            // 6
    (0 500 50)             // 7   

    // (0 0 25)                // 12
    // (0 500 25)              // 13
);
blocks
(
    hex (0 1 2 3 4 5 6 7) (500 125 25) simpleGrading (1 1 1)
    // hex (4 5 6 7 8 9 10 11) (500 125 35) simpleGrading (1 1 1) 
    
);

edges
(
);

boundary
(

    walls
    {
        type wall;
        faces
        (
            (0 1 5 4)
            (2 3 7 6)
            (0 4 7 3)
            (1 2 6 5)
        );
    }
   
    atmosphere
    {
        type patch;
        faces
        (
            (4 5 6 7)
            (0 1 2 3)
        );
    }
mergePatchPairs
(
);

我搞这个破玩意好几天了，一点头绪没有，还得求着各位大佬帮帮忙瞅一眼。。。。

CFD中文网

16c

帖子