一個標(biāo)準(zhǔn)的PISO循環(huán)解析

CFD模型理論對于不可壓流體的控制方程主要包括連續(xù)性方程及動量守恒方程，對于CFD-DEM耦合模型中的液相，其連續(xù)性方程及動量守恒方程如下：連續(xù)性方程： (STYLEREF1\s4-SEQ公式\*ARABIC\s144)動量守恒方程： (STYLEREF1\s4-SEQ公式\*ARABIC\s145)在CFD-DEM模型中我們采用下標(biāo)f來表示連續(xù)相，下標(biāo)p表示離散相。根據(jù)Boussinesq假設(shè)，應(yīng)力項τf包括基于分子粘性的平均應(yīng)力τm及湍流應(yīng)力τt兩部分，不可壓縮流體消掉密度，其計算公式變換為如下所示：基于分子粘性的平均應(yīng)力： (STYLEREF1\s4-SEQ公式\*ARABIC\s146)將雷諾應(yīng)力分解為兩個部分： 偏分量： (STYLEREF1\s4-SEQ公式\*ARABIC\s147) 各向同性分量： (STYLEREF1\s4-SEQ公式\*ARABIC\s148)將式4-48歸入壓力項，令壓力： (STYLEREF1\s4-SEQ公式\*ARABIC\s149)則不可壓縮流動的有效應(yīng)力表示為： (STYLEREF1\s4-SEQ公式\*ARABIC\s150) 其中有效粘度定義： (STYLEREF1\s4-SEQ公式\*ARABIC\s151)動量方程中相間動量交換定義如下： (STYLEREF1\s4-SEQ公式\*ARABIC\s152)相間動量交換液相部分做隱式處理，動量方程根據(jù)顯隱式劃分變換為： (STYLEREF1\s4-SEQ公式\*ARABIC\s153)其中右側(cè)部分為隱式求解部分，顯示部分為顯示求解部分。對于液相流場的計算采用PISO算法，速度預(yù)測時動量方程左側(cè)可以Cuf*： (STYLEREF1\s4-SEQ公式\*ARABIC\s154)對于速度預(yù)測方程左側(cè)代碼如下所示：fvVectorMatrixUEqnfvVectorMatrixUEqn(fvm::ddt(voidfraction,U)+fvm::div(phi,U)//+turbulence->divDevReff(U)+particleCloud.divVoidfractionTau(U,voidfraction)==-fvm::Sp(Ksl/rho,U));對于方程組Ax=b，源碼中fvm表示用于隱式操作，生成方程左側(cè)的系數(shù)矩陣，與此相對的fvc則用于顯式運算，用于方程的右側(cè)矩陣。代碼中的particleCloud.divVoidfractionTau(U,voidfraction)用來計算應(yīng)力部分，該代碼的解析在cfdemCloud類中實現(xiàn)，具體代碼如下所示：tmp<fvVectorMatrix>cfdemCloud::divVoidfractionTau(volVectorField&U,volScalarField&voidfraction)consttmp<fvVectorMatrix>cfdemCloud::divVoidfractionTau(volVectorField&U,volScalarField&voidfraction)const{return(-fvm::laplacian(voidfractionNuEff(voidfraction),U)-fvc::div(voidfractionNuEff(voidfraction)*dev(fvc::grad(U)().T())));}解析為張量運算形式： (STYLEREF1\s4-SEQ公式\*ARABIC\s155) 其中： (STYLEREF1\s4-SEQ公式\*ARABIC\s156)速度預(yù)測的右側(cè)方程用一下代碼求解：solvesolve(UEqn==fvc::reconstruct(( fvc::interpolate(voidfraction)*(g&mesh.Sf())-fvc::interpolate(voidfraction)*fvc::snGrad(p)*mesh.magSf()+(fvc::interpolate(Ksl/rho*Us)&mesh.Sf()))));對于方程右側(cè)顯式部分的處理并沒有直接使用體積場，例如對于壓力的梯度并沒有使用fvc::grad(p)，而是應(yīng)用了fvc::snGrad(p)*mesh.magSf()實現(xiàn)，然后通過fvc::reconstruct命令將面通量場創(chuàng)建體積場，由于面通量具有連續(xù)性，因此該法實現(xiàn)基于擬交錯網(wǎng)格的數(shù)值求解，具有更高的插值效率，并且不存在同位網(wǎng)格標(biāo)量棋盤狀分布時的問題。矩陣C為系數(shù)矩陣，uf*為速度預(yù)測值，其中矩陣C可以分解為對角線矩陣A及剩余部分矩陣H’，動量方程則可以進一步寫作： (STYLEREF1\s4-SEQ公式\*ARABIC\s157)PISO算法中對速度預(yù)測后需要進行速度校正，校正后的速度uf**是由預(yù)測速度uf*及校正壓力p*計算得來，對式（4-55）代入校正速度及校正壓力： (STYLEREF1\s4-SEQ公式\*ARABIC\s158)由于矩陣A為對角矩陣，因此可以方便地求其逆矩陣，另外引入H’uf**=-H，對速度校正方程則可以寫為： (STYLEREF1\s4-SEQ公式\*ARABIC\s159)校正速度的目的是為了使校正后的速度滿足連續(xù)性方程，因此校正后的速度uf**應(yīng)當(dāng)滿足液相的連續(xù)性方程： (STYLEREF1\s4-SEQ公式\*ARABIC\s160)將校正方程式（4-57）代入連續(xù)性方程式（4-58），我們可以得到： (STYLEREF1\s4-SEQ公式\*ARABIC\s161)進一步重寫為泊松方程： (STYLEREF1\s4-SEQ公式\*ARABIC\s162)由式（4-60）計算得到校正壓力p*就可以計算校正速度uf**，這樣就可以根據(jù)壓力方程及速度校正方程進行迭代，直至滿足殘差條件。速度校正的源碼如下所示：for(intcorr=0;corr<nCorrSoph;corr++)for(intcorr=0;corr<nCorrSoph;corr++){volScalarFieldrUA=1.0/UEqn.A();surfaceScalarFieldrUAvoidfractionf("(1|A(U))",fvc::interpolate(rUA*voidfraction));U=rUA*UEqn.H();phi=fvc::interpolate(U*voidfraction)&mesh.Sf();//+fvc::ddtPhiCorr(rUA,U,phi)surfaceScalarFieldphiS(fvc::interpolate(Us)&mesh.Sf());surfaceScalarFieldphiGes=phi+rUAvoidfractionf*(fvc::interpolate(Ksl/rho)*phiS+fvc::interpolate(voidfraction)*(g&mesh.Sf()));volScalarFieldrUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);//Non-orthogonalpressurecorrectorloopfor(intnonOrth=0;nonOrth<=nNonOrthCorr;nonOrth++){//PressurecorrectorfvScalarMatrixpEqn(fvm::laplacian(rUAvoidfraction,p)==fvc::div(phiGes)+fvc::ddt(voidfraction));pEqn.setReference(pRefCell,pRefValue);if(corr==nCorr-1&&nonOrth==nNonOrthCorr){pEqn.solve(mesh.solver("pFinal"));}else{pEqn.solve();}if(nonOrth==nNonOrthCorr){phiGes-=pEqn.flux();}}//endnon-orthogonalcorrectorloop#include"continuityErrs.H"U+=rUA*fvc::reconstruct((phiGes-phi)/rUAvoidfractionf);U.correctBoundaryConditions();}//endpisoloop全部的源碼（包括modeltype=B，壓力場共享）//Pressure-velocityPISOcorrector{//MomentumpredictorfvVectorMatrixUEqn(fvm::ddt(voidfraction,U)+fvm::div(phi,U)//+turbulence->divDevReff(U)+particleCloud.divVoidfractionTau(U,voidfraction)==-fvm::Sp(Ksl/rho,U));UEqn.relax();if(momentumPredictor){//solveUEqnif(modelType=="B"){ solve ( UEqn == fvc::reconstruct ( ( fvc::interpolate(voidfraction)*(g&mesh.Sf()) -fvc::snGrad(p)*mesh.magSf() +(fvc::interpolate(Ksl/rho*Us)&mesh.Sf()) ) ) ); }else { solve ( UEqn == fvc::reconstruct ( ( fvc::interpolate(voidfraction)*(g&mesh.Sf()) -fvc::interpolate(voidfraction)*fvc::snGrad(p)*mesh.magSf() +(fvc::interpolate(Ksl/rho*Us)&mesh.Sf()) ) ) ); }}//PISOloop//for(intcorr=0;corr<nCorr;corr++)intnCorrSoph=nCorr+5*pow((1-particleCloud.dataExchangeM().timeStepFraction()),1);for(intcorr=0;corr<nCorrSoph;corr++){volScalarFieldrUA=1.0/UEqn.A();surfaceScalarFieldrUAvoidfractionf("(1|A(U))",fvc::interpolate(rUA*(voidfraction+scalar(0.0001))));U=rUA*UEqn.H();phi=fvc::interpolate(U*voidfraction)&mesh.Sf();//+fvc::ddtPhiCorr(rUA,U,phi)surfaceScalarFieldphiS(fvc::interpolate(Us)&mesh.Sf());surfaceScalarFieldphiGes=phi+rUAvoidfractionf*(fvc::interpolate(Ksl/rho)*phiS+fvc::interpolate(voidfraction)*(g&mesh.Sf()));//surfaceScalarFieldphiGes=phi+rUAvoidfractionf*(fvc::interpolate(Ksl/rho)*phiS*;volScalarFieldrUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);if(modelType=="A")rUAvoidfraction=volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);//Non-orthogonalpressurecorrectorloopfor(intnonOrth=0;nonOrth<=nNonOrthCorr;nonOrth++){//PressurecorrectorfvScalarMatrixpEqn(fvm::laplacian(rUAvoidfraction,p)==fvc::div(phiGes)+fvc::ddt(voidfraction));pEqn.se