如何在多孔介质模型中使用自定义阻力分布？
-
最近根据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]); } } -
我又回来了，这次发现虽然写了用坐标去计算不同点的阻力系数，但是实际计算中发现还是均一的，以为是定义的变量被重新赋值了，重写了所有的中间变量发现也没办法解决，目前认为是坐标计算的不对或者openfoam porosityModel不支持非均一的阻力系数设置
如果想用场自定义阻力需要如何操作呢？直接改求解器加入源项，然后在0文件夹写初始场？forAll(cells, i) { gamma_[zoneI][i] = Zero; scalar x = mesh_.C()[i].x(); scalar y = mesh_.C()[i].y(); scalar z = mesh_.C()[i].z(); scalar dir = y; if (xxx) { gamma_[zoneI][i].xx() = xxx; gamma_[zoneI][i].yy() = xxx; gamma_[zoneI][i].zz() = xxx; } D_[zoneI][i] = Zero; D_[zoneI][i].xx() += pow(gamma_[zoneI][i].xx(), 1/3) * dXYZ_.value().x(); D_[zoneI][i].yy() += pow(gamma_[zoneI][i].yy(), 1/3) * dXYZ_.value().y(); D_[zoneI][i].zz() += pow(gamma_[zoneI][i].zz(), 1/3) * 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 * fXYZ_.value().x(); F_[zoneI][i].yy() += 0.5 * fXYZ_.value().y(); F_[zoneI][i].zz() += 0.5 * fXYZ_.value().z(); F_[zoneI][i] = coordSys_.R().transformTensor(F_[zoneI][i]); }