Euler-Lagrange的一些解析,sprayFoam
-
@星星星星晴 您好,您推荐的2. 输出文件方法我基本弄明白了,特别适合输出一个平面上不同时刻的统计数据,然后用python数据处理。
我遇到的问题主要是输出三维规则网格空间点上的水滴质量,把这个数据导入到别的程序中。我尝试过用paraview输出粒子所在空间点数据,因为很不规则,所以之后用python三维插值,但是效果很差。
所以我觉得您提出的1. 将拉格朗日场转到EULER场比较适合我的情况。但是关于这个代码使用还想请教您几个问题:
1.template<class CloudType> Foam::volScalarField& Foam::IVT8_StochasticCollisionModel_2022<CloudType>::PPC() { if (!PPCPtr_.valid() ) { const fvMesh& mesh = this->owner().mesh(); PPCPtr_.reset ( new volScalarField ( IOobject ( this->owner().name() + ":PPC", mesh.time().timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), mesh, dimensionedScalar ("zero", dimless, 0.0), zeroGradientFvPatchScalarField::typeName ) ); } return PPCPtr_(); }
这段代码应该是定义了一个函数,应该放在求解器的哪个文件中呢?
或者是自己重新定义一个头文件?-
这个函数中的PPC应该是一个通用模板吧?比如我想输出 p.mass()质量的数据到EULER场,这个PPC直接用mass替换?
-
调用的时候您说用下面的代码,我想知道这段代码又是放在对应的哪个文件中呢?
scalar& PPC = PPCPtr_->primitiveFieldRef() [cellI]; PPC += 1.0;
-
-
@chengan-wang
同样可以放到cloudfunction中,对应的改名什么的
第一步写这个function
第二步在constructor里面调用,生成这个field
第三步在你需要的地方调用这个方法比较复杂,需要多个地方定义,我这仅仅是提供一个思路,可以在任何地方生成这个field,只要在计算的时候能够赋值即可,你也可以在你的solver里面生成这个field,只要在拉格朗日部分能调用到这个field就行。
为什么要插值?我没get到你的点。插值也有很多方法的,另外当你的数据量不够的时候自然会存在很大的误差,这是统计学的问题了。
方法一天然的会损失位置信息,因为网格中储存的数据都是在网格中心,如果做复杂的运算的话,这样会损失掉很多信息,anyway这是你的计算,你可以试试。
我个人觉得方法二足矣应对所有的统计问题,无非就是后处理计算复杂一点,无法paraview中直接连续性的可视化
-
@chengan-wang 这也简单啊,同理,你不需要判断是否碰到面,而是做一个判断是否在你的区域内即可,比如你需要判断在两个面之间,z-dir
if(p.position(2)<1 & p.position(2)>0) { 输出 }
但是这种方法会出现什么问题,比如一个parcel在这个区域待了3个$\Delta t_L$, 那就输出3次,就是你的输出文件几何倍数的增长。。
-
@星星星星晴 在 Euler-Lagrange的一些解析,sprayFoam 中说:
@chengan-wang 这也简单啊,同理,你不需要判断是否碰到面,而是做一个判断是否在你的区域内即可,比如你需要判断在两个面之间,z-dir
if(p.position(2)<1 & p.position(2)>0) { 输出 }
但是这种方法会出现什么问题,比如一个parcel在这个区域待了3个$\Delta t_L$, 那就输出3次,就是你的输出文件几何倍数的增长。。
您好,如果加判断条件的话应该是在particleCollector代码中加,然后编译,对吧?这个圆柱面或者矩形面的定义就无所谓了。也就是说单独命名一个输出文件。
还有,我的网格数量挺大,200x100x100,我得定义好多个长方体区域,采用循环的办法,判断粒子与每一个区域的相对位置,如果符合条件就保存数据,跳出循环。但如果我的计算时间很长,每个时间点都保存数据,数据量肯定很大,能否只保存最后一个时间步的数据?也就是说我需要传递controlDict的参数,比如startFrom startTime; startTime 0; stopAt endTime; endTime 5; deltaT 0.01;
-
-
-
@chengan-wang 自然可以 这就需要你自己研究啦,也有许多其他的方法,比如你获得你要查的目标的cell number,然后只搞这部分就行。
-
@chengan-wang
不能给你我的code,这个不是我自己的私产,我只能大致的说一下方法。
祝你好运。 -
@chengan-wang
正如我前面所说的,你可以在任何Euler/Lagrange场的计算的地方建立这个field,solver部分也可以。
然后在任何Lagrange计算的地方对这个field进行操作。
只要是这个计算的方法是符合你的要求的,Euler 时间步还是Lagrange 时间步。
你可以参考具体下面的拉格朗日部分的流程图,在你需要的地方放置计算即可。 -
@chengan-wang
ivt-xxxx 这个是我们自己的cloudfunction的名字而已,
ppc也就是这个cloudfunction中的一个function的名字,可以被调用。你可以尝试在solver中建立一个field,叫任何名字都可以。
-
@星星星星晴 前辈,您好!我在使用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();感觉有问题,从结果上看感觉没有把相同位置的粒子质量加到一起。
想麻烦您帮忙看看,谢谢