CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

    sprayFoam计算等压环境下点火延迟时间

    OpenFOAM
    4
    8
    288
    正在加载更多帖子
    • 从旧到新
    • 从新到旧
    • 最多赞同
    回复
    • 在新帖中回复
    登录后回复
    此主题已被删除。只有拥有主题管理权限的用户可以查看。
    • 香柏树
      香柏树 最后由 编辑

      目前在尝试使用sprayFoam等压环境下的C12H26点火延迟时间,由于sprayFoam中通过泊松方程求解压力(见原始pEqn.H),网格内监测压力显示不能保持等压状态。因此参考chemFoam的源代码,想把泊松方程改为完全气体状态方程P=rhoRT(见修改的pEqn.H),同时在sprayFoam.C中注释掉压力修正的过程(见修改的sprayFoam.C)。另外能量方程也改为求解Enthalpy的形式。在OpenFOAM6环境下编译没问题,但是计算两个时间步后即浮点溢出了,计算设置见下楼的图1,错误提示下楼的图2。请问各位前辈有没有解决方法啊?
      原始pEqn.H

      rho = thermo.rho();
      rho = max(rho, rhoMin);
      rho = min(rho, rhoMax);
      rho.relax();
      
      volScalarField rAU(1.0/UEqn.A());
      surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
      volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
      
      if (pimple.nCorrPISO() <= 1)
      {
          tUEqn.clear();
      }
      
      if (pimple.transonic())
      {
          surfaceScalarField phid
          (
              "phid",
              fvc::interpolate(psi)
             *(
                  fvc::flux(HbyA)
                + MRF.zeroFilter
                  (
                      rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
                  )
              )
          );
      
          MRF.makeRelative(fvc::interpolate(psi), phid);
      
          while (pimple.correctNonOrthogonal())
          {
              fvScalarMatrix pEqn
              (
                  fvm::ddt(psi, p)
                + fvm::div(phid, p)
                - fvm::laplacian(rhorAUf, p)
               ==
                  parcels.Srho()
                + fvOptions(psi, p, rho.name())
              );
      
              pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
      
              if (pimple.finalNonOrthogonalIter())
              {
                  phi == pEqn.flux();
              }
          }
      }
      else
      {
          surfaceScalarField phiHbyA
          (
              "phiHbyA",
              (
                  fvc::flux(rho*HbyA)
                + rhorAUf*fvc::ddtCorr(rho, U, phi)
              )
          );
      
          MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
      
          // Update the pressure BCs to ensure flux consistency
          constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
      
          while (pimple.correctNonOrthogonal())
          {
              fvScalarMatrix pEqn
              (
                  fvm::ddt(psi, p)
                + fvc::div(phiHbyA)
                - fvm::laplacian(rhorAUf, p)
               ==
                  parcels.Srho()
                + fvOptions(psi, p, rho.name())
              );
      
              pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
      
              if (pimple.finalNonOrthogonalIter())
              {
                  phi = phiHbyA + pEqn.flux();
              }
          }
      }
      
      #include "rhoEqn.H"
      #include "compressibleContinuityErrs.H"
      
      // Explicitly relax pressure for momentum corrector
      p.relax();
      
      // Recalculate density from the relaxed pressure
      rho = thermo.rho();
      rho = max(rho, rhoMin);
      rho = min(rho, rhoMax);
      rho.relax();
      Info<< "rho max/min : " << max(rho).value()
          << " " << min(rho).value() << endl;
      
      U = HbyA - rAU*fvc::grad(p);
      U.correctBoundaryConditions();
      fvOptions.correct(U);
      K = 0.5*magSqr(U);
      
      if (thermo.dpdt())
      {
          dpdt = fvc::ddt(p);
      }
      
      

      修改的pEqn.H

      {
          rho = thermo.rho();
          p = rho*Rspecific*thermo.T();  
      }
      

      修改的sprayFoam.C

      /*---------------------------------------------------------------------------*\
        =========                 |
        \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
         \\    /   O peration     | Website:  https://openfoam.org
          \\  /    A nd           | Copyright (C) 2011-2018 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
          sprayFoam
      
      Description
          Transient solver for compressible, turbulent flow with a spray particle
          cloud.
      
      \*---------------------------------------------------------------------------*/
      
      #include "fvCFD.H"
      #include "turbulentFluidThermoModel.H"
      #include "basicSprayCloud.H"
      #include "psiReactionThermo.H"
      #include "CombustionModel.H"
      #include "radiationModel.H"
      #include "SLGThermo.H"
      #include "pimpleControl.H"
      #include "fvOptions.H"
      
      // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
      
      int main(int argc, char *argv[])
      {
          #include "postProcess.H"
      
          #include "setRootCaseLists.H"
          #include "createTime.H"
          #include "createMesh.H"
          #include "createControl.H"
          #include "createTimeControls.H"
          #include "createFields.H"
          #include "createFieldRefs.H"
          #include "compressibleCourantNo.H"
          #include "setInitialDeltaT.H"
          
          turbulence->validate();
      
          // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
      
          Info<< "\nStarting time loop\n" << endl;
      
          while (runTime.run())
          {
              #include "readTimeControls.H"
              #include "compressibleCourantNo.H"
              #include "setDeltaT.H"
      
              runTime++;
      
              Info<< "Time = " << runTime.timeName() << nl << endl;
      
              parcels.evolve();
      
              #include "rhoEqn.H"
              #include "YEqn.H"
              #include "EEqn.H"
              #include "pEqn.H"
      
              /*
              if (!pimple.frozenFlow())
              {
                  #include "rhoEqn.H"
      
                  // --- Pressure-velocity PIMPLE corrector loop
                  while (pimple.loop())
                  {
                      // #include "UEqn.H"
                      #include "YEqn.H"
                      #include "EEqn.H"
      
                     
                      // --- Pressure corrector loop
                      while (pimple.correct())
                      {
                          #include "pEqn.H"
                      }
      
                      if (pimple.turbCorr())
                      {
                          turbulence->correct();
                      }
                  }
      
      
                  rho = thermo.rho();
      
                  if (runTime.write())
                  {
                      combustion->Qdot()().write();
                  }
              }
              else
              {
                  if (runTime.writeTime())
                  {
                      parcels.write();
                  }
              }
               */
      
              Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
                  << "  ClockTime = " << runTime.elapsedClockTime() << " s"
                  << nl << endl;
          }
      
          Info<< "End\n" << endl;
      
          return 0;
      }
      
      
      // ************************************************************************* //
      
      
      香柏树 1 条回复 最后回复 回复 引用
      • 香柏树
        香柏树 @香柏树 最后由 编辑

        68816195-0e44-46d3-9f95-df533596c80d-图片.png
        f834dba5-f097-4c55-ac86-66b125a60c38-图片.png

        1 条回复 最后回复 回复 引用
        • 李东岳
          李东岳 管理员 最后由 编辑

          @香柏树 在 sprayFoam计算等压环境下点火延迟时间 中说:

          因此参考chemFoam的源代码,想把泊松方程改为完全气体状态方程P=rhoRT(见修改的pEqn.H)

          你说的点火之类我不太熟悉。但针对你的这方面描述,chemFoam不涉及到整个场的连续性。sprayFoam要满足质量守恒。你这样做从算法上有点奇怪。

          后续的截图看起来发散了,除以0了,当然原因肯定是非常本质的。

          CFD高性能服务器 http://dyfluid.com/servers.html

          香柏树 1 条回复 最后回复 回复 引用
          • 香柏树
            香柏树 @李东岳 最后由 编辑

            @李东岳 谢谢李老师指导,这么修改确实不太正确。我从算法上看看如何能够合理些。

            星 2 条回复 最后回复 回复 引用
            • 星
              星星星星晴 @香柏树 最后由 编辑

              @香柏树 除以0可能是最大的问题,但是0在哪呢,请问你怎么处理parcel 极小的部分的,会不会是parcel部分,因为你点火了,就出现了这种情况,在clockTime以后就是计算Co 然后solving parcel了,也有可能是你parcel的maxCo=1, parcel在你计算的时候已经飞出了对应的cell 也有可能报错,尝试分析分析这俩地方,看看你的parcel都什么样,看看你的MaxCo减小到0.1或者0.01看看,还是turbulence爆了

              m.sui20@foxmail.com

              1 条回复 最后回复 回复 引用
              • 星
                星星星星晴 @香柏树 最后由 编辑

                @香柏树 另外你直接把整个PIMPLE循环给关掉了,是不是保留一下然后PISO模式运行就行呢, 瞬态

                m.sui20@foxmail.com

                香柏树 1 条回复 最后回复 回复 引用
                • 香柏树
                  香柏树 @星星星星晴 最后由 编辑

                  @星星星星晴 谢谢前辈!我尝试下。

                  1 条回复 最后回复 回复 引用
                  • 流
                    流体中二少年 最后由 编辑

                    前辈,你好!请问sprayFoam中的液滴蒸发率怎么获得?

                    1 条回复 最后回复 回复 引用
                    • First post
                      Last post