如何在多孔介质模型中使用自定义阻力分布?



  • 最近根据OpenFOAM中DarcyForchheimer模型想要实现自定义的阻力分布曲线,但是写好之后在porousSimpleFoam中运行会报错,查看流场发现阻力明显过大了,不知道能不能有简单点的定义方法呢?

    这是运行的报错信息

    Time = 166.5
    
    GAMG:  Solving for p, Initial residual = 0.0725239, Final residual = 4.0494e-05, No Iterations 6
    time step continuity errors : sum local = 7.31179e-06, global = -6.23766e-08, cumulative = 2.94484e-05
    smoothSolver:  Solving for epsilon, Initial residual = 0.038307, Final residual = 1.89737e-05, No Iterations 5
    smoothSolver:  Solving for k, Initial residual = 0.0462341, Final residual = 3.10477e-05, No Iterations 5
    ExecutionTime = 2.25 s  ClockTime = 3 s
    
    Time = 167
    
    #0  Foam::error::printStack(Foam::Ostream&) at ??:?
    #1  Foam::sigFpe::sigHandler(int) at ??:?
    #2  ? at ??:?
    #3  Foam::GAMGSolver::scale(Foam::Field<double>&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field<double> const&, unsigned char) const at ??:?
    #4  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
    #5  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
    #6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
    #7  Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/veen/Disk/OpenFOAM/OpenFOAM-8/platforms/linux64GccDPInt32Opt/bin/porousSimpleFoam"
    #8  Foam::fvMatrix<double>::solve() in "/home/veen/Disk/OpenFOAM/OpenFOAM-8/platforms/linux64GccDPInt32Opt/bin/porousSimpleFoam"
    #9  ? in "/home/veen/Disk/OpenFOAM/OpenFOAM-8/platforms/linux64GccDPInt32Opt/bin/porousSimpleFoam"
    #10  __libc_start_main at ??:?
    #11  ? in "/home/veen/Disk/OpenFOAM/OpenFOAM-8/platforms/linux64GccDPInt32Opt/bin/porousSimpleFoam"
    [1]    14023 floating point exception (core dumped)  porousSimpleFoam
    
    

    有可能是我的代码写的有点问题,在调用网格坐标的时候使用mesh.C()提示没有声明,但是已经在前边有写const fvMesh& mesh 这是在model中定义的计算方式

    	forAll(cellZoneIDs_, zoneI)
    		{
    			const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
    
    			D_[zoneI].setSize(cells.size());
    			F_[zoneI].setSize(cells.size());
    
    			scalar a1 = geo_[0]*(1-betaXYZ_.value().x());
    			scalar b1 = geo_[1]*(1-betaXYZ_.value().y());
    			scalar c1 = geo_[2]*(1-betaXYZ_.value().z());
    			scalar c2 = geo_[3]*(1-betaXYZ_.value().z());
    			scalar d_ = geo_[3]+geo_[4];
    			scalar pox = pow(poXYZ_.value().x(), alphaXYZ_.value().x());
    			scalar poy = pow(poXYZ_.value().y(), alphaXYZ_.value().y());
    			scalar poz = pow(poXYZ_.value().z(), alphaXYZ_.value().z());
    			scalar gammax;
    			scalar gammay;
    			scalar gammaz;
    	 
    			forAll(cells, i)
    			{
    			    scalar x = mesh_.C()[i].x();
    			    scalar y = mesh_.C()[i].y();
    			    scalar z = mesh_.C()[i].z();
    			    
    			    scalar dis1 = pow(x,2)/pow(a1,2)+pow(y,2)/pow(b1,2)+pow(z-d_,2)/pow(c1,2);
    			    scalar dis2 = pow(x,2)/pow(a1,2)+pow(y,2)/pow(b1,2)+pow(z-d_,2)/pow(c2,2);
    			    scalar dis1max = max(pow(geo_[0],2)/pow(a1,2), max(pow(geo_[1],2)/pow(b1,2), pow(geo_[2]-d_,2)/pow(c1,2)));
    			    scalar dis2max = max(pow(geo_[0],2)/pow(a1,2), max(pow(geo_[1],2)/pow(b1,2), pow(geo_[2]-d_,2)/pow(c2,2)));
    			    
    			    if (z>geo_[3] && dis1<1)
    			    {
    			        //scalar dis = 1-dis1;
    			        gammax = pox;
    			        gammay = poy;
    			        gammaz = poz;
    			    }
    			    else if (z<geo_[3] && dis2<1)
    			    {
    			        //scalar dis = 1-dis2;
    			        gammax = pox;
    			        gammay = poy;
    			        gammaz = poz;
    			    }
    			    else if (z>=geo_[3] && dis1>=1)
    			    {
    			        scalar dis = dis1-1;
    			        scalar dism = dis1max-1;
    			        gammax = (1-dis/dism)*pox;
    			        gammay = (1-dis/dism)*poy;
    			        gammaz = (1-dis/dism)*poz;
    			    }
    			    else
    			    {
    			        scalar dis = dis2-1;
    			        scalar dism = dis2max-1;
    			        gammax = (1-dis/dism)*pox;
    			        gammay = (1-dis/dism)*poy;
    			        gammaz = (1-dis/dism)*poz;
    			    }
    			    
    			    D_[zoneI][i] = Zero;
    			    D_[zoneI][i].xx() = gammax*dXYZ_.value().x();
    			    D_[zoneI][i].yy() = gammay*dXYZ_.value().y();
    			    D_[zoneI][i].zz() = gammaz*dXYZ_.value().z();
    			    
    			    D_[zoneI][i] = coordSys_.R().transformTensor(D_[zoneI][i]);
    
    			    // leading 0.5 is from 1/2*rho
    			    F_[zoneI][i] = Zero;
    			    F_[zoneI][i].xx() = 0.5*pow(gammax,2)*fXYZ_.value().x();
    			    F_[zoneI][i].yy() = 0.5*pow(gammay,2)*fXYZ_.value().y();
    			    F_[zoneI][i].zz() = 0.5*pow(gammaz,2)*fXYZ_.value().z();
    			    
    			    F_[zoneI][i] = coordSys_.R().transformTensor(F_[zoneI][i]);
    			}
    		}
    


  • @veen 破案了,求解器发散是因为我把坐标写错了,F_[zoneI]的值不对造成了流场阻力过大,主要思路是没有问题的。



  • 最近几天又仔细算了一下,发现算出来阻力还是不对,本来定义的是个阻力项,但是实际算出来是加速作用。
    想请教一下各位老师有没有什么解决的办法?


Log in to reply
 


CFD中文网 | 东岳流体学术 | 东岳流体商业 | 吉ICP备20003622号-1