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



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



  • 把网格抠除?


  • OpenFOAM教授

    @Samuel-Tu 试过,可行。



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



  • @dzw05 :qinqin: 感谢



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



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



  • @dzw05 请问您具体是如何处理的?:xiexie:


  • OpenFOAM教授

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



  • @Samuel-TuOpenFOAM有方法能够使一部分网格不参与计算吗? 中说:

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

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



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



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



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


  • OpenFOAM副教授

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


  • 离散相副教授

    #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;
        }
    );
    

    UEqn.png pEqn.png
    :wocao:



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



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


  • 离散相副教授

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


  • OpenFOAM副教授

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



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


  • OpenFOAM副教授

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



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



  • @马乔OpenFOAM有方法能够使一部分网格不参与计算吗? 中说:

    set the value of p and U as 0.

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


  • OpenFOAM副教授

    (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:

    69e5a95f-8963-430f-aded-3ebc52e3aed5-image.png

    This is after the "dig":
    97acbf4d-9894-4c21-88e7-34f655ab3fd7-image.png



  •         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等,



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



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


Log in to reply
 

CFD中文网 2016 - 2020 | 京ICP备15017992号-2