OpenFOAM有方法能够使一部分网格不参与计算吗？

Samuel-Tu

想在计算过程中使得一部分网格不参与计算，有没有啥近似的方法呢？我想的是对那部分网格赋予一个很大的粘度，使他们近似于固体，不知道可不可行。

七辆战车

把网格抠除？

dzw05

@Samuel-Tu 试过，可行。

Samuel-Tu

@七辆战车对，就是想达到扣除的效果

Samuel-Tu

@dzw05 感谢

我是河滩

请问你的网格是怎么扣除的，嵌套网格有挖洞的操作，洞单元就不参与计算，你了解这个吗？

Samuel-Tu

@我是河滩不了解.......

林之流风

@dzw05 请问您具体是如何处理的？

dzw05

@林之流风我计算的是两相流问题，把感兴趣区域中某一相的粘度设的特别大（比如1e5），该相就相当于凝固了。

东岳

@Samuel-Tu 在 OpenFOAM有方法能够使一部分网格不参与计算吗？中说：

想在计算过程中使得一部分网格不参与计算，有没有啥近似的方法呢？我想的是对那部分网格赋予一个很大的粘度，使他们近似于固体，不知道可不可行。

这是要模拟什么工况？凝固？

Samuel-Tu

@东岳类似吧，单向流问题，想随着流场计算，我可以控制流体域内的固体形状发展

东岳

看起来你还需要考虑边界条件的问题

Samuel-Tu

@东岳嗯嗯，就是考虑怎么在流固交界面上把壁面剪切力修正一下

cccrrryyy

借鉴一下沉浸边界的思路？

马乔

#include "fvCFD.H"
#include "pisoControl.H"

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

int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Transient solver for incompressible, laminar flow"
        " of Newtonian fluids."
    );

    #include "postProcess.H"

    #include "addCheckCaseOptions.H"
    #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();

    if(mesh.cellZones().size())
    {
	const cellZone& 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)
        );
	
	//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(mesh.solver(p.select(piso.finalInnerIter())));

                if (piso.finalNonOrthogonalIter())
                {
                    phi = phiHbyA - pEqn.flux();
                }
            }

            #include "continuityErrs.H"

            U = HbyA - rAU*fvc::grad(p);
            U.correctBoundaryConditions();
        }

        runTime.write();

        runTime.printExecutionTime(Info);
    }

    Info<< "End\n" << endl;

    return 0;
}

toposetDict
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

actions
(
    // porousBlockage
    {
        name    porousBlockageCellSet;
        type    cellSet;
        action  new;
        source  boxToCell;
        box     (0.04 0.04 -1) (0.06 0.06 1);
    }
    {
        name    porousBlockage;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        set     porousBlockageCellSet;
    }
);

Samuel-Tu

@cccrrryyy 对，现在我用的就是浸入边界法

Samuel-Tu

@马乔这是建立了一个集合，然后对这个集合赋予不同的粘度吗

马乔

modify the matrix, set the value of p and U as 0.

cccrrryyy

@Samuel-Tu 沉浸边界法流体域和固体域都计算吧？我记得它的核心思想就是把固体域也囊括到流体域的计算里面，通过修改计算域的控制方程（类似于NS方程的一个东西？），是的计算域内某一部分的表征是固体某一部分的表征是流体。如果你是想让计算域中其中一部分不参与计算，那和普通的只计算流体域不计算固体域（固体域作为流场边界）的区别是什么呢？挺好奇的。

Samuel-Tu

@cccrrryyy 浸入边界法其中的一种方法是加入附加体积力，确实如你所说是所有计算域都会计算，只是原本的固体域由于附加体积力的影响形成旋涡，不与流体域发生质量交换了，相当于就是固体了。

cccrrryyy

所以说你这个做的应该不是这个体积力的东西，是另外一种了。有什么参考文献可以分享一下么？我很好奇别的处理是什么样子的

Samuel-Tu

@cccrrryyy 我现在做的就是体积力的方法，我在想有没有其他办法。。

我是河滩

@马乔在 OpenFOAM有方法能够使一部分网格不参与计算吗？中说：

set the value of p and U as 0.

可以给出一些有关矩阵修改的参考文献吗？目前也需要对一些给定单元赋值（通过周围正常求解网格单元的值进行插值），OpenFoam中的fixedinternalValue边界条件也有类似功能，但不是很理解。

random_ran

(Sorry, I don't have chinese input right now)

I once had the exactly same issue!

Check my thesis Chapter 6.3 :

https://scholar.uwindsor.ca/cgi/viewcontent.cgi?article=8585&context=etd

I wrote a paper in this: https://ascelibrary.org/doi/pdf/10.1061/9780784415153.ch06

Figure 2 is the original plate. (Sorry due to copy-right reason, I cannot put it here.)

In my thesis, I had to open the plate with some holes in ICEM. So, I maual dig out the block out of the computation domain.

It's really a painful process.

But it works!

This is the original plate:

This is after the "dig":

东岳

        fvVectorMatrix UEqn
        (   
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

        // 固定值
        labelList cells = mesh.cellZones()["null"];
        vectorField value(cells.size(), Zero);
        UEqn.setValues(cells, value);

        if (piso.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));
        }

可以通过设定一个null区域，然后对区域的网格点强行赋值，但就像我之前说的，总觉得边界处需要处理。也就像 @cccrrryyy 说的，变成了浸没边界法。这个代码只是大体上是个思路，需要更深的研究。可以看看fvOptions里面一些源项的处理，都是操纵矩阵的，比如fixedTempresure，meanVelocityForce等，

Samuel-Tu

@东岳嗯嗯，我试一下看看能不能行。我记得setVAlue是切断了相邻面上的联系。我处理下剪切力试试

Samuel-Tu

@random_ran 主要的，我的边界随时要改变。这样的话就会一直改网格，有点像自适应网格了。。