Euler-Lagrange的一些解析,sprayFoam
-
@星星星星晴 前辈,您好!我在使用sprayFoam喷射液滴的时候,在两个不同进气流量的case进行计算,case1为100%空气,case2为200%空气,所有的constant,system里的文件设置都一样,仅仅改变了进气流量。进行液滴后处理时发现,Case2的液滴数目(nParticles)比Case1少了近一半,请问您知道是什么原因么?液滴喷射模型设置具体如下:
model1 { type coneNozzleInjection; SOI 0.029; massTotal 0.0044166; parcelBasisType mass; injectionMethod disc; flowType constantVelocity; UMag 33.0; outerDiameter 0.0007; innerDiameter 0; duration 1.0; position (0 0 -0.03584797); direction (0 0 1); parcelsPerSecond 2000000; massFlowRate 0.0044166; flowRateProfile table ( (0 0.0044166) (1 0.0044166) (2 0.0044166) (3 0.0044166) (4 0.0044166) (5 0.0044166) (6 0.0044166) (7 0.0044166) (8 0.0044166) (9 0.0044166) (10 0.0044166) ); Cd constant 0.9; thetaInner constant 25; thetaOuter constant 30; sizeDistribution { type RosinRammler; RosinRammlerDistribution { minValue 1e-05; maxValue 6e-05; d 4e-05; n 3; } } }
-
@星星星星晴 谢谢前辈!100%和200%代表计算域进口空气流量的相对值,比如100%是1kg/s的话,200%就是2kg/s。流场内空气速度增大了两倍左右,没有开启碰撞模型,sprayCloudProperties文件内容具体如下:
\*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format binary; class dictionary; location "constant"; object SprayCloudProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solution { active true; coupled true; transient yes; cellValueSourceCorrection on; maxCo 0.3; sourceTerms { schemes { rho explicit 1; U explicit 1; Yi explicit 1; h explicit 1; radiation explicit 1; } } interpolationSchemes { rho cell; U cellPoint; thermo:mu cell; T cell; Cp cell; kappa cell; p cell; } integrationSchemes { U Euler; T analytical; } } constantProperties { T0 300; // place holders for rho0 and Cp0 // - reset from liquid properties using T0 rho0 744; Cp0 2090; Tvap 298; Tbp 489; constantVolume false; } subModels { particleForces { sphereDrag; } injectionModels { model1 { type coneNozzleInjection; SOI 0.029; massTotal 0.0044166; parcelBasisType mass; injectionMethod disc; flowType constantVelocity; UMag 33.0; outerDiameter 0.0007; innerDiameter 0; duration 1.0; position (0 0 -0.03584797); direction (0 0 1); parcelsPerSecond 2000000; massFlowRate 0.0044166; flowRateProfile table ( (0 0.0044166) (1 0.0044166) (2 0.0044166) (3 0.0044166) (4 0.0044166) (5 0.0044166) (6 0.0044166) (7 0.0044166) (8 0.0044166) (9 0.0044166) (10 0.0044166) ); Cd constant 0.9; thetaInner constant 25; thetaOuter constant 30; sizeDistribution { type RosinRammler; RosinRammlerDistribution { minValue 1e-05; maxValue 6e-05; d 4e-05; n 3; } } } dispersionModel none; patchInteractionModel standardWallInteraction; heatTransferModel RanzMarshall; compositionModel singlePhaseMixture; phaseChangeModel liquidEvaporationBoil; surfaceFilmModel none; atomizationModel none; breakupModel ReitzDiwakar; // ReitzKHRT; stochasticCollisionModel none; radiation off; standardWallInteractionCoeffs { type rebound; } RanzMarshallCoeffs { BirdCorrection true; } singlePhaseMixtureCoeffs { phases ( liquid { C12H26 1; } ); } liquidEvaporationBoilCoeffs { enthalpyTransfer enthalpyDifference; activeLiquids ( C12H26 ); } ReitzDiwakarCoeffs { solveOscillationEq yes; Cbag 6; Cb 0.785; Cstrip 0.5; Cs 10; } /* ReitzKHRTCoeffs { solveOscillationEq yes; B0 0.61; B1 40; Ctau 1; CRT 0.1; msLimit 0.2; WeberLimit 6; } */ TABCoeffs { y0 0; yDot0 0; Cmu 10; Comega 8; WeCrit 12; } } cloudFunctions { voidFraction { type voidFraction; } } }
-
-
首先你Euler那边速度发生了变化,自然Lagrange这边也会收到影响,我看你的parcel size 是10-60um, 然后parcel的初始速度就是33m/s,不知道你空气速度如何,自然drag的区别很大,然后parcel就可能很快的飞出去了。
-
在此之上,你有破碎模型,如果气液相对速度过大,自然会破碎,你要具体看一下你的破碎模型是怎么处理质量的。
-
其次,你的问题中“Case2的液滴数目(nParticles)比Case1少了近一半” 我没有get到,如果指的是particle per parcel的话,自然是受到其他模型的影响, 我看你这边有涉及到热,破碎等模型。这些模型也会影响nParticles,在使用这些模型的时候,你也是需要知道这些模型是怎么算的,算什么的,怎么处理mass/momentum balance的。
-
另外flowrateprofile constant就行,我看你这个flowrate不是随时间变化的,而且就喷1秒。
-
-
@星星星星晴 感谢前辈悉心指导!关于您说的问题,我今天查阅了一些公式,关于您提的四个问题,具体回复如下:
- 根据图1中液滴阻力的公式,在100%和200%air中的阻力分别为3.7e-6N和2.3e-5N,也就是液滴受到的阻力变大了,是不是说明液滴更倾向于保持目前的球形状态呢?
- 破碎模型用的是ReitzDiwakar model,具体公式如图中所示,计算得到100%和200%air中的rstable分别为18μm和1.2μm,也就是说后者的液滴粒径更小。这里您说的破碎模型如何处理质量是什么意思呢?是指公式中的密度么?
- “Case2中的液滴数目(nParticles)”说错了,应该是流域内的parcels的数目,也就是图2中的8370.计算结果显示200%air case中的parcels数目比100%air case少了近一半。
- Flowrateprofile constant已经改正。
再次感谢前辈指点!
-
-
spheredrag的前提就是假设所有的parcel 保持球形,所以根据阻力计算公式你可以自己推导一下,到底怎么回事。。严格来说drag 不一定是阻力吧,要看相对速度方向。
-
你需要查一下关于破碎模型的原始论文,破碎同样也涉及到计算相对速度。在破碎模型中到底是如何假设破碎的,ReitzDiwakar 我不太了解,但是比如TAB,要计算一个y,KHRT也要计算波动的,而且是否生成了child parcel,生成的值是什么,你要看code的。单纯问,没办法解释,看code看原始论文,看of是如何implement的
-
这个不难理解,如果你流场的速度增加了,那你喷雾的parcel也是收到流场影响的啊,自然可能会发生parcel被加速,如果被加速了的话,有些parcel离开的快也可以理解,你也没有具体描述你的case,所以我只能根据你当前的描述说了。
-
-
@星星星星晴 在 Euler-Lagrange的一些解析,sprayFoam 中说:
@chengan-wang 这也简单啊,同理,你不需要判断是否碰到面,而是做一个判断是否在你的区域内即可,比如你需要判断在两个面之间,z-dir
if(p.position(2)<1 & p.position(2)>0) { 输出 }
但是这种方法会出现什么问题,比如一个parcel在这个区域待了3个$\Delta t_L$, 那就输出3次,就是你的输出文件几何倍数的增长。。
您好,我按照您的提示,思路如下:
先是定义了三维空间的范围,以及坐标质量的数组。scalar xmin = -0.01; scalar xmax = 0.01; scalar ymin = 0; scalar ymax = 0.1; scalar zmin = -0.01; scalar zmax = 0.01; scalar dxyz = 5e-3; const int nx = floor((xmax - xmin)/dxyz); const int ny = floor((ymax - ymin)/dxyz); const int nz = floor((zmax - zmin)/dxyz); int nxc; int nyc; int nzc; scalar xc[nx]; scalar yc[ny]; scalar zc[nz]; scalar mass3D[nx][ny][nz] = {0};
我只想最后时刻输出结果,所以先用if条件判断是否要输出数据
if ((this->owner().db().time().value()) == 0.02) { for (int i=0; i<nx; i++) { xc[i] = i*dxyz+xmin; } for (int j=0; j<ny; j++) { yc[j] = j*dxyz+ymin; } for (int k=0; k<nz; k++) { zc[k] = k*dxyz+zmin; } nxc = floor((p.position().component(0) - xmin)/dxyz); nyc = floor((p.position().component(1) - ymin)/dxyz); nzc = floor((p.position().component(2) - zmin)/dxyz); mass3D[nxc][nyc][nzc] += p.nParticle()*p.mass(); std::ofstream outfile; outfile.setf(ios_base::fixed, ios_base::floatfield); outfile.precision(8); outfile.open("postProcessingWCA/Mass3D.txt", ios_base::app); for (int i=0; i<nx; i++) { for (int j=0; j<ny; j++) { for (int k=0; k<nz; k++) { outfile << xc[i] << tab << yc[j] << tab << zc[k] << mass3D[i][j][k] << tab << tab << nl; } } } }
之后我的想法就是根据粒子的位置计算索引位置 nxc nyc nzc,带入到mass3D[nxc][nyc][nzc] 相应的位置。
但是存在两个问题:- 每一个粒子统计之后都输出到了Mass3D.txt,数据量非常大,而我是想把所有的粒子质量统计到相应的空间点上累积;
- 另外就是定义方式scalar mass3D[nx][ny][nz] = {0}; 以及累积方法mass3D[nxc][nyc][nzc] += p.nParticle()*p.mass();感觉有问题,从结果上看感觉没有把相同位置的粒子质量加到一起。
想麻烦您帮忙看看,谢谢
-
- 数据量极大就是这个方法最大的问题,没办法。
- 你的mass3D没有输出,在结束这一步的迭代以后,在下一步迭代会被清空的。
你这一步已经很靠近方法1了,你不要用什么mass3D这个变量,你在solver中建立一个field, 比如叫summass,然后在cloudfunction的postmove 这个function中直接通过,就可以对一下field进行赋值了。
cellI = p.cell(); scalar& PPC = summass_->primitiveFieldRef() [cellI]; PPC += p.nParticle()*p.mass();
我的技能也是有限,所以没有办法给你解决所有问题,你可以自己试试,找找资料,如果发现更好的办法,也希望你分享出来。
对于方法2,我个人有个想法,你获得你所要区域的cell number, 然后仅仅对这些cell进行输出也可以的。我自己没做过,大家都是摸着石头过河,可能每个人的解决方法并不通用,很case- sensitive,这也是开源的一个弊端。。
祝你好运~
-
@星星星星晴 在 Euler-Lagrange的一些解析,sprayFoam 中说:
- 数据量极大就是这个方法最大的问题,没办法。
- 你的mass3D没有输出,在结束这一步的迭代以后,在下一步迭代会被清空的。
你这一步已经很靠近方法1了,你不要用什么mass3D这个变量,你在solver中建立一个field, 比如叫summass,然后在cloudfunction的postmove 这个function中直接通过,就可以对一下field进行赋值了。
cellI = p.cell(); scalar& PPC = summass_->primitiveFieldRef() [cellI]; PPC += p.nParticle()*p.mass();
我的技能也是有限,所以没有办法给你解决所有问题,你可以自己试试,找找资料,如果发现更好的办法,也希望你分享出来。
对于方法2,我个人有个想法,你获得你所要区域的cell number, 然后仅仅对这些cell进行输出也可以的。我自己没做过,大家都是摸着石头过河,可能每个人的解决方法并不通用,很case- sensitive,这也是开源的一个弊端。。
祝你好运~
谢谢。我在createFields.H建立了一个field:
volScalarField summass ( IOobject ( "summass", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh, dimensionedScalar ("zero", dimMass, 0.0) );
我没有理解“在cloudfunction的postmove 这个function中直接通过”是啥意思呢?麻烦了
-
@chengan-wang 前辈想问一下,GOFUN后来再没出particle simulation的第二讲了么,我看2019年本来有但是取消了
-
我刚看到,相信一年多过去你肯定已经解决了。
- 在0/ 文件夹中,添加lagrangian的文件夹,提供parcel所需的相关信息即可。setFields 应该也行,但是我的第一反应是这个方案。
- 随inlet边界?就是patchinjection, PatchFlowRateInjection 等喷射模型吧。
- 在spraycloudproperties中设置的C7H16是parcel的,在外面的/constant文件夹中是外面的气体。在paraview中,parcel是parcel,流场是流场。你可以轻松的看出来。在时间步中的lagrangian文件夹中也能看到每个物质的量的
-
@星星星星晴 您好,大佬,我最近潜心修心了一个月,了解了拉格朗日源文件里面一些大概的东西,我目前打算根据一个方程更新我的拉格朗日颗粒的直径,所以我根据一篇名叫“水质对梢涡空化初生的影响研究”的博士论文在修改拉格朗日的源文件库,最近遇到了很多问题,不知道您能否给我指引一点方向来学习。
这个代码原理就是在KinematicParcel.C中按照他原本的calcVelocity,自己定义一个气泡长大方程取名为calcDimater,但是这篇博士论文的代码似乎有点不全,我一直在修改。目前改的文件有以下这6个基础文件KinematicCloud.C KinematicCloud.H KinematicCloudI.H KinematicParcel.C KinematicParcel.H KinematicParcelI.H
我能编译成功源文件库,但是在编译求解器时候会有报错,于是我又简单更改了这三个文件夹
CollidingCloud.C CollidingCloud.H CollidingCloudI.H
目前我的库能够编译成功,就是存在warning。但是当我编译我的新的newDPMFoam时候总是有很多奇奇怪怪报错,目前我的报错就很奇怪。总而言之就是说我有一些参数我定义错误
g++ -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -O3 -DNoRepository -ftemplate-depth-100 -I./DPMTurbulenceModels/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/lagrangian/basic/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/lagrangian/intermediate/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/thermophysicalModels/specie/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/transportModels/compressible/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/thermophysicalModels/basic/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/thermophysicalModels/reactionThermo/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/thermophysicalModels/radiation/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/transportModels -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/transportModels/incompressible/singlePhaseTransportModel -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/TurbulenceModels/turbulenceModels/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/TurbulenceModels/incompressible/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/TurbulenceModels/phaseIncompressible/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/regionModels/regionModel/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/regionModels/surfaceFilmModels/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/finiteVolume/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/meshTools/lnInclude -IlnInclude -I. -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/OpenFOAM/lnInclude -I/home/zly/OpenFOAM/OpenFOAM-3.0.0/src/OSspecific/POSIX/lnInclude -fPIC -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64GccDPInt32Opt/ppDPMFoam.o -L/home/zly/OpenFOAM/OpenFOAM-3.0.0/platforms/linux64GccDPInt32Opt/lib \ -llagrangian -llagrangianIntermediate -llagrangianTurbulence -lthermophysicalFunctions -lspecie -lradiationModels -lincompressibleTransportModels -lturbulenceModels -lincompressibleTurbulenceModels -lDPMTurbulenceModels -lregionModels -lsurfaceFilmModels -lsampling -lfiniteVolume -lmeshTools -lOpenFOAM -ldl \ -lm -o /home/zly/OpenFOAM/zly-3.0.0/platforms/linux64GccDPInt32Opt/bin/ppDPMFoam Make/linux64GccDPInt32Opt/ppDPMFoam.o: In function `Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > >::setParcelThermoProperties(Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> >&, double) [clone .isra.656]': ppDPMFoam.C:(.text+0x1377): undefined reference to `Foam::KinematicParcel<Foam::particle>::Ro()' ppDPMFoam.C:(.text+0x13c9): undefined reference to `Foam::KinematicParcel<Foam::particle>::pgo()' Make/linux64GccDPInt32Opt/ppDPMFoam.o: In function `void Foam::KinematicParcel<Foam::particle>::setCellValues<Foam::KinematicParcel<Foam::particle>::TrackingData<Foam::CollidingCloud<Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > > > > >(Foam::KinematicParcel<Foam::particle>::TrackingData<Foam::CollidingCloud<Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > > > >&, double, int)': ppDPMFoam.C:(.text._ZN4Foam15KinematicParcelINS_8particleEE13setCellValuesINS2_12TrackingDataINS_14CollidingCloudINS_14KinematicCloudINS_5CloudINS_15CollidingParcelIS2_EEEEEEEEEEEEvRT_di[_ZN4Foam15KinematicParcelINS_8particleEE13setCellValuesINS2_12TrackingDataINS_14CollidingCloudINS_14KinematicCloudINS_5CloudINS_15CollidingParcelIS2_EEEEEEEEEEEEvRT_di]+0x23c): undefined reference to `Foam::KinematicParcel<Foam::particle>::TrackingData<Foam::CollidingCloud<Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > > > >::PInterp() const' collect2: error: ld returned 1 exit status make: *** [/home/zly/OpenFOAM/zly-3.0.0/platforms/linux64GccDPInt32Opt/bin/ppDPMFoam] Error 1
我这个原理上我感觉我的代码基本上都能说得通,但是报错经常就提醒我缺少一些东西,我现在就有一点迷茫,不知道怎么前进了,请大佬指点下迷津。我这个方法是否有可行性 附上一部分我的代码(//- myadd部分即为我添加部分)
KinematicParcel.C(这是方程的核心部分在里面,其他文件我只是按照原有的模板照猫画虎定义了一些参数,例如压力场的插值)
template<class ParcelType> template<class TrackData> void Foam::KinematicParcel<ParcelType>::setCellValues ( TrackData& td, const scalar dt, const label cellI ) { tetIndices tetIs = this->currentTetIndices(); rhoc_ = td.rhoInterp().interpolate(this->position(), tetIs); if (rhoc_ < td.cloud().constProps().rhoMin()) { if (debug) { WarningIn ( "void Foam::KinematicParcel<ParcelType>::setCellValues" "(" "TrackData&, " "const scalar, " "const label" ")" ) << "Limiting observed density in cell " << cellI << " to " << td.cloud().constProps().rhoMin() << nl << endl; } rhoc_ = td.cloud().constProps().rhoMin(); } Uc_ = td.UInterp().interpolate(this->position(), tetIs); muc_ = td.muInterp().interpolate(this->position(), tetIs); // Apply dispersion components to carrier phase velocity Uc_ = td.cloud().dispersion().update ( dt, cellI, U_, Uc_, UTurb_, tTurb_ ); Pc_ = td.PInterp().interpolate(this->position(), tetIs); //- myadd } template<class ParcelType> template<class TrackData> void Foam::KinematicParcel<ParcelType>::cellValueSourceCorrection ( TrackData& td, const scalar dt, const label cellI ) { Uc_ += td.cloud().UTrans()[cellI]/massCell(cellI); } template<class ParcelType> template<class TrackData> void Foam::KinematicParcel<ParcelType>::calc ( TrackData& td, const scalar dt, const label cellI ) { // Define local properties at beginning of time step // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ const scalar np0 = nParticle_; const scalar mass0 = mass(); // Reynolds number const scalar Re = this->Re(U_, d_, rhoc_, muc_); // Sources //~~~~~~~~ // Explicit momentum source for particle vector Su = vector::zero; // Linearised momentum source coefficient scalar Spu = 0.0; // Momentum transfer from the particle to the carrier phase vector dUTrans = vector::zero; //- myadd if(Ro_>0) { vector storeddRt = calcDiameter(td, dt, cellI, muc_, rhoc_, Pc_); this->dRt_=storeddRt[0]; this->d_=storeddRt[1]; } // Motion // ~~~~~~ // Calculate new particle velocity this->U_ = calcVelocity(td, dt, cellI, Re, muc_, mass0, Su, dUTrans, Spu); // Accumulate carrier phase source terms // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (td.cloud().solution().coupled()) { // Update momentum transfer td.cloud().UTrans()[cellI] += np0*dUTrans; // Update momentum transfer coefficient td.cloud().UCoeff()[cellI] += np0*Spu; } } template<class ParcelType> template<class TrackData> const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity ( TrackData& td, const scalar dt, const label cellI, const scalar Re, const scalar mu, const scalar mass, const vector& Su, vector& dUTrans, scalar& Spu ) const { typedef typename TrackData::cloudType cloudType; typedef typename cloudType::parcelType parcelType; typedef typename cloudType::forceType forceType; const forceType& forces = td.cloud().forces(); // Momentum source due to particle forces const parcelType& p = static_cast<const parcelType&>(*this); const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu); const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu); const forceSuSp Feff = Fcp + Fncp; const scalar massEff = forces.massEff(p, mass); // New particle velocity //~~~~~~~~~~~~~~~~~~~~~~ // Update velocity - treat as 3-D const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/massEff; const scalar bp = Feff.Sp()/massEff; Spu = dt*Feff.Sp(); IntegrationScheme<vector>::integrationResult Ures = td.cloud().UIntegrator().integrate(U_, dt, abp, bp); vector Unew = Ures.value(); // note: Feff.Sp() and Fc.Sp() must be the same dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su()); // Apply correction to velocity and dUTrans for reduced-D cases const polyMesh& mesh = td.cloud().pMesh(); meshTools::constrainDirection(mesh, mesh.solutionD(), Unew); meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans); return Unew; } //- myadd template<class ParcelType> template<class TrackData> const Foam::vector Foam::KinematicParcel<ParcelType>::calcDiameter ( TrackData& td, const scalar dt, const label cellI, const scalar mu, const scalar rhoc, const scalar Pc ) const { // Saturated vapor pressure const scalar pv=0; //surface tension const scalar sigma=0.0712; scalar R=d_/2; //The fourth-order Rungekuta method solves second-order differential equations const scalar k1=dRt_; const scalar h1=-3*pow(k1,2)/2/R-(4*mu*k1+2*sigma)/rhoc/pow(R,2)+(pv-Pc_+pgo_*pow(Ro_,4.2)*pow(R,-4.2))/rhoc/R; const scalar k2=k1+0.5*dt*h1; const scalar h2=-3*pow(k2,2)/2/(R+0.5*dt*k1)-(4*mu*k2+2*sigma)/rhoc/pow(R+0.5*dt*k1,2)+(pv-Pc_+pgo_*pow(Ro_,4.2)*pow(R+0.5*dt*k1,-4.2))/rhoc/(R+0.5*dt*k1); const scalar k3=k1+0.5*dt*h2; const scalar h3=-3*pow(k3,2)/2/(R+0.5*dt*k2)-(4*mu*k3+2*sigma)/rhoc/pow(R+0.5*dt*k2,2)+(pv-Pc_+pgo_*pow(Ro_,4.2)*pow(R+0.5*dt*k2,-4.2))/rhoc/(R+0.5*dt*k2); const scalar k4=k1+dt*h3; const scalar h4=-3*pow(k4,2)/2/(R+0.5*dt*k3)-(4*mu*k4+2*sigma)/rhoc/pow(R+0.5*dt*k3,2)+(pv-Pc_+pgo_*pow(Ro_,4.2)*pow(R+0.5*dt*k3,-4.2))/rhoc/(R+0.5*dt*k3); vector Rtvector(0,0,0); Rtvector[0]=k1+dt/6*(h1+2*h2+2*h3+h4); Rtvector[1]=2*(R+dt/6*(k1+2*k2+2*k3+k4)); return Rtvector; }
KinematicParcel.H
template<class ParcelType> class KinematicParcel : public ParcelType { // Private data //- Size in bytes of the fields static const std::size_t sizeofFields_; //- Number of particle tracking attempts before we assume that it stalls static label maxTrackAttempts; public: //- Class to hold kinematic particle constant properties class constantProperties { protected: // Protected data //- Constant properties dictionary const dictionary dict_; private: // Private data //- Parcel type id - used for post-processing to flag the type // of parcels issued by this cloud demandDrivenEntry<label> parcelTypeId_; //- Minimum density [kg/m3] demandDrivenEntry<scalar> rhoMin_; //- Particle density [kg/m3] (constant) demandDrivenEntry<scalar> rho0_; //- Minimum parcel mass [kg] demandDrivenEntry<scalar> minParcelMass_; //- myadd Particle initial radius [m] demandDrivenEntry<scalar> Ro0_; //- myadd Particle initial pressure inside [kg/m/s2] demandDrivenEntry<scalar> pgo0_; public: // Constructors //- Null constructor constantProperties(); //- Copy constructor constantProperties(const constantProperties& cp); //- Construct from dictionary constantProperties(const dictionary& parentDict); // Member functions //- Return const access to the constant properties dictionary inline const dictionary& dict() const; //- Return const access to the parcel type id inline label parcelTypeId() const; //- Return const access to the minimum density inline scalar rhoMin() const; //- Return const access to the particle density inline scalar rho0() const; //- Return const access to the minimum parcel mass inline scalar minParcelMass() const; //- myadd Return const access to the particle intitial radius inline scalar Ro0() const; //- myadd Return const access to the particle initial pressure inline scalar pgo0() const; }; template<class CloudType> class TrackingData : public ParcelType::template TrackingData<CloudType> { public: enum trackPart { tpVelocityHalfStep, tpLinearTrack, tpRotationalTrack }; private: // Private data // Interpolators for continuous phase fields //- Density interpolator autoPtr<interpolation<scalar> > rhoInterp_; //- Velocity interpolator autoPtr<interpolation<vector> > UInterp_; //- Dynamic viscosity interpolator autoPtr<interpolation<scalar> > muInterp_; //-myadd Pressure interpolator autoPtr<interpolation<scalar> > PInterp_; //- Local gravitational or other body-force acceleration const vector& g_; // label specifying which part of the integration // algorithm is taking place trackPart part_; public: // Constructors //- Construct from components inline TrackingData ( CloudType& cloud, trackPart part = tpLinearTrack ); // Member functions //- Return conat access to the interpolator for continuous // phase density field inline const interpolation<scalar>& rhoInterp() const; //- Return conat access to the interpolator for continuous // phase velocity field inline const interpolation<vector>& UInterp() const; //- Return conat access to the interpolator for continuous // phase dynamic viscosity field inline const interpolation<scalar>& muInterp() const; // Return const access to the gravitational acceleration vector inline const vector& g() const; //- Return the part of the tracking operation taking place inline trackPart part() const; //- Return access to the part of the tracking operation taking place inline trackPart& part(); //- myadd Return const access to the interpolator for continuous phase pressure field inline const interpolation<scalar>& PInterp() const; }; protected: // Protected data // Parcel properties //- Active flag - tracking inactive when active = false bool active_; //- Parcel type id label typeId_; //- Number of particles in Parcel scalar nParticle_; //- Diameter [m] scalar d_; //- Target diameter [m] scalar dTarget_; //- Velocity of Parcel [m/s] vector U_; //- Density [kg/m3] scalar rho_; //- Age [s] scalar age_; //- Time spent in turbulent eddy [s] scalar tTurb_; //- Turbulent velocity fluctuation [m/s] vector UTurb_; //- myadd change rate of bubble radius [m/s] scalar dRt_; //- myadd initial bubble radius [m] scalar Ro_; //- myadd initial bubble pressuer [Pa] scalar pgo_; // Cell-based quantities //- Density [kg/m3] scalar rhoc_; //- Velocity [m/s] vector Uc_; //- Viscosity [Pa.s] scalar muc_; //- myadd Pressure [Pa] scalar Pc_; // Protected Member Functions //- Calculate new particle velocity template<class TrackData> const vector calcVelocity ( TrackData& td, const scalar dt, // timestep const label cellI, // owner cell const scalar Re, // Reynolds number const scalar mu, // local carrier viscosity const scalar mass, // mass const vector& Su, // explicit particle momentum source vector& dUTrans, // momentum transfer to carrier scalar& Spu // linearised drag coefficient ) const; //- myadd Calculate new particle Diameter template<class TrackData> const vector calcDiameter ( TrackData& td, const scalar dt, const label cellI, const scalar mu, const scalar rhoc, const scalar Pc ) const; public: // Static data members //- Runtime type information TypeName("KinematicParcel"); //- String representation of properties AddToPropertyList ( ParcelType, " active" + " typeId" + " nParticle" + " d" + " dTarget " + " (Ux Uy Uz)" + " rho" + " age" + " tTurb" + " (UTurbx UTurby UTurbz)" //- myadd + " dRt" + " Ro" + " pgo" );
-
@dxl 首先要看你基于的哪个cloud做的,并不是所有的cloud都真实的包含了颗粒碰撞。要查一下code确认你的程序是不是真正的引用了碰撞模型。
其次,OF中的碰撞模型是O'Rourke的碰撞模型,是计算碰撞概率的,不是deteministic的(不知道中文怎么说)。即不是追踪所有颗粒,计算两个颗粒的路径,判断两个颗粒是否会发生碰撞的方法。
基于OF10的话,你可以在下面找到具体的源程序。
OpenFOAM-10/src/lagrangian/parcel/submodels/Spray/StochasticCollision /ORourkeCollision/