使用openfoam但不采用贴体网格计算圆柱绕流的方法，请各位同学老师指点一下，谢谢
-
各位老师和同学好。之前一直使用课题组代码，最近在转向学习OpenFOAM，在做一个小东西来练习，想不采用贴体网格，而是浸入边界法或者cut cell的方法来计算圆柱绕流。
-
方法1：浸入边界法。拟采用论文1：Ji C, Munjiza A, Williams J J R. A novel iterative direct-forcing immersed boundary method and its finite volume applications[J]. Journal of Computational Physics, 2012, 231(4): 1797-1821.
-
方法2：cut cell方法。拟采用论文2：Lin P. A fixed-grid model for simulation of a moving body in free surface flows[J]. Computers & fluids, 2007, 36(3): 549-561.
首先采用方法1：论文1中的方法是一个direct-forcing immersed boundary method，固体的作用通过体积力施加给NS方程，固体内部相当于存在虚拟流体，虚拟流体和固体边界外流体刚好在固体边界处满足无滑移边界条件。具体的实现在论文1中是基于两步映射法，在第一步求得中间速度后将流体速度插值到固体点边界，固体速度与流体插值过来的速度应相等满足无滑移，则可计算出速度差得额外体积力f，并将该体积力通过分布函数施加给邻近的影响的流体网格，更新中间速度，然后采用中间速度求压力泊松方程得到该步压强再更新最终速度。(注：论文中增加了一个迭代循环，但是循环1次也即等同于传统直接应力的浸入边界法)
那么在OpenFOAM中，采用的piso或pimple算法。我应该在每一步将额外体积力f施加在什么部分呢？
首先我试了采用上一步计算得到的流体速度来插值到固体边界先求出f，然后施加在Ueqn里//do something to calculate f using U of last step, //then : tmp<fvVectorMatrix> tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevSigma(U) == fvOptions(U) + f );用这种方法计算无法收敛。
然后我在UEqn之后，采用该中间速度算出f，施加在momentumPredictor的部分：tmp<fvVectorMatrix> tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevSigma(U) == fvOptions(U) ); fvVectorMatrix& UEqn = tUEqn.ref(); UEqn.relax(); fvOptions.constrain(UEqn); //do something to calculate f using the intermediate U , //then : if (pimple.momentumPredictor()) { solve(UEqn == -fvc::grad(p) + f); fvOptions.correct(U); }这种方法计算是可以收敛的，而且可以正常算出涡，我计算的是一个Re=100的圆柱绕流，但是这个涡是不对的，和正确位置比滞后了很多，而且拖曳力系数是正确拖曳力系数的好几倍。并且这种施加方法，结果会和piso或pimple算法的correct步数相关了(nCorrectors)。correct步数次数在openfoam中默认2次，可以增加次数，当correct次数改变时候，我计算得到的拖曳力系数成倍增长，涡的形态也发生变化。
这是Re=100圆柱绕流正确的形态，为我使用课题组代码计算得到的：
这是当前我采用Openfoam利用我上述施加f在momentumPredictor()部分的方法算出的结果：
涡滞后了很多，而且拖曳力大了几倍甚至几十倍，随着correct次数的增加而增加。在谷歌一些相关资料后，很多方法是将f施加在Peqn里，看起来类似interFoam里重力施加的方式，初步尝试总是不收敛，暂时我应该遗漏了某些细节。
除此之外，也尝试了方法2的cut cell方法：在论文2采用一个系数表征固体，流体和流固交界cell。固体部分系数为0，流体系数为1，相交部分为0-1，这类似VOF的方法。论文2核心其实是在当固体运动时采用一个相对速度来修正连续性方程保证质量守恒。我假设最最简单的情况，忽略部分切割(系数在0-1之间)的网格，且固体不动，就暂时不需要论文2对连续性方程的修正部分了。同样计算Re=100的圆柱绕流。按照论文2的方法，完全被固体占据的网格（系数=0）不参与计算要跳过（如果不跳过，会出现除以0的情况）。为了实现固体完全占据的部分不参与计算且速度设为0（有些等同于于扣掉网格了），我采用了这个帖子中@马乔 的方法：https://www.cfd-china.com/topic/3138/openfoam有方法能够使一部分网格不参与计算吗/28
\*---------------------------------------------------------------------------*/ #include "fvCFD.H" #include "pisoControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { //Info<< "\nStarting time loop\n" << endl; #include "setRootCaseLists.H" #include "createTime.H" #include "createMesh.H" pisoControl piso(mesh); #include "createFields.H" #include "initContinuityErrs.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; DynamicList<label> lowerFaceCells(mesh.C().size()); DynamicList<label> upperFaceCells(mesh.C().size()); //mesh LDU address const labelUList& lowAddr = mesh.lduAddr().lowerAddr(); const labelUList& upAddr = mesh.lduAddr().upperAddr(); //label refCell = 4; //labelList cellBlock (20); //DynamicList<label> cellBlock; if(mesh.cellZones().size()) { const labelList& cellBlock = mesh.cellZones()[0]; Info << "cell block size: " << cellBlock.size() << endl; //Info << cellBlock << endl; //low efficient forAll(lowAddr, faceI) { for(auto cellI : cellBlock) { if(lowAddr[faceI] == cellI) { lowerFaceCells.append(faceI); } if(upAddr[faceI] == cellI) { upperFaceCells.append(faceI); } } } } lowerFaceCells.shrink(); upperFaceCells.shrink(); Info << "low Cells size: " << lowerFaceCells.size() << ", low Cells size: " << upperFaceCells.size() << endl; while (runTime.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; #include "CourantNo.H" // Momentum predictor fvVectorMatrix UEqn ( fvm::ddt(U) + fvm::div(phi, U) - fvm::laplacian(nu, U) ); /*label refCell = 4; labelList refCells (1,refCell); vectorField refValues (1, vector(0,0,0)); UEqn.setValues(refCells, refValues);*/ //modify matrix if(mesh.cellZones().size()) { const cellZone& cellBlock = mesh.cellZones()[0]; //should select cellZone accor name scalarField& lower = UEqn.lower(); scalarField& upper = UEqn.upper(); scalarField& diag = UEqn.diag(); for(auto cellI : cellBlock) { diag[cellI] = 1.; UEqn.source()[cellI] = vector::zero; //vector(0.,0.,0.); you can set your deserved value here, or read from dict, ref. fvoptions } for(auto faceI : lowerFaceCells) { lower[faceI] = 0.; } for(auto faceI : upperFaceCells) { upper[faceI] = 0.; } } else { Info << "There is no cell selected" << endl; } if (piso.momentumPredictor()) { solve(UEqn == -fvc::grad(p)); } // --- PISO loop while (piso.correct()) { volScalarField rAU(1.0/UEqn.A()); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); surfaceScalarField phiHbyA ( "phiHbyA", fvc::flux(HbyA) + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi) ); adjustPhi(phiHbyA, U, p); // Update the pressure BCs to ensure flux consistency constrainPressure(p, U, phiHbyA, rAU); // Non-orthogonal pressure corrector loop while (piso.correctNonOrthogonal()) { // Pressure corrector fvScalarMatrix pEqn ( fvm::laplacian(rAU, p) == fvc::div(phiHbyA) ); pEqn.setReference(pRefCell, pRefValue); //modify p matrix if(mesh.cellZones().size()) { const cellZone& cellBlock = mesh.cellZones()[0]; //p matrix symmetry scalarField& offDiag = pEqn.hasUpper() ? pEqn.upper() : pEqn.lower(); scalarField& diag = pEqn.diag(); for(auto cellI : cellBlock) { diag[cellI] = 1.; pEqn.source()[cellI] = 0.; } //for(auto faceI : upperFaceCells) //should find faces again? //{ // offDiag[faceI] = 0.; //} const labelUList& offAddr = pEqn.hasUpper() ? pEqn.lduAddr().upperAddr() : pEqn.lduAddr().lowerAddr(); forAll(offAddr, faceI) { for(auto cellI : cellBlock) { if(offAddr[faceI] == cellI) { offDiag[faceI] = 0.; } } } } pEqn.solve(); if (piso.finalNonOrthogonalIter()) { phi = phiHbyA - pEqn.flux(); } } #include "continuityErrs.H" U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); } runTime.write(); //const vectorField& UCells = U.internalField(); //dimensionedVector v1 ("v1", dimVelocity, vector (1.0,0.0,0.0)); //U.internalField()[4] = v1; /*Info << U << nl << endl; Info << mesh.C() << nl << endl;*/ Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; }确实可以计算，但是会一开始出现一点点涡，同样也非常滞后看起来不正常。然后随着计算时间的增加，涡完全消失不产生绕流。不知是否在这种情况还需要对跳过cell处的边界做一些特殊处理？下面的图是刚计算了几秒的时候，看起来是正常的，算一段时间后产生涡，同样滞后。再算一段时间，涡消失完全没涡了。
小弟接触openfoam时间不长，希望各位老师和同学能帮忙给点建议或者方向。解决这个采用IBM或cut cell来计算圆柱绕流的问题。不采用贴体网格是因为我后续要做的东西更适合和IBM或cut cell相结合。谢谢大家！
-
-
-
@李东岳 李老师 是的 我尝试了体积力插到面上 代码是这样的 f是算好的。主要是三句：
surfaceScalarField phiIB(fvc::flux(rAtU*f)); fvScalarMatrix pEqn ( fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA) + fvc::div(phiIB) ); U = HbyA - rAtU*fvc::grad(p) + rAtU*f;即插值得到phiIB， 加到pEqn，然后下面是U
= HbyA - rAtU*fvc::grad(p) + rAtU*f， 这样做无法收敛，时间步会无限小下去。不知道错在哪了，这是详细的Peqn.C：volScalarField rAU(1.0/UEqn.A()); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); surfaceScalarField phiHbyA ( "phiHbyA", fvc::flux(HbyA) + MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf)) ); MRF.makeRelative(phiHbyA); if (p.needReference()) { fvc::makeRelative(phiHbyA, U); adjustPhi(phiHbyA, U, p); fvc::makeAbsolute(phiHbyA, U); } tmp<volScalarField> rAtU(rAU); if (pimple.consistent()) { rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU); phiHbyA += fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf(); HbyA -= (rAU - rAtU())*fvc::grad(p); } surfaceScalarField phiIB(fvc::flux(rAtU*f)); if (pimple.nCorrPiso() <= 1) { tUEqn.clear(); } // Update the pressure BCs to ensure flux consistency constrainPressure(p, U, phiHbyA, rAtU(), MRF); // Non-orthogonal pressure corrector loop while (pimple.correctNonOrthogonal()) { Info << "!herehrehrehrehrheh" << nl; fvScalarMatrix pEqn ( fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA) + fvc::div(phiIB) ); pEqn.setReference(pRefCell, pRefValue); pEqn.solve(); if (pimple.finalNonOrthogonalIter()) { phi = phiHbyA - pEqn.flux(); } } #include "continuityErrs.H" // Explicitly relax pressure for momentum corrector p.relax(); U = HbyA - rAtU*fvc::grad(p) + rAtU*f; U.correctBoundaryConditions(); fvOptions.correct(U); -
@李东岳
李老师 按您说的尝试了下 还是时间步无限减小surfaceScalarField phiIB(fvc::flux(rAtU*f)); ... fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA) + fvc::div(phiIB) ... phi = phiHbyA - pEqn.flux() + phiIB; ... U = HbyA - rAtU*fvc::grad(p) + rAtU*f;然后只能算步长0.02， 算2步到0.04秒就不收敛了，结果是错误的。全域一开始1m/s,入口1m/s,然后0.04s时速度和压强云图是这个样子，速度在圆柱一圈变得很大，到了20多米每秒，压强是一个很奇怪的完全一半一半对称的形状：
然后我也试着直接把f加到HbyA里。一样是不收敛的volVectorField HbyA(constrainHbyA(rAU*UEqn.H() + rAU*f, U, p));奇怪了，感觉加个f而已应该很简单才对呀。
-
