CFD中文网

把粒子按x坐标排序，i是序号，形成xorder(i,1)存x坐标，xorder(i,2)存粒子编号；
同理，把粒子按y、z坐标排序，形成yorder(i,1..2)，zorder(i,1..2)
计算网格顶点的xmin，xmax，ymin，ymax，zmin，zmax

loop:循环网格{
      
       找到xmin，xmax在xorder(i,1)中的位置
       同理，处理ymin，ymax，zmin，zmax
       new 一个数组 possible(i)
     
       loop: 处理上面缩小的范围{
              把以上范围内的粒子编号xorder(i,2)，yorder(i,2),
              zorder(i,2)插入排 序到possible(i)
       }

       loop: 确定真正在网格中的粒子{
             possible(i)中有连续三个相同编号，判断是否真的在网格内,
             如果真在，就存起来
       }
}

至于每个网格里的粒子编号怎么存，搞一个数组叫data(i)一个接一个地存上面伪代码得到的编号，另一个数组叫position(j)存单元j在data(i)中的起始位置
不知道问题理解对不对，抛个砖试试看

@dzw05 学习了，div(phid,p) 确实不用 bounded 为好
@东岳 把 z01~02 后面的zip去掉，然后对着 rhoSimple.zip 选择解压

非气动专业，尝试和大佬讨论，想问一下：
（1）div(phid,p) 的格式为什么不用 bounded Gauss upwind 呢？经验来看，bounded格式稳定很多；
（2）grad用最小二乘，以及interpolation用迎风也是为了数值稳定性吧？但貌似当前算例都用中心差分也可以，所以为了精度给改了。
按照以上两项修改试着跑了几万步，不过不是特别收敛，连续性方程不满足，升阻力系数还在震荡。是不是这个物理问题本身是非定常的，或者是网格太粗糙呢？
还有就是挺好奇这个网格使用sHM画的吗？像这类混合网格如何处理trailing edge这种特变尖锐的地方？

可能你的操作系统位数和别人不一样。
试下装一个和你目前系统位数不一样的虚拟机（或者用超算、别人的电脑）
然后把case搬到这台机器上，在controlDict 里面writeFormat 改成 ascii，终端键入foamFormatConvert看看能不能转换

抱歉说漏了，rho、rho.water、rho.air这三个文件都要放进0文件夹。funkySetFields是能够使用函数表达式的，用来设置场内分布；边界上的值可以用funkySetBoundaryField来设置。在system下添加funkySetFieldsDict和funkySetBoundaryDict来使用。里面的内容，随便举个例子：

//funkySetFields
expressions
 (
    pressureWater
    {
     field p_rgh; //需要设置的物理量
     expression "100000 - 99800*pos().y";  //这行写你要的函数
     condition  "(pos().x <= 5) && (pos().x >= 0)"; //起作用的范围、可以注释掉，不影响使用
     keepPatches 1; //填1就好，要不然可能把边界条件改成zeroGradient
    }
 );

//funkySetBoundaryField
pressureWater
 {
    field p_rgh;
    expressions
    (
        {
            target value;
            patchName INLET; //选择边界面
            //variables "maxY=max(pts().y);";  //可以注释掉，不影响使用
            expression "100000 - 99800*pos().y"; //这行写你要的函数
        }
         //可以再expression下同时写多个边界的设置
        {
            target value;
            patchName OUTLET;
            //variables "maxY=max(pts().y);";
            expression "100000 - 99800*pos().y";
        }
    );
}

更多细节参考这里链接文本，以及swak4foam的tutorials
另外，你的这个算例里面压力也应当要设置，需要和密度相匹配。

1. 终端输入mpirun -np n pvserver （n根据情况调整），保持该终端不关闭
2. 打开paraFoam，file -> connect -> add server : 其中，Name随意， Server type 用 client/server，Host 用 localhost，port 默认即可。
3. 加载 .foam文件

感觉这种方法要比reconstructPar快

compressibleInterFoam读入rho的方式是READ_IF_PRESENT，也就是说存在rho文件的话，就会被读进去。所以可以在0文件夹下自己添加一个rho，然后用funkySetField 设置你要的初始场。
至于“如何设置边界密度随空间变化”，compressibleInterFoam是基于压力的求解器，也就是密度由压力决定。所以只要把p_rgh设置对了就可以吧？

@maoyanjun_dut 非常感谢你的工作！另外想再请教个问题，如果要把overset用到compressibleInterFoam这类可压缩求解器里面的话，可压缩带来的dgdt，pEqnComp应该怎么处理？

谢谢东岳老师考古！讲解非常清楚！

如果是除0的话，应该报 #2 Foam::divide at ??:? 这样的错误。
你的这个错误像是缺少了自己写的库文件。FOAM_USER_APPBIN下能找到吗？试着重新编译一下？

举个相关的例子，比如我们要给a赋值，公式是a等于b乘以c，那么有两种方式：
第一种，最普通的：

a = b*c;

还有第二种，调用lookUpObject：

a = b*b.db().lookUpObject("c");

这样，lookUpObject可以从已经注册的变量中去寻找需要的那一个，如果找不到就会报你看到的那种错误。
可想而知，如果用的是第一种方式，没有定义变量c就会导致编译时提示错误。而如果用的是第二种方式，就算没有定义变量c，编译时也不会报错，但计算时就会出错。
totalPressureFvPatchScalarField.C 中虽然调用的是lookupPatchField<surfaceScalarField, scalar>("xxxxxx")，但情况是类似的

把压力的边界条件贴出了看看吧

某个边界的压力边条用了totalPressure吧？你的代码里面是不是没有psi这个变量？

@李东岳
谢谢东岳老师！明白是什么意思了！

1.它对矩阵进行了哪些操作？
2.比如，在fvSolution中设置：

        solver          smoothSolver;
        smoother        GaussSeidel;

那么这和教科书上最普通的GaussSeidel迭代有什么区别呢？

按照我浅薄的理解，有界往往是通过差分格式（如TVD格式），或者对流项特殊处理方法（如FCT）来完成的。所以没有想通 alpha1*divU 究竟对保证有界究竟有何贡献？如何用数学语言分析它的作用？
另外，问题的来源是这样的：我结合compressibleInterFoam和interPhaseChangeFoam改编了一个包括水、水蒸气、不可凝结气体的三相带空化的可压缩求解器。目前有一个困扰了我很长时间的bug，如果存在不可凝结气体设置为fixedValue uniform 1的边界（如带有自由面问题的最上方的大气），水中本应该发生空化的地方莫名其妙产生不可凝结气体。很肯定空化模型没有写错，所以认为问题很可能出现在alphaEqn的源项处理上不得当，不得不深挖一下其中的细节。
如果前辈对引入divU 的作用解释得再详细一点就好了！求甩我一脸文献！ @李东岳
欢迎研究类似问题的朋友们一起来讨论！

@李东岳 谢谢东岳老师的回复！
但是为什么说原方程的对流项不能保证有界呢？

第一相的体积分数输运方程为：

但实际求解的时候要利用：

将方程改写为：

我的疑问是，为什么要将alpha1*divU引入原方程中呢？它的作用如何用数学语言表达？

像interFoam、interPhaseChangeFoam采用了均相流模型的多相流求解器必然要对各相粘性系数取平均。而粘性项在Uqen中用turbulence->divDevRhoReff 来调用。我想看看 turbulence->divDevRhoReff 究竟是经历哪些成员函数的调用过程，才求得各相的粘性系数的平均值。但由于湍流类采用了Run Time Selection 机制，代码比较复杂一些，所以没搞懂。
大家能否以laminar湍流模型为例，解答一下这个问题？

@影川风
twoPhaseMixtureThermo.C里面确实有，但是，，

void Foam::twoPhaseMixtureThermo::correct()
{
    thermo1_->he() = thermo1_->he(p_, T_);
    thermo1_->correct();

    thermo2_->he() = thermo2_->he(p_, T_);
    thermo2_->correct();

    psi_ = alpha1()*thermo1_->psi() + alpha2()*thermo2_->psi();
    mu_ = alpha1()*thermo1_->mu() + alpha2()*thermo2_->mu();
    alpha_ = alpha1()*thermo1_->alpha() + alpha2()*thermo2_->alpha();
}

这里的thermo1_是指向rhoThermo类的，真的可以调用twoPhaseMixtureThermo类里面的函数吗？我在twoPhaseMixtureThermo::he()里面加了一句Info，运行求解器的时候貌似没有在终端里看到那句Info。c++基础不牢求指导