OpenFOAM的一個簡單並行計算(計算房間內的流量通量)
來自專欄技術鄰CAE學院
11 人贊了文章
最近在利用OpenFOAM計算室內對流問題，其中通風量的大小是一個重要物理量，如果採用LES來計算，沒法利用後處理軟體得到每個時間步瞬態的值，我查找了OpenFOAM的tutorials，似乎也沒有計算某個打開的窗戶上通風量的case。因此編寫了一個簡單的並行計算通風量的程序，如果發現裡面有問題，歡迎大家指正，也歡迎大家隨意轉載。
計算通風量過程主要碰到的問題就是並行問題，為了演示，畫了3d建築物如下所示，前後兩個面上有矩形方框，示意為窗戶，在OpenFOAM中可以採用faceZone來標記這些面.(我是用topoSet生成的面，不知道各位有沒有什麼好的方法來取個faceZone).OpenFOAM進行分解時候會把這個faceSet分解到不同的processor文件夾中，也就是流量需要在不同的核心上進行計算，然後合在一起。
Info<< "Initializing the patch zone for calculating flow rate" << endl; IOdictionary flowRatePatchDict ( IOobject ( "flowRatePatchDict", "system", mesh, IOobject::MUST_READ, IOobject::NO_WRITE ) ); word faceZoneName_; flowRatePatchDict.lookup("faceZone") >> faceZoneName_; const vector windowDr(flowRatePatchDict.lookup("windowDirection")); label zoneId = mesh.faceZones().findZoneID(faceZoneName_); if (zoneId < 0) { FatalErrorInFunction << " Unknown face zone name: " << faceZoneName_ << ". Valid face zones are: " << mesh.faceZones().names() << nl << exit(FatalError); } const faceZone& fZone = mesh.faceZones()[zoneId];
上面這段代碼是faceZone的讀取與在網格中的識別，就不多敘述了，這裡主要介紹並行計算時候的策略。
以計算整個faceZone上面的總面積為例:
- 首先要獲得這個faceZone上面的個數，來生成一個總的數組儲存每個面的面積，建立一個數組faceFieldSizePerProc來存儲每個核心上分配得到的在faceZone上的面的個數，這個數組的大小就是核心數，可以調用Pstream::nProcs()來得到，再設置每個核心對應faceZone上的面的個數，Pstream::myProcNo()指的是當前核心。最後通過reduce函數進行合併數組。這樣就得到faceZone上總的面個數
labelList faceFieldSizePerProc(Pstream::nProcs(), 0.0); faceFieldSizePerProc[Pstream::myProcNo()] = fZone.size(); reduce( faceFieldSizePerProc, sumOp<labelList>() ); int totalSize = sum( faceFieldSizePerProc); Info << faceZoneName_ << "s size is " << totalSize << nl;
2. 下面就是計算faceZone上面的個數之和。
方法有點類似步驟一，只是每個核心上對應的faceZone上所有面的面積要做一個循環，最後再一組合就可以得到。
scalar totalArea(0); int n(0); scalar flowRateNew,flowRateMean(0); scalarList areaPerProc(Pstream::nProcs(), 0.0); forAll(fZone, i) { areaPerProc[Pstream::myProcNo()] += mesh.magSf()[fZone[i]] ; } reduce(areaPerProc, sumOp<scalarList>()); totalArea = sum( areaPerProc);
3. 最終就是在每個計算步的循環過程中，計算流量，其中面上的流量我用的是自己生成，而不是phi,主要因為topoSet產生的faceZone在分解之後，每個面的法向不一樣，最終有正有負，需要對面的法向進行統一。
最後放出全部的代碼
/*---------------------------------------------------------------------------* ========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox \ / O peration | \ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \/ M anipulation |-------------------------------------------------------------------------------License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.Application pimpleFoamDescription Large time-step transient solver for incompressible, turbulent flow, using the PIMPLE (merged PISO-SIMPLE) algorithm. Sub-models include: - turbulence modelling, i.e. laminar, RAS or LES - run-time selectable MRF and finite volume options, e.g. explicit porosity*---------------------------------------------------------------------------*/#include "fvCFD.H"#include "singlePhaseTransportModel.H"#include "turbulentTransportModel.H"#include "pimpleControl.H"#include "fvOptions.H"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //int main(int argc, char *argv[]){ #include "postProcess.H" #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" #include "createControl.H" #include "createTimeControls.H" #include "createFields.H" #include "createFvOptions.H" #include "initContinuityErrs.H" turbulence->validate(); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Initialize the patch zone for calculating flow rate // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "Initializing the patch zone for calculating flow rate" << endl; IOdictionary flowRatePatchDict ( IOobject ( "flowRatePatchDict", "system", mesh, IOobject::MUST_READ, IOobject::NO_WRITE ) ); word faceZoneName_; flowRatePatchDict.lookup("faceZone") >> faceZoneName_; const vector windowDr(flowRatePatchDict.lookup("windowDirection")); label zoneId = mesh.faceZones().findZoneID(faceZoneName_); if (zoneId < 0) { FatalErrorInFunction << " Unknown face zone name: " << faceZoneName_ << ". Valid face zones are: " << mesh.faceZones().names() << nl << exit(FatalError); } const faceZone& fZone = mesh.faceZones()[zoneId]; labelList faceFieldSizePerProc(Pstream::nProcs(), 0.0); faceFieldSizePerProc[Pstream::myProcNo()] = fZone.size(); reduce( faceFieldSizePerProc, sumOp<labelList>() ); int totalSize = sum( faceFieldSizePerProc); Info << faceZoneName_ << "s size is " << totalSize << nl; label labelListBefore(0); for(int i=0; i < Pstream::myProcNo(); i++) { labelListBefore += faceFieldSizePerProc[i]; } scalar totalArea(0); int n(0); scalar flowRateNew,flowRateMean(0); scalarList areaPerProc(Pstream::nProcs(), 0.0); forAll(fZone, i) { areaPerProc[Pstream::myProcNo()] += mesh.magSf()[fZone[i]] ; } reduce(areaPerProc, sumOp<scalarList>()); totalArea = sum( areaPerProc); Info << "Total area is " << totalArea << nl; std::fstream fout; std::fstream ftest; vector wlT(0,0,1); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "
Starting time loop
" << endl; while (runTime.run()) { #include "readTimeControls.H" #include "CourantNo.H" #include "setDeltaT.H" runTime++; flowRateNew = 0; Info<< "Time = " << runTime.timeName() << nl << endl; // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) { #include "UEqn.H" // --- Pressure corrector loop while (pimple.correct()) { #include "pEqn.H" } if (pimple.turbCorr()) { laminarTransport.correct(); turbulence->correct(); } } runTime.write(); surfaceVectorField u_interpolate ( IOobject ( "u_interpolate", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), fvc::interpolate(U) ); scalarList flowRateList(totalSize, 0.0); if( faceFieldSizePerProc[Pstream::myProcNo()] > 0 ) { for( int i = 0; i < faceFieldSizePerProc[Pstream::myProcNo()]; i++) { if(scalar(mesh.Sf()[fZone[i]]&windowDr) > 0) { flowRateList[ labelListBefore + i ] = mesh.Sf()[fZone[i]] & u_interpolate[fZone[i]]; } else { flowRateList[ labelListBefore + i ] = -mesh.Sf()[fZone[i]] & u_interpolate[fZone[i]]; } } } reduce( flowRateList, sumOp<scalarList>() ); flowRateNew = sum(flowRateList) ; ftest.open("test.dat",std::ios::trunc|std::ios::out); forAll(flowRateList, i) { ftest << flowRateList[i] << std::endl; } ftest.close(); Info << "Instantaneous flow rate is " << flowRateNew << nl; n++; flowRateMean = ((n-1)*flowRateMean+flowRateNew)/n; Info << "Mean flow rate is " << flowRateMean << nl; fout.open("flowRate.dat",std::ios::app); fout << runTime.timeName() << " " << flowRateNew << " " << flowRateMean << "
"; fout.close(); /**/ Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End
" << endl; return 0;}
另外友情提醒:在分解時候記得perserve那些faceZone，不然會產生重複的面的現象。
推薦閱讀：