30 m_log <<
" Initialising spray model" << endl;
32 s_backPtr = _particle;
33 s_lengthFactor = s_backPtr->
m_partList[0].s_lengthFactor;
35 m_activePrimaryBUp = s_backPtr->m_activePrimaryBUp;
36 m_activeSecondaryBUp = s_backPtr->m_activeSecondaryBUp;
38 if(!s_backPtr->m_restart) {
39 s_backPtr->m_particleResiduum = 0.0;
43 mAlloc(m_needleLiftTime, 2,
"needleLiftTime", AT_);
44 mAlloc(m_sprayAngle, 2,
"m_sprayAngle", F0, AT_);
45 mAlloc(m_injData, m_injDataSize,
"m_injData", AT_);
47 readSprayProperties();
49 logDistributionFunction();
61 m_injStopTimeStep = -1;
62 m_injStopTimeStep = Context::getSolverProperty<MInt>(
"injStopTimeStep", solverId(), AT_, &m_injStopTimeStep);
64 m_injStartTimeStep = -1;
65 m_injStartTimeStep = Context::getSolverProperty<MInt>(
"injStartTimeStep", solverId(), AT_, &m_injStartTimeStep);
67 if(m_activePrimaryBUp) {
68 m_needleLiftTime[0] = 0.01;
69 m_needleLiftTime[1] = 0.02;
79 m_needleLiftTime[0] = Context::getSolverProperty<MFloat>(
"injectorNeedleLiftTime", solverId(), AT_, 0);
80 m_needleLiftTime[1] = Context::getSolverProperty<MFloat>(
"injectorNeedleLiftTime", solverId(), AT_, 1);
81 m_useNeedleLiftTime =
true;
91 m_injDuration = 0.00055;
92 m_injDuration = Context::getSolverProperty<MFloat>(
"injectorInjectionTime", solverId(), AT_, &m_injDuration);
102 m_injectionCrankAngle = -1;
103 m_injectionCrankAngle =
104 Context::getSolverProperty<MFloat>(
"injectionCrankAngle", solverId(), AT_, &m_injectionCrankAngle);
108 m_Strouhal = Context::getSolverProperty<MFloat>(
"Strouhal", solverId(), AT_, &m_Strouhal);
110 m_initialCad = Context::getSolverProperty<MFloat>(
"initialCrankAngle", solverId(), AT_, &m_initialCad);
119 m_primRosinRammler =
true;
121 Context::getSolverProperty<MBool>(
"sprayPrimRosinRammler", solverId(), AT_, &m_primRosinRammler);
132 MString injectorType =
"FULLCONE";
133 injectorType = Context::getSolverProperty<MString>(
"injectorType", solverId(), AT_, &injectorType);
137 m_multiHoleInjector =
false;
138 m_singleHoleInjector =
false;
139 m_hollowConeInjector =
false;
141 switch(m_injectorType) {
145 m_multiHoleInjector =
true;
149 m_singleHoleInjector =
true;
153 m_hollowConeInjector =
true;
157 mTerm(1, AT_,
"Unknown injector-type!");
168 m_injectorNozzleDiameter = Context::getSolverProperty<MFloat>(
"injectorNozzleDiameter", solverId(), AT_);
171 if(m_multiHoleInjector) {
182 mAlloc(m_injectorDiameter, 1,
"injectorDiameter", AT_);
183 m_injectorDiameter[0] = Context::getSolverProperty<MFloat>(
"injectorDiameter", solverId(), AT_);
186 mAlloc(m_injectorDiameter, length + 1,
"injectorDiameter", AT_);
187 m_injectorDiameter[0] = 0.0;
188 for(
MInt i = 0; i < length; i++) {
189 m_injectorDiameter[i + 1] = Context::getSolverProperty<MFloat>(
"injectorDiameter", solverId(), AT_, i);
190 m_injectorDiameter[0] =
mMax(m_injectorDiameter[0], m_injectorDiameter[i + 1]);
203 mAlloc(m_injectorOrificeAngle, 1,
"injectorOrificeAngle", AT_);
204 m_injectorOrificeAngle[0] = 37.0;
205 m_injectorOrificeAngle[0] =
206 Context::getSolverProperty<MFloat>(
"injectorOrificeAngle", solverId(), AT_, &m_injectorOrificeAngle[0]);
209 mAlloc(m_injectorOrificeAngle, length,
"injectorOrificeAngle", AT_);
210 for(
MInt i = 0; i < length; i++) {
211 m_injectorOrificeAngle[i] = Context::getSolverProperty<MFloat>(
"injectorOrificeAngle", solverId(), AT_, i);
224 mAlloc(m_orificePositionAngle, length,
"orificePositionAngle", AT_);
225 for(
MInt i = 0; i < length; i++) {
226 m_orificePositionAngle[i] = Context::getSolverProperty<MFloat>(
"orificePositionAngle", solverId(), AT_, i);
238 m_injectorOrficeSize = Context::getSolverProperty<MFloat>(
"injectorOrficeSize", solverId(), AT_);
241 if(m_primRosinRammler && m_singleHoleInjector) {
248 m_RosinRammlerMin = 90;
249 m_RosinRammlerMean = 15;
250 m_RosinRammlerMax = 5;
251 m_RosinRammlerSpread = 3;
253 }
else if(m_primRosinRammler && m_hollowConeInjector) {
256 m_RosinRammlerMin = 3.33333333;
257 m_RosinRammlerMean = 1;
258 m_RosinRammlerMax = 0.33333333;
259 m_RosinRammlerSpread = 3;
260 }
else if(m_primRosinRammler && m_multiHoleInjector) {
261 m_RosinRammlerMin = 3.33333333;
262 m_RosinRammlerMean = 1;
263 m_RosinRammlerMax = 0.33333333;
264 m_RosinRammlerSpread = 3;
268 if(m_primRosinRammler) {
276 m_RosinRammlerSpread =
277 Context::getSolverProperty<MFloat>(
"RosinRammlerSpread", solverId(), AT_, &m_RosinRammlerSpread);
286 m_RosinRammlerMin = Context::getSolverProperty<MFloat>(
"RosinRammlerMin", solverId(), AT_, &m_RosinRammlerMin);
295 m_RosinRammlerMax = Context::getSolverProperty<MFloat>(
"RosinRammlerMax", solverId(), AT_, &m_RosinRammlerMax);
304 m_RosinRammlerMean = Context::getSolverProperty<MFloat>(
"RosinRammlerMean", solverId(), AT_, &m_RosinRammlerMean);
314 m_injectionDistSigmaCoeff =
315 Context::getSolverProperty<MFloat>(
"spawnDistSigmaCoeff", solverId(), AT_, &m_injectionDistSigmaCoeff);
324 MString partDist =
"PART_EMITT_DIST_NONE";
325 m_partEmittDist =
string2enum(Context::getSolverProperty<MString>(
"partEmittDist", solverId(), AT_, &partDist));
327 MFloat diameter = m_injectorNozzleDiameter;
328 if(m_multiHoleInjector) {
329 diameter = m_injectorDiameter[0];
330 }
else if(m_singleHoleInjector) {
331 diameter = m_injectorOrficeSize;
334 if(diameter < (2 * s_backPtr->c_cellLengthAtLevel(s_backPtr->minLevel()))) {
335 m_broadcastInjected =
false;
336 if(s_backPtr->domainId() == 0) {
337 cerr <<
"Injection broadcast not necessary!" << endl;
348 m_injectionSpeed = Context::getSolverProperty<MFloat>(
"injectionSpeed", solverId(), AT_, &m_injectionSpeed);
350 ASSERT(m_injectionSpeed > 0,
"");
355 if(numInjections == 1) {
356 m_currentInjectionRate = Context::getSolverProperty<MFloat>(
"sprayInjectionRate", solverId(), AT_);
359 m_injectionRateList.resize(numInjections);
360 m_injectionTimings.resize(numInjections);
362 for(
MInt i = 0; i < numInjections; i++) {
370 m_injectionRateList[i] = Context::getSolverProperty<MFloat>(
"sprayInjectionRate", solverId(), AT_, i);
379 m_injectionTimings[i] = Context::getSolverProperty<MFloat>(
"sprayInjectionTiming", solverId(), AT_, i);
383 ASSERT(!m_activePrimaryBUp,
"");
394 for(
MInt i = 0; i < nDim; i++) {
395 m_injectorDir[i] = 0;
396 if(i == nDim - 1) m_injectorDir[i] = 1.0;
401 TERMM(1,
"Need to give a Coordinate for every dimension");
404 for(
MInt i = 0; i < nDim; i++) {
405 m_injectorDir[i] = Context::getSolverProperty<MFloat>(
"injectorDir", solverId(), AT_, i);
409 if(m_hollowConeInjector) {
417 m_maxNoPrimParcels = 1;
419 Context::getSolverProperty<MInt>(
"sprayPrimaryMaxPrimParcels", solverId(), AT_, &m_maxNoPrimParcels);
428 m_primParcelsPerInj = 1000000;
429 m_primParcelsPerInj =
430 Context::getSolverProperty<MInt>(
"sprayPrimaryParcelsPerInj", solverId(), AT_, &m_primParcelsPerInj);
439 m_angularGap = Context::getSolverProperty<MFloat>(
"injectorAngularGap", solverId(), AT_, &m_angularGap);
442 if(m_singleHoleInjector || m_multiHoleInjector) {
452 m_primBrkUpParcelSize = 1;
453 m_primBrkUpParcelSize =
454 Context::getSolverProperty<MInt>(
"sprayPrimaryBUpParcelSize", solverId(), AT_, &m_primBrkUpParcelSize);
455 ASSERT(m_primBrkUpParcelSize > 0,
"ERROR: Invalid parcel size.");
465 m_sprayAngle[0] = Context::getSolverProperty<MFloat>(
"sprayAngle", solverId(), AT_, &m_sprayAngle[0]);
466 m_sprayAngle[1] = Context::getSolverProperty<MFloat>(
"sprayAngle", solverId(), AT_, &m_sprayAngle[0]);
468 m_sprayAngle[0] = Context::getSolverProperty<MFloat>(
"sprayAngle", solverId(), AT_, 0);
469 m_sprayAngle[1] = Context::getSolverProperty<MFloat>(
"sprayAngle", solverId(), AT_, 1);
473 m_sprayConeAngle = Context::getSolverProperty<MFloat>(
"sprayConeAngle", solverId(), AT_, &m_sprayConeAngle);
484 m_spraySeed = (
MLong)Context::getSolverProperty<MInt>(
"spawnSeed", solverId(), AT_);
487 if(m_activePrimaryBUp) m_PRNGPrimBreakUp.seed(m_spraySeed);
488 if(m_activeSecondaryBUp) m_PRNGSecBU.seed(m_spraySeed);
493 if(m_activeSecondaryBUp) {
501 m_secBUDisplayMessage =
false;
502 m_secBUDisplayMessage =
503 Context::getSolverProperty<MBool>(
"spraysecBUPOutput", solverId(), AT_, &m_secBUDisplayMessage);
505 m_RTsecBreakUp =
true;
506 m_RTsecBreakUp = Context::getSolverProperty<MBool>(
"RayleighTaylorBreakUp", solverId(), AT_, &m_RTsecBreakUp);
508 m_KHsecBreakUp =
true;
509 m_KHsecBreakUp = Context::getSolverProperty<MBool>(
"KelvinHelmholtzBreakUp", solverId(), AT_, &m_KHsecBreakUp);
511 m_RTsecBreakUpLength =
true;
512 m_RTsecBreakUpLength =
513 Context::getSolverProperty<MBool>(
"RayleighTaylorBreakUpLength", solverId(), AT_, &m_RTsecBreakUpLength);
527 m_predictivePRNG =
true;
528 m_predictivePRNG = Context::getSolverProperty<MBool>(
"predictablePRNGSecBU", solverId(), AT_, &m_predictivePRNG);
545 m_sprayAngleKH = Context::getSolverProperty<MFloat>(
"sprayAngleKH", solverId(), AT_, &m_sprayAngleKH);
561 m_weberLimitRT = Context::getSolverProperty<MFloat>(
"sprayWeberLimitRT", solverId(), AT_, &m_weberLimitRT);
570 m_B0 = Context::getSolverProperty<MFloat>(
"sprayB0", solverId(), AT_, &m_B0);
580 m_B1 = Context::getSolverProperty<MFloat>(
"sprayB1", solverId(), AT_, &m_B1);
592 m_CRT = Context::getSolverProperty<MFloat>(
"sprayCRT", solverId(), AT_, &m_CRT);
601 m_CT = Context::getSolverProperty<MFloat>(
"sprayCT", solverId(), AT_, &m_CT);
610 m_weberLimitKH = Context::getSolverProperty<MFloat>(
"sprayWeberLimitKH", solverId(), AT_, &m_weberLimitKH);
621 m_massLimit = Context::getSolverProperty<MFloat>(
"sprayMassLimit", solverId(), AT_, &m_massLimit);
630 m_maxRTChildParcels = Context::getSolverProperty<MInt>(
"sprayMaxRTChilds", solverId(), AT_, &m_maxRTChildParcels);
645 m_RTDiameterMode = Context::getSolverProperty<MInt>(
"sprayRTDiameterMode", solverId(), AT_, &m_RTDiameterMode);
659 m_Cbu = Context::getSolverProperty<MFloat>(
"sprayCbu", solverId(), AT_, &m_Cbu);
661 MFloat breakUpLength = m_Cbu * m_injectorNozzleDiameter * sqrt(material().ambientDensityRatio());
663 cerr0 <<
"Liquid break-up length for KH-RT model is " << breakUpLength << endl;
675 ASSERT(s_backPtr->a_noParticles() == 0,
"");
676 if(domainId() == 0) {
677 cerr <<
"Before injection" << m_injectionCrankAngle <<
" "
680 for(
MInt i = 0; i < m_injDataSize; i++) {
686 m_timeSinceSOI += dt;
689 MBool endOfInjection =
false;
690 if(m_injStopTimeStep > -1 &&
globalTimeStep > m_injStopTimeStep) endOfInjection =
true;
691 if(m_timeSinceSOI > m_injDuration) endOfInjection =
true;
692 if(m_injStartTimeStep > -1 &&
globalTimeStep < m_injStartTimeStep) endOfInjection =
true;
696 cerr0 <<
"After injection " << m_timeSinceSOI <<
" " << m_injDuration << endl;
698 for(
MInt i = 0; i < m_injDataSize; i++) {
701 if(s_backPtr->m_spawnCellId > -1) {
702 m_injData[0] = m_timeSinceSOI;
709 if(s_backPtr->m_spawnCellId < 0) {
710 ASSERT(s_backPtr->noDomains() > 1,
"");
713 if(m_broadcastInjected) {
714 s_backPtr->recvInjected();
716 for(
MInt i = 0; i < m_injDataSize; i++) {
722 ASSERT(domainId() == s_backPtr->m_spawnDomainId,
"");
725 updateInjectionRate();
728 const MInt prevNo = s_backPtr->a_noParticles();
733 if(m_broadcastInjected) {
734 s_backPtr->broadcastInjected(prevNo);
747 ASSERT(domainId() == s_backPtr->m_spawnDomainId,
"");
748 ASSERT(s_backPtr->m_spawnCellId > -1,
"");
751 MInt noInjParticles = -1;
752 MInt noInjDroplets = -1;
754 MFloat injParticleDiameter = NAN;
759 MFloat injectionProgress = m_timeSinceSOI / m_injDuration;
766 if(injectionProgress < m_needleLiftTime[0]) {
767 rampFactor = injectionProgress / m_needleLiftTime[0];
768 }
else if(injectionProgress > (1 - m_needleLiftTime[1])) {
769 rampFactor = 1 - (injectionProgress - (1 - m_needleLiftTime[1])) / m_needleLiftTime[1];
774 injMass = m_currentInjectionRate * rampFactor * dt + s_backPtr->m_particleResiduum;
777 MInt noDroplets =
mMax(1, (
MInt)(m_primBrkUpParcelSize * rampFactor));
781 if(m_singleHoleInjector) {
783 if(m_primRosinRammler) {
786 const MFloat diameterMin = m_injectorNozzleDiameter / m_RosinRammlerMin;
787 const MFloat diameterMean = m_injectorNozzleDiameter / m_RosinRammlerMean;
788 const MFloat diameterMax = m_injectorNozzleDiameter / m_RosinRammlerMax;
791 const MFloat minMass = sphericalMass(diameterMin, material().T());
796 MFloat remainingMass = injMass;
799 while(remainingMass > minMass) {
801 injParticleDiameter = lastD;
803 injParticleDiameter =
804 rosinRammler(diameterMin, diameterMean, diameterMax, m_RosinRammlerSpread, randomPrimBreakUp(1));
807 const MFloat particleMass = noDroplets * sphericalMass(injParticleDiameter, material().T());
808 const MFloat tempMass = remainingMass - particleMass;
812 lastD = injParticleDiameter;
813 s_backPtr->m_particleResiduum = remainingMass;
819 const MFloat coneDiameter = m_injectorOrficeSize - 0.5 * injParticleDiameter;
821 injectParticles(1, s_backPtr->m_spawnCoord, injParticleDiameter, material().densityRatio(),
822 m_injectionSpeed * rampFactor, m_sprayConeAngle, coneDiameter,
false, 0.0, noDroplets, 1);
825 noInjDroplets += noDroplets;
826 injMass += particleMass;
827 remainingMass -= particleMass;
833 injParticleDiameter = m_injectorNozzleDiameter;
836 const MFloat initParticleMass = sphericalMass(injParticleDiameter, material().T());
839 MFloat noNewParticles = (m_currentInjectionRate * dt) / initParticleMass + s_backPtr->m_particleResiduum;
840 noInjParticles = floor(noNewParticles);
842 s_backPtr->m_particleResiduum = noNewParticles - noInjParticles;
844 injectParticles(noInjParticles, s_backPtr->m_spawnCoord, injParticleDiameter, material().densityRatio(),
845 m_injectionSpeed, m_sprayConeAngle, 0.0,
false, 0.0, 1, 1);
847 }
else if(m_hollowConeInjector) {
853 const MInt parcelTargetNo = m_primParcelsPerInj / m_injDuration * dt;
856 noInjParticles = std::min(parcelTargetNo, m_maxNoPrimParcels);
859 const MFloat scaledInjectionV = rampFactor * m_injectionSpeed;
862 MFloat initialLiquidSheetThickness = NAN;
864 if(m_angularGap > -1) {
865 initialLiquidSheetThickness = max(m_angularGap * rampFactor, m_primMinDropletSize);
867 initialLiquidSheetThickness =
868 max(0.5 * m_injectorNozzleDiameter
869 - sqrt(M_PI *
POW2(m_injectorNozzleDiameter) * material().density() * scaledInjectionV
870 - 4.0 * m_currentInjectionRate * rampFactor)
871 / (2.0 * sqrt(M_PI) * sqrt(material().density() * scaledInjectionV)),
872 m_primMinDropletSize);
877 injParticleDiameter = 4.0 / M_PI * initialLiquidSheetThickness;
878 if(m_primRosinRammler) {
879 const MFloat diameterMin = injParticleDiameter / m_RosinRammlerMin;
880 const MFloat diameterMean = injParticleDiameter / m_RosinRammlerMean;
881 const MFloat diameterMax = injParticleDiameter / m_RosinRammlerMax;
882 injParticleDiameter =
883 rosinRammler(diameterMin, diameterMean, diameterMax, m_RosinRammlerSpread, randomPrimBreakUp(1));
887 const MFloat initialDropletMass = sphericalMass(injParticleDiameter, material().T());
889 const MFloat newDroplets = injMass / initialDropletMass;
891 noInjDroplets = floor(newDroplets / noInjParticles);
893 if(noInjDroplets > 0) {
894 const MFloat injectorGapDiameter = m_injectorNozzleDiameter - 0.5 * injParticleDiameter;
895 s_backPtr->m_particleResiduum = (newDroplets - noInjParticles * noInjDroplets) * initialDropletMass;
897 injectParticles(noInjParticles, s_backPtr->m_spawnCoord, injParticleDiameter, material().densityRatio(),
898 scaledInjectionV, m_sprayConeAngle, injectorGapDiameter * s_lengthFactor,
true, m_angularGap,
901 injMass = noInjParticles * noInjDroplets * initialDropletMass;
908 }
else if(m_multiHoleInjector) {
915 switch(m_injectorType) {
929 mTerm(1, AT_,
"Unknown multi-hole injector-type!");
932 static constexpr MFloat C_a = 0.65;
936 const MFloat initialConeAngle = m_sprayConeAngle;
944 ASSERT(m_primRosinRammler,
"");
946 const MFloat effectiveDiam = sqrt(C_a) * m_injectorNozzleDiameter;
947 const MFloat diameterMin = effectiveDiam / m_RosinRammlerMin;
948 const MFloat diameterMean = effectiveDiam / m_RosinRammlerMean;
949 const MFloat diameterMax = effectiveDiam / m_RosinRammlerMax;
951 MFloat remainingMass = injMass;
952 const MFloat minMass = sphericalMass(diameterMin, material().T());
959 const MFloat angleDelta = 2 * M_PI / 8;
960 switch(m_injectorType) {
962 return (
id * angleDelta);
972 return id * angleDelta;
975 return m_orificePositionAngle[
id] / 180 * M_PI;
978 mTerm(1, AT_,
"Unknown multi-hole injector-type!");
983 switch(m_injectorType) {
986 return m_injectorDiameter[0];
989 return m_injectorDiameter[
id + 1];
992 mTerm(1, AT_,
"Unknown multi-hole injector-type!");
998 switch(m_injectorType) {
1004 return m_injectorOrificeAngle[
id];
1007 mTerm(1, AT_,
"Unknown multi-hole injector-type!");
1011 while(remainingMass > minMass) {
1013 injParticleDiameter = lastD;
1015 injParticleDiameter =
1016 rosinRammler(diameterMin, diameterMean, diameterMax, m_RosinRammlerSpread, randomPrimBreakUp(1));
1019 const MFloat particleMass = noHoles * noDroplets * sphericalMass(injParticleDiameter, material().T());
1021 const MFloat tempMass = remainingMass - particleMass;
1024 lastD = injParticleDiameter;
1025 s_backPtr->m_particleResiduum = remainingMass;
1030 injectParticles(noHoles, s_backPtr->m_spawnCoord, injParticleDiameter, material().densityRatio(),
1031 m_injectionSpeed, initialConeAngle, m_injectorOrficeSize, holePosition, injectorDiameter,
1032 holeAngle,
false, noDroplets, noHoles);
1034 injMass += particleMass;
1035 noInjParticles += noHoles;
1036 noInjDroplets += noHoles * noDroplets;
1037 remainingMass -= particleMass;
1042 m_injData[0] = m_timeSinceSOI;
1043 m_injData[1] = injectionProgress;
1044 m_injData[2] = rampFactor;
1045 m_injData[3] = noInjParticles;
1046 m_injData[4] = noInjDroplets;
1047 m_injData[5] = injParticleDiameter;
1048 m_injData[6] = injMass;
1050 array<MFloat, 3> injMomentum = {};
1051 for(
MInt i = 0; i < 3; i++) {
1055 for(
MInt i = 0; i < noInjParticles; i++) {
1056 const MInt id = s_backPtr->a_noParticles() - 1 - i;
1057 const MFloat mass = sphericalMass(s_backPtr->m_partList[
id].m_diameter, s_backPtr->m_partList[
id].m_temperature)
1058 * s_backPtr->m_partList[
id].m_noParticles;
1059 MFloat velMagSquared = 0;
1060 for(
MInt j = 0; j < nDim; j++) {
1061 injMomentum[j] += mass * s_backPtr->m_partList[
id].m_velocity[j];
1062 velMagSquared +=
POW2(s_backPtr->m_partList[
id].m_velocity[j]);
1064 injEnergy += mass * s_backPtr->m_material->cp(s_backPtr->m_partList[
id].m_temperature)
1065 * s_backPtr->m_partList[
id].m_temperature * 1 / s_backPtr->m_material->gammaMinusOne();
1066 injEnergy += 0.5 * mass * velMagSquared;
1068 m_injData[7] = injMomentum[0];
1069 m_injData[8] = injMomentum[1];
1070 m_injData[9] = injMomentum[2];
1071 m_injData[10] = injEnergy;
1073 cerr <<
globalTimeStep <<
" " << m_timeSinceSOI <<
" injMass " << injMass <<
" injParticles " << noInjParticles
1084 m_noRTsecBreakUp = 0;
1085 m_noKHsecBreakUp = 0;
1087 if(!m_RTsecBreakUp && !m_KHsecBreakUp)
return;
1090 const MInt noPart = s_backPtr->a_noParticles();
1093 if(0.5 * s_backPtr->m_partList.capacity() < s_backPtr->a_noParticles()) {
1094 s_backPtr->m_partList.reserve(10 * s_backPtr->m_partList.capacity());
1102 const MFloat breakupLength =
1103 m_RTsecBreakUpLength ? m_Cbu * m_injectorNozzleDiameter * sqrt(material().ambientDensityRatio()) : 0;
1111 ASSERT(m_injectorNozzleDiameter > 0,
"ERROR: No valid nozzle diameter set.");
1112 ASSERT(breakupLength > -MFloatEps,
"ERROR: Invalid breakup length.");
1115 for(
MInt i = 0; i < noPart; i++) {
1116 auto& droplet = s_backPtr->m_partList.at(i);
1118 if(droplet.firstStep())
continue;
1119 if(droplet.isInvalid())
continue;
1120 if(droplet.fullyEvaporated())
continue;
1121 ASSERT(!isnan(droplet.m_diameter),
"Invalid diameter! ");
1122 if(droplet.m_diameter < s_backPtr->m_sizeLimit)
continue;
1123 if(droplet.hadWallColl())
continue;
1125 ASSERT(droplet.m_cellId > 0,
"ERROR: Invalid cellId");
1128 droplet.m_breakUpTime += dt;
1130 ASSERT(!isnan(droplet.m_shedDiam) && droplet.m_shedDiam > -MFloatEps,
"");
1137 MFloat liquidLength = MFloatMax;
1138 if(m_RTsecBreakUpLength) {
1153 const MFloat dropRadius = 0.5 * droplet.m_diameter;
1154 const MFloat relV = droplet.magRelVel(&droplet.m_fluidVel[0], &droplet.m_velocity[0]);
1158 cerr << droplet.m_fluidVel[0] <<
" " << droplet.m_fluidVel[1] <<
" " << droplet.m_fluidVel[nDim - 1] << endl;
1165 droplet.WeberNumber(s_backPtr->m_material->density(droplet.m_temperature), pow(relV, 2), droplet.m_temperature)
1169 if(We_l < m_weberLimitKH)
continue;
1174 * droplet.particleRe(material().density(droplet.m_temperature), relV,
1175 material().dynamicViscosity(droplet.m_temperature))
1179 const MFloat We_g = droplet.WeberNumber(droplet.m_fluidDensity, pow(relV, 2), droplet.m_temperature) * droplet.s_We;
1181 const MFloat Z = sqrt(We_l) / Re_l;
1183 const MFloat T = Z * sqrt(We_g);
1186 const MFloat KH_waveL = 9.02 * dropRadius * (1.0 + 0.45 * sqrt(Z)) * (1.0 + 0.4 * pow(T, 0.7))
1187 / pow(1.0 + 0.865 * pow(We_g, 1.67), 0.6);
1191 const MFloat KH_diameter = 2.0 * m_B0 * KH_waveL;
1195 if(KH_diameter < 2 * MFloatEps) {
1196 cerr << KH_diameter <<
" d_KH " << dropRadius <<
" " << Z <<
" " << T <<
" " << We_g <<
" " << We_l << endl;
1200 const MFloat KH_growthR = ((0.34 + 0.38 * pow(We_g, 1.5)) / ((1.0 + Z) * (1.0 + 1.4 * pow(T, 0.6))))
1201 * sqrt(material().spraySurfaceTension(droplet.m_temperature)
1202 / (material().density(droplet.m_temperature) * pow(dropRadius, 3)))
1203 * sqrt(1 / droplet.s_We);
1207 const MFloat KH_time = 3.726 * m_B1 * dropRadius / (KH_waveL * KH_growthR);
1213 const MFloat magVel = droplet.magVel();
1214 array<MFloat, nDim> parentTrajectory{};
1216 for(
MInt v = 0; v < nDim; v++) {
1217 parentTrajectory.at(v) = droplet.m_velocity.at(v) / magVel;
1221 if(magVel < MFloatEps) {
1240 for(
MInt dir = 0; dir < nDim; dir++) {
1241 drop_accel +=
POW2(droplet.m_accel.at(dir));
1243 drop_accel = sqrt(drop_accel);
1246 const MFloat oldMass = droplet.sphericalMass() * droplet.m_noParticles;
1247 vector<MFloat> oldMom(nDim);
1248 for(
MInt n = 0; n < nDim; n++) {
1249 oldMom[n] = oldMass * droplet.m_velocity[n];
1251 const MFloat oldDiameter = droplet.m_diameter;
1261 * sqrt(3.0 * material().spraySurfaceTension(droplet.m_temperature)
1262 / (drop_accel * (material().density(droplet.m_temperature) - droplet.m_fluidDensity)))
1263 / sqrt(droplet.s_We);
1266 const MFloat RT_frequency =
1267 sqrt(2.0 / (3.0 * sqrt(3.0 * material().spraySurfaceTension(droplet.m_temperature)))
1268 * pow(drop_accel * (material().density(droplet.m_temperature) - droplet.m_fluidDensity), 3.0 / 2.0)
1269 / (material().density(droplet.m_temperature) + droplet.m_fluidDensity))
1270 * pow(droplet.s_We, 1 / 4);
1272 const MFloat RT_time = m_CT / RT_frequency;
1275 const MFloat RT_diameter = m_CRT * RT_waveL;
1277 ASSERT(RT_diameter > 0,
"");
1283 if(m_RTsecBreakUp && liquidLength >= breakupLength && droplet.m_breakUpTime >= RT_time
1284 && droplet.m_diameter > RT_diameter && We_l < m_weberLimitRT) {
1286 MFloat childDiameter = RT_diameter;
1287 if(m_RTDiameterMode == 2) {
1288 initPRNG(
globalTimeStep, s_backPtr->c_globalId(droplet.m_cellId));
1289 const MFloat diameterMin = RT_diameter / m_RosinRammlerMin;
1290 const MFloat diameterMean = RT_diameter / m_RosinRammlerMean;
1291 const MFloat diameterMax = RT_diameter / m_RosinRammlerMax;
1292 childDiameter =
rosinRammler(diameterMin, diameterMean, diameterMax, m_RosinRammlerSpread, randomSecBreakUp());
1293 if(droplet.m_diameter < childDiameter) {
1297 }
else if(m_RTDiameterMode == 1) {
1298 childDiameter = pow(oldDiameter, 2 / 3) * pow(RT_diameter, 1 / 3);
1301 if(m_maxRTChildParcels == 0) {
1303 auto noDroplets =
static_cast<MInt>(droplet.m_noParticles *
POW3(oldDiameter / childDiameter));
1308 if(noDroplets > 0) {
1309 childDiameter = pow(6.0 / PI, 1.0 / 3.0)
1310 * pow(oldMass / (noDroplets * material().density(droplet.m_temperature)), 1.0 / 3.0);
1311 droplet.m_diameter = childDiameter;
1312 droplet.m_shedDiam = childDiameter;
1313 droplet.m_noParticles = noDroplets;
1314 droplet.m_breakUpTime = 0.0;
1317 const MFloat newMass = sphericalMass(droplet.m_diameter, droplet.m_temperature) * droplet.m_noParticles;
1318 if(fabs(newMass - oldMass) > MFloatEps) {
1319 cerr <<
"Missing mass is RT-breakup " << oldMass <<
" " << newMass << endl;
1324 if(m_RTDiameterMode != 2) {
1325 initPRNG(
globalTimeStep, s_backPtr->c_globalId(droplet.m_cellId));
1329 const MFloat dropletVolume = pow(droplet.m_diameter, 3.0);
1330 const MFloat dropletRTVolume = pow(childDiameter, 3.0);
1332 MInt noBrokenUpDroplets = dropletVolume / dropletRTVolume;
1335 if(noBrokenUpDroplets == 1) {
1338 noBrokenUpDroplets *= droplet.m_noParticles;
1341 auto dropsPerParcel =
static_cast<MInt>(
1342 noBrokenUpDroplets > m_maxRTChildParcels ? floor(noBrokenUpDroplets / m_maxRTChildParcels) : 1);
1345 MInt noNewParcels = noBrokenUpDroplets;
1346 if(dropsPerParcel > 1) {
1347 noNewParcels = m_maxRTChildParcels - 1;
1352 const MFloat RT_dropletMass = sphericalMass(childDiameter, droplet.m_temperature);
1353 const MFloat parentMass = oldMass - dropsPerParcel * RT_dropletMass * noNewParcels;
1354 const MFloat parentDiam =
1355 pow(6.0 / PI, 1.0 / 3.0)
1356 * pow(parentMass / (dropsPerParcel * material().density(droplet.m_temperature)), 1.0 / 3.0);
1358 ASSERT(parentMass > 0,
"Invalid mass!");
1359 ASSERT(parentDiam > 0 && !isnan(parentDiam),
"");
1360 ASSERT(dropsPerParcel > 0,
"");
1363 droplet.m_diameter = parentDiam;
1364 droplet.m_shedDiam = parentDiam;
1365 droplet.m_noParticles = dropsPerParcel;
1366 droplet.m_breakUpTime = 0.0;
1372 if(m_secBUDisplayMessage) {
1373 cerr <<
"#################################################" << endl;
1374 cerr <<
"Break-Up type: RT" << endl;
1375 cerr <<
"parent diameter: " << droplet.m_diameter << endl;
1376 cerr <<
"Time: " << m_timeSinceSOI << endl;
1377 cerr <<
"old number of particles: " << droplet.m_noParticles << endl;
1378 cerr <<
"new particles/parcels: " << noNewParcels << endl;
1379 cerr <<
"parcel size: " << dropsPerParcel << endl;
1380 cerr <<
"Weber number: " << We_l << endl;
1381 cerr <<
"#################################################" << endl;
1384 MFloat sprayAngle = m_sprayAngle[0];
1385 if(m_hollowConeInjector) {
1388 static constexpr MFloat angularGapAngle = 45.0;
1389 sprayAngle = m_sprayAngle[0] - 2 * angularGapAngle;
1393 vector<MFloat> mom(nDim);
1397 array<MFloat, nDim> childVelocity{};
1398 for(
MInt s = 0; s < noNewParcels; s++) {
1399 randomVectorInCone(&childVelocity[0], &parentTrajectory[0], magVel, sprayAngle, m_partEmittDist,
1400 randomSecBreakUp(), m_injectionDistSigmaCoeff);
1402 const MInt childId = s_backPtr->addParticle(droplet.m_cellId, childDiameter, droplet.m_densityRatio, 0, 3,
1403 &childVelocity[0], &droplet.m_position[0], dropsPerParcel);
1406 s_backPtr->m_partList[childId].m_temperature = droplet.m_temperature;
1407 s_backPtr->m_partList[childId].m_heatFlux = 0.0;
1408 s_backPtr->m_partList[childId].m_dM = 0.0;
1414 for(
MInt n = 0; n < nDim; n++) {
1415 magVelChild +=
POW2(childVelocity[n]);
1417 magVelChild = sqrt(magVelChild);
1418 if(fabs(magVel - magVelChild) > MFloatEps) {
1419 cerr <<
"RT kin. energy loss: " << magVel <<
" " << magVelChild << endl;
1421 for(
MInt n = 0; n < nDim; n++) {
1422 mom[n] += sphericalMass(childDiameter, droplet.m_temperature) * dropsPerParcel * childVelocity[n];
1429 for(
MInt n = 0; n < nDim; n++) {
1431 droplet.m_velocity[n] * sphericalMass(droplet.m_diameter, droplet.m_temperature) * droplet.m_noParticles;
1433 if(fabs(oldMom[n] - mom[n]) > MFloatEps) {
1434 cerr <<
"RT momentum change: " << oldMom[n] <<
" " << mom[n] << endl;
1439 const MFloat newMass = sphericalMass(droplet.m_diameter, droplet.m_temperature) * droplet.m_noParticles;
1441 for(
MInt n = 0; n < noNewParcels; n++) {
1442 const MInt id = s_backPtr->a_noParticles() - 1 - n;
1443 childMass += (sphericalMass(s_backPtr->m_partList[
id].m_diameter, s_backPtr->m_partList[
id].m_temperature)
1444 * s_backPtr->m_partList[
id].m_noParticles);
1446 const MFloat massLoss_RT = oldMass - newMass - childMass;
1447 if(fabs(massLoss_RT) > MFloatEps) {
1448 cerr <<
"RT Mass loss! Old " << oldMass <<
" parent " << newMass <<
" childs " << childMass <<
" diff "
1449 << massLoss_RT << endl;
1459 if(m_KHsecBreakUp && droplet.m_diameter > KH_diameter) {
1460 const MFloat KH_mass = sphericalMass(KH_diameter, droplet.m_temperature);
1461 if(KH_mass < MFloatEps)
continue;
1467 droplet.m_shedDiam -= (droplet.m_shedDiam - KH_diameter) / KH_time * dt;
1472 if(droplet.m_shedDiam < KH_diameter) {
1473 MInt noDroplets = floor(oldMass / KH_mass);
1474 if(noDroplets > 0) {
1475 const MFloat diameter = pow(6.0 / PI, 1.0 / 3.0)
1476 * pow(oldMass / (noDroplets * material().density(droplet.m_temperature)), 1.0 / 3.0);
1480 droplet.m_diameter = diameter;
1481 droplet.m_shedDiam = diameter;
1482 droplet.m_noParticles = noDroplets;
1485 const MFloat newMass = sphericalMass(droplet.m_diameter, droplet.m_temperature) * droplet.m_noParticles;
1486 if(fabs(newMass - oldMass) > MFloatEps) {
1487 cerr <<
"Missing mass at KH-breakup 1 is " << oldMass <<
" " << newMass << endl;
1494 const MFloat sheddedMass =
1495 oldMass - sphericalMass(droplet.m_shedDiam, droplet.m_temperature) * droplet.m_noParticles;
1497 if(m_massLimit > sheddedMass / oldMass) {
1501 MInt noSheddedOffDrops = floor(sheddedMass / KH_mass);
1502 if(noSheddedOffDrops <= 0)
continue;
1503 const MFloat parentMass = oldMass - KH_mass * noSheddedOffDrops;
1504 ASSERT(parentMass > 0,
"");
1506 const MFloat parentDiam =
1507 pow(6.0 / PI * parentMass / (droplet.m_noParticles * material().density(droplet.m_temperature)), 1.0 / 3.0);
1513 initPRNG(
globalTimeStep, s_backPtr->c_globalId(droplet.m_cellId));
1516 droplet.m_diameter = parentDiam;
1517 droplet.m_shedDiam = parentDiam;
1520 MFloat angle = m_sprayAngle[0];
1521 if(m_sprayAngleKH > 1 && m_sprayAngleKH < 2.9) {
1522 angle = m_sprayAngle[0];
1523 if(liquidLength < breakupLength) {
1524 angle = m_sprayConeAngle;
1526 }
else if(m_sprayAngleKH > 0 && m_sprayAngleKH < 1) {
1527 angle = 2 * atan(m_sprayAngleKH * KH_growthR * KH_waveL / m_injectionSpeed);
1528 angle = min(m_sprayAngle[0], angle);
1529 }
else if(m_sprayAngleKH > 2.9) {
1530 if(liquidLength < breakupLength) {
1531 angle = m_sprayAngle[0];
1533 angle = m_sprayAngle[1];
1537 if(angle < 0 || angle > 90) {
1538 cerr <<
"Large KH-spray angle : " << angle << endl;
1541 array<MFloat, nDim> childVelocity{};
1543 randomSecBreakUp(), m_injectionDistSigmaCoeff);
1546 const MInt childId = s_backPtr->addParticle(droplet.m_cellId, KH_diameter, droplet.m_densityRatio, 0, 3,
1547 &childVelocity[0], &droplet.m_position[0], noSheddedOffDrops);
1550 s_backPtr->m_partList[childId].m_temperature = droplet.m_temperature;
1551 s_backPtr->m_partList[childId].m_heatFlux = 0.0;
1552 s_backPtr->m_partList[childId].m_dM = 0.0;
1560 if(m_secBUDisplayMessage) {
1561 cerr <<
"#################################################" << endl;
1562 cerr <<
"Break-Up type: KH" << endl;
1563 cerr <<
"parent diameter: " << droplet.m_diameter <<
" " << parentDiam << endl;
1564 cerr <<
"KH size: " << KH_diameter << endl;
1565 cerr <<
"KH time: " << KH_time << endl;
1566 cerr <<
"Time: " << m_timeSinceSOI << endl;
1567 cerr <<
"old #droplets: " << droplet.m_noParticles << endl;
1568 cerr <<
"# new droplets: " << noSheddedOffDrops << endl;
1569 cerr <<
"Weber number: " << We_l << endl;
1570 cerr <<
"#################################################" << endl;
1591 vector<MFloat> mom(nDim);
1592 for(
MInt n = 0; n < nDim; n++) {
1594 droplet.m_velocity[n] * sphericalMass(droplet.m_diameter, droplet.m_temperature) * droplet.m_noParticles
1595 + sphericalMass(KH_diameter, droplet.m_temperature) * noSheddedOffDrops * childVelocity[n];
1596 if(fabs(oldMom[n] - mom[n]) > MFloatEps) {
1597 cerr <<
"KH momentum change: " << oldMom[n] <<
" " << mom[n] <<
" "
1598 << droplet.m_velocity[n] * sphericalMass(droplet.m_diameter, droplet.m_temperature)
1599 * droplet.m_noParticles
1606 for(
MInt n = 0; n < nDim; n++) {
1607 magVelChild +=
POW2(childVelocity[n]);
1609 magVelChild = sqrt(magVelChild);
1610 if(fabs(magVel - magVelChild) > MFloatEps) {
1611 cerr <<
"KH kin. energy loss: " << magVel <<
" " << magVelChild << endl;
1613 MFloat magVelNew = droplet.magVel();
1614 if(fabs(magVel - magVelNew) > MFloatEps) {
1615 cerr <<
"KH kin. energy loss: " << magVel <<
" " << magVelNew << endl;
1618 const MFloat newMass = sphericalMass(droplet.m_diameter, droplet.m_temperature) * droplet.m_noParticles;
1619 const MInt childId = s_backPtr->a_noParticles() - 1;
1621 sphericalMass(s_backPtr->m_partList[childId].m_diameter, s_backPtr->m_partList[childId].m_temperature)
1622 * s_backPtr->m_partList[childId].m_noParticles;
1623 const MFloat massLoss_KH = oldMass - newMass - childMass;
1624 if(fabs(massLoss_KH) > MFloatEps) {
1625 cerr <<
"KH Mass loss! Old " << oldMass <<
" parent " << newMass <<
" childs " << childMass <<
" diff "
1626 << massLoss_KH << endl;
1659 const MInt numInjections = m_injectionRateList.dim0();
1660 if(numInjections == 0)
return;
1662 if(m_injectionTimings[0] > m_timeSinceSOI) {
1663 m_currentInjectionRate = (m_timeSinceSOI + 1) / (m_injectionTimings[0] + 1) * m_injectionRateList[0];
1665 for(
MInt i = 0; i <= numInjections; i++) {
1666 if(m_injectionTimings[i] > m_timeSinceSOI) {
1668 (m_timeSinceSOI - m_injectionTimings[i - 1]) / (m_injectionTimings[i] - m_injectionTimings[i - 1]);
1669 m_currentInjectionRate = m_injectionRateList[i - 1] + t * (m_injectionRateList[i] - m_injectionRateList[i - 1]);
1697 const MFloat particleDiameter,
const MFloat particleDensityRatio,
1698 const MFloat particleVelocity,
const MFloat sprayConeAngle,
1700 std::function<
MFloat(
MInt)> injectorDiameter,
1702 const MInt parcelSize,
const MInt noInjectionHoles) {
1703 std::array<MFloat, nDim> spawnParticlesInitVelo;
1704 array<MFloat, 3> injectionLocation{};
1705 array<MFloat, 3> holeLocation{};
1707 if(spawnParticlesCount > 0) {
1711 if(particleDiameter < s_backPtr->m_sizeLimit) {
1712 cerr <<
"Inj. small particle " << particleDiameter <<
" " << s_backPtr->m_sizeLimit << endl;
1715 for(
MInt i = 0; i < spawnParticlesCount; i++) {
1716 MInt randomShift = 1;
1719 m_PRNGPrimBreakUpCount +=
1720 randomVectorInCone(&spawnParticlesInitVelo[0], &m_injectorDir[0], particleVelocity, sprayConeAngle,
1721 m_partEmittDist, randomPrimBreakUp(0), m_injectionDistSigmaCoeff, 0);
1723 if(coneDiameter > numeric_limits<MFloat>::epsilon() && !hull) {
1725 randomPointInCircle(&holeLocation[0], &m_injectorDir[0], coneDiameter, randomPrimBreakUp(2));
1727 if(noInjectionHoles > 1) {
1729 pointOnCircle(&injectionLocation[0], &m_injectorDir[0], injectorDiameter(i), holePositionAngle(i));
1732 m_PRNGPrimBreakUpCount +=
1733 randomVectorInCone(&spawnParticlesInitVelo[0], &m_injectorDir[0], particleVelocity, sprayConeAngle,
1734 m_partEmittDist, randomPrimBreakUp(0), m_injectionDistSigmaCoeff, holeNozzleAngle(i));
1738 static constexpr MBool randomDist =
false;
1740 randomPointOnCircle(&injectionLocation[0], &m_injectorDir[0], coneDiameter, randomPrimBreakUp(1),
1741 spawnParticlesCount, i);
1744 const static MFloat angleBetween = 360.0 / spawnParticlesCount;
1745 const MFloat phi = (i * angleBetween + ((m_injStep - 1) %
static_cast<MInt>(angleBetween))) / 180.0 * M_PI;
1746 pointOnCircle(&injectionLocation[0], &m_injectorDir[0], coneDiameter, phi);
1752 if(fabs(holeNozzleAngle(i)) > numeric_limits<MFloat>::epsilon()) {
1753 const MFloat holeAngleRad = -holeNozzleAngle(i) / 180 * M_PI;
1756 array<MFloat, nDim> crossP{};
1757 array<MFloat, nDim> revInjLoc{};
1758 for(
MInt j = 0; j < nDim; j++) {
1759 revInjLoc[j] = -injectionLocation[j] / absDist;
1766 for(
MInt j = 0; j < nDim; j++) {
1767 injLocBasis(j, 0) = m_injectorDir[j];
1768 injLocBasis(j, 1) = revInjLoc[j];
1769 injLocBasis(j, 2) = crossP[j];
1774 for(
MInt j = 0; j < nDim; j++) {
1775 for(
MInt k = 0; k < nDim; k++) {
1776 inverse(j, k) = injLocBasis(j, k);
1787 R(0, 0) = cos(holeAngleRad);
1788 R(1, 0) = sin(holeAngleRad);
1789 R(0, 1) = -sin(holeAngleRad);
1790 R(1, 1) = cos(holeAngleRad);
1796 for(
MInt j = 0; j < nDim; j++) {
1797 tempV[j] = spawnParticlesInitVelo[j];
1802 for(
MInt j = 0; j < nDim; j++) {
1803 spawnParticlesInitVelo[j] = 0;
1804 for(
MInt k = 0; k < nDim; k++) {
1805 spawnParticlesInitVelo[j] += result(j, k) * tempV[k];
1812 if(
scalarProduct(&spawnParticlesInitVelo[0], &m_injectorDir[0], nDim) < 0.0) {
1813 for(
MInt j = 0; j < nDim; j++) {
1814 spawnParticlesInitVelo[j] *= -1;
1817 for(
MInt j = 0; j < nDim; j++) {
1818 spawnParticlesInitVelo[j] *= particleVelocity;
1822 for(
MInt j = 0; j < nDim; j++) {
1823 injectionLocation[j] += spawnCoord[j] + holeLocation[j];
1826 s_backPtr->addParticle(s_backPtr->m_spawnCellId, particleDiameter, particleDensityRatio, randomShift, 2,
1827 &spawnParticlesInitVelo[0], &injectionLocation[0], parcelSize);
1838 if(domainId() != 0)
return;
1840 const MInt noDraws = 1000000;
1841 MFloat minValue = m_injectorNozzleDiameter / m_RosinRammlerMin;
1842 MFloat maxValue = m_injectorNozzleDiameter / m_RosinRammlerMax;
1843 MFloat meanValue = m_injectorNozzleDiameter / m_RosinRammlerMean;
1845 std::mt19937_64 randomNumberGenerator;
1846 randomNumberGenerator.seed(m_spraySeed);
1848 const MInt noBinns = 1000;
1849 MFloat binWidth = (maxValue - minValue) / noBinns;
1853 for(
MInt i = 0; i < noDraws; i++) {
1854 const MFloat drawValue =
rosinRammler(minValue, meanValue, maxValue, m_RosinRammlerSpread, randomNumberGenerator);
1855 const MInt binId = std::floor((drawValue - minValue) / binWidth);
1859 m_log <<
"IDSD Binning: " << endl;
1860 m_log <<
"Rosin-Rammler distribution: " << endl;
1861 m_log <<
"Min-value : " << minValue << endl;
1862 m_log <<
"Max-value : " << maxValue << endl;
1863 m_log <<
"Mean-value: " << meanValue << endl;
1864 m_log <<
"Spread : " << m_RosinRammlerSpread << endl;
1865 m_log <<
"Number-draws: " << noDraws << endl;
1867 for(
MInt i = 0; i < noBinns; i++) {
1868 m_log << minValue + i * binWidth <<
" " << binValue[i] << endl;
1872 maxValue = 5 * meanValue;
1873 binWidth = (maxValue - minValue) / noBinns;
1877 for(
MInt i = 0; i < noDraws; i++) {
1879 const MInt binId = std::floor((drawValue - minValue) / binWidth);
1883 m_log <<
"Nukiyama-Tanasawa distribution: " << endl;
1884 m_log <<
"Mean-value: " << meanValue << endl;
1885 for(
MInt i = 0; i < noBinns; i++) {
1886 m_log << minValue + i * binWidth <<
" " << binValue[i] << endl;
void mAlloc(T *&a, const MLong N, const MString &objectName, MString function)
allocates memory for one-dimensional array 'a' of size N
static MInt propertyLength(const MString &name, MInt solverId=m_noSolvers)
Returns the number of elements of a property.
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
std::vector< LPTSpherical< nDim > > m_partList
This class is a ScratchSpace.
void fill(T val)
fill the scratch with a given value
void injectParticles(const MInt spawnParticlesCount, const MFloat *spawnCoord, const MFloat particleDiameter, const MFloat particleDensityRatio, const MFloat particleVelocity, const MFloat sprayConeAngle, const MFloat coneDiameter, std::function< MFloat(MInt)> holePositionAngle, std::function< MFloat(MInt)> injectorDiameter, std::function< MFloat(MInt)> holeNozzleAngle, const MBool hull, const MInt parcelSize, const MInt noInjectionHoles)
Create new random particles using the provided options.
void logDistributionFunction()
create log of distribution functions for plots and debugging
void injection(const MFloat dt)
Perform a spray injection.
void primaryBreakUp(const MFloat dt)
Inject new particles from the injector exit.
void updateInjectionRate()
Update injection rate.
void init(LPT< nDim > *_particle)
Spray model initialisation.
void readSprayProperties()
void secondaryBreakUp(const MFloat dt)
Atomization of particles.
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
void mTerm(const MInt errorCode, const MString &location, const MString &message)
constexpr Real POW3(const Real x)
constexpr Real POW2(const Real x)
constexpr T mMax(const T &x, const T &y)
std::basic_string< char > MString
MFloat NTDistribution(const MFloat x_mean, std::mt19937_64 &PRNG)
MFloat scalarProduct(const MFloat *a, const MFloat *b, const MInt length)
void randomPointInCircle(MFloat *vec, const MFloat *normalDirection, const MFloat diameter, std::mt19937_64 &PRNG)
void randomPointOnCircle(MFloat *vec, const MFloat *normalDirection, const MFloat diameter, std::mt19937_64 &PRNG, const MInt circleSplit=1, const MInt splitNo=0)
MInt randomVectorInCone(MFloat *vec, const MFloat *coneAxis, const MFloat length, const MFloat openingAngle, const MInt dist, std::mt19937_64 &PRNG, const MFloat distCoeff=0.0, const MFloat nozzleAngle=0.0)
Generate a random vector in a cone defined by its opening angle.
MFloat rosinRammler(const MFloat min, const MFloat mean, const MFloat max, const MFloat spread, std::mt19937_64 &PRNG)
void pointOnCircle(MFloat *vec, const MFloat *normalDirection, const MFloat diameter, MFloat phi)
void cross(const T *const u, const T *const v, T *const c)
MFloat distance(const MFloat *a, const MFloat *b)
void invert(MFloat *A, const MInt m, const MInt n)
void normalize(std::array< T, N > &u)
MFloat crankAngle(const MFloat time, const MFloat Strouhal, const MFloat offset, const MInt mode)
help-function for engine calculations which returns the crank-angle for a given time,...
MFloat norm(const std::array< T, N > &u)
void multiplyMatricesSq(MFloatScratchSpace &m1, MFloatScratchSpace &m2, MFloatScratchSpace &result, MInt dim)
Namespace for auxiliary functions/classes.