OpenFOAM有方法能够使一部分网格不参与计算吗？
-
@Samuel-Tu 在 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; } );
感谢
