MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
lptspray.cpp
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
7#include "lptspray.h"
8#include <cmath>
9#include "IO/tomlutils.h"
10#include "lpt.h"
11#include "lptlib.h"
12
13using namespace maia::lpt;
14using namespace std;
15using namespace maia;
16
17template <MInt nDim>
19template <MInt nDim>
21
28template <MInt nDim>
30 m_log << " Initialising spray model" << endl;
31
32 s_backPtr = _particle;
33 s_lengthFactor = s_backPtr->m_partList[0].s_lengthFactor;
34
35 m_activePrimaryBUp = s_backPtr->m_activePrimaryBUp;
36 m_activeSecondaryBUp = s_backPtr->m_activeSecondaryBUp;
37
38 if(!s_backPtr->m_restart) {
39 s_backPtr->m_particleResiduum = 0.0;
40 }
41
42 // allocate spray data
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_);
46
47 readSprayProperties();
48
49 logDistributionFunction();
50}
51
52
58template <MInt nDim>
60 // debug parameter to validate injection/conservation
61 m_injStopTimeStep = -1;
62 m_injStopTimeStep = Context::getSolverProperty<MInt>("injStopTimeStep", solverId(), AT_, &m_injStopTimeStep);
63
64 m_injStartTimeStep = -1;
65 m_injStartTimeStep = Context::getSolverProperty<MInt>("injStartTimeStep", solverId(), AT_, &m_injStartTimeStep);
66
67 if(m_activePrimaryBUp) {
68 m_needleLiftTime[0] = 0.01;
69 m_needleLiftTime[1] = 0.02;
78 if(Context::propertyExists("injectorNeedleLiftTime", solverId())) {
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;
82 }
83
91 m_injDuration = 0.00055;
92 m_injDuration = Context::getSolverProperty<MFloat>("injectorInjectionTime", solverId(), AT_, &m_injDuration);
93
102 m_injectionCrankAngle = -1;
103 m_injectionCrankAngle =
104 Context::getSolverProperty<MFloat>("injectionCrankAngle", solverId(), AT_, &m_injectionCrankAngle);
105
106 // NOTE: the Strouhal number must be different from the FV Strouhal-number
107 // if a different non-dimensionalisation is used in the LPT-solver!
108 m_Strouhal = Context::getSolverProperty<MFloat>("Strouhal", solverId(), AT_, &m_Strouhal);
109
110 m_initialCad = Context::getSolverProperty<MFloat>("initialCrankAngle", solverId(), AT_, &m_initialCad);
111
119 m_primRosinRammler = true;
120 m_primRosinRammler =
121 Context::getSolverProperty<MBool>("sprayPrimRosinRammler", solverId(), AT_, &m_primRosinRammler);
122
132 MString injectorType = "FULLCONE";
133 injectorType = Context::getSolverProperty<MString>("injectorType", solverId(), AT_, &injectorType);
134 m_injectorType = string2enum(injectorType);
135
136 // injector categories:
137 m_multiHoleInjector = false;
138 m_singleHoleInjector = false;
139 m_hollowConeInjector = false;
140
141 switch(m_injectorType) {
142 case MULTIHOLE:
143 case MULTIHOLE_OPT:
144 case MULTIHOLE_TME: {
145 m_multiHoleInjector = true;
146 break;
147 }
148 case FULLCONE: {
149 m_singleHoleInjector = true;
150 break;
151 }
152 case HOLLOWCONE: {
153 m_hollowConeInjector = true;
154 break;
155 }
156 default: {
157 mTerm(1, AT_, "Unknown injector-type!");
158 }
159 }
160
168 m_injectorNozzleDiameter = Context::getSolverProperty<MFloat>("injectorNozzleDiameter", solverId(), AT_);
169
170 // read injector type specific properties
171 if(m_multiHoleInjector) {
181 if(Context::propertyLength("injectorDiameter") == 1) {
182 mAlloc(m_injectorDiameter, 1, "injectorDiameter", AT_);
183 m_injectorDiameter[0] = Context::getSolverProperty<MFloat>("injectorDiameter", solverId(), AT_);
184 } else {
185 const MInt length = Context::propertyLength("injectorDiameter");
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]);
191 }
192 }
193
202 if(!Context::propertyExists("injectorOrificeAngle") || Context::propertyLength("injectorOrificeAngle") == 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]);
207 } else {
208 const MInt length = Context::propertyLength("injectorOrificeAngle");
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);
212 }
213 }
222 if(Context::propertyExists("orificePositionAngle")) {
223 const MInt length = Context::propertyLength("orificePositionAngle");
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);
227 }
228 }
229 }
230
238 m_injectorOrficeSize = Context::getSolverProperty<MFloat>("injectorOrficeSize", solverId(), AT_);
239
240
241 if(m_primRosinRammler && m_singleHoleInjector) {
242 // determine diameter factors for the rosin rammler initial droplet size distriution (IDSD)
243 // default values are fraction of the nominal injector nozzle diameter for spray-A as in:
244 // LARGE EDDY SIMULATION OF HIGH-VELOCITY FUEL SPRAYS:
245 // STUDYING MESH RESOLUTION AND BREAKUP MODEL EFFECTS FOR SPRAY A
246 // A. Wehrfritz, V. Vuorinen, O. Kaario, & M. Larmi
247 // Atomization and Sprays, 23 (5): 419–442 (2013)
248 m_RosinRammlerMin = 90;
249 m_RosinRammlerMean = 15;
250 m_RosinRammlerMax = 5;
251 m_RosinRammlerSpread = 3;
252
253 } else if(m_primRosinRammler && m_hollowConeInjector) {
254 // determine diameter factors for the rosin rammler initial droplet size distriution (IDSD)
255 // default values are fraction of the liquid sheet length for hollow-cone injectors
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;
265 }
266
267 // read specified values from the property file
268 if(m_primRosinRammler) {
276 m_RosinRammlerSpread =
277 Context::getSolverProperty<MFloat>("RosinRammlerSpread", solverId(), AT_, &m_RosinRammlerSpread);
278
286 m_RosinRammlerMin = Context::getSolverProperty<MFloat>("RosinRammlerMin", solverId(), AT_, &m_RosinRammlerMin);
287
295 m_RosinRammlerMax = Context::getSolverProperty<MFloat>("RosinRammlerMax", solverId(), AT_, &m_RosinRammlerMax);
296
304 m_RosinRammlerMean = Context::getSolverProperty<MFloat>("RosinRammlerMean", solverId(), AT_, &m_RosinRammlerMean);
305 }
306
314 m_injectionDistSigmaCoeff =
315 Context::getSolverProperty<MFloat>("spawnDistSigmaCoeff", solverId(), AT_, &m_injectionDistSigmaCoeff);
316
324 MString partDist = "PART_EMITT_DIST_NONE";
325 m_partEmittDist = string2enum(Context::getSolverProperty<MString>("partEmittDist", solverId(), AT_, &partDist));
326
327 MFloat diameter = m_injectorNozzleDiameter;
328 if(m_multiHoleInjector) {
329 diameter = m_injectorDiameter[0];
330 } else if(m_singleHoleInjector) {
331 diameter = m_injectorOrficeSize;
332 }
333
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;
338 }
339 }
340
348 m_injectionSpeed = Context::getSolverProperty<MFloat>("injectionSpeed", solverId(), AT_, &m_injectionSpeed);
349
350 ASSERT(m_injectionSpeed > 0, "");
351
352 if(Context::propertyExists("sprayInjectionRate", solverId())) {
353 const MInt numInjections = Context::propertyLength("sprayInjectionRate", solverId());
354
355 if(numInjections == 1) {
356 m_currentInjectionRate = Context::getSolverProperty<MFloat>("sprayInjectionRate", solverId(), AT_);
357
358 } else {
359 m_injectionRateList.resize(numInjections);
360 m_injectionTimings.resize(numInjections);
361
362 for(MInt i = 0; i < numInjections; i++) {
370 m_injectionRateList[i] = Context::getSolverProperty<MFloat>("sprayInjectionRate", solverId(), AT_, i);
371
379 m_injectionTimings[i] = Context::getSolverProperty<MFloat>("sprayInjectionTiming", solverId(), AT_, i);
380 }
381 }
382 } else {
383 ASSERT(!m_activePrimaryBUp, "");
384 }
385
394 for(MInt i = 0; i < nDim; i++) {
395 m_injectorDir[i] = 0;
396 if(i == nDim - 1) m_injectorDir[i] = 1.0;
397 }
398
399 if(Context::propertyExists("injectorDir", solverId())) {
400 if(Context::propertyLength("injectorDir", solverId()) != nDim) {
401 TERMM(1, "Need to give a Coordinate for every dimension");
402 }
403
404 for(MInt i = 0; i < nDim; i++) {
405 m_injectorDir[i] = Context::getSolverProperty<MFloat>("injectorDir", solverId(), AT_, i);
406 }
407 }
408
409 if(m_hollowConeInjector) {
417 m_maxNoPrimParcels = 1;
418 m_maxNoPrimParcels =
419 Context::getSolverProperty<MInt>("sprayPrimaryMaxPrimParcels", solverId(), AT_, &m_maxNoPrimParcels);
420
428 m_primParcelsPerInj = 1000000;
429 m_primParcelsPerInj =
430 Context::getSolverProperty<MInt>("sprayPrimaryParcelsPerInj", solverId(), AT_, &m_primParcelsPerInj);
431
439 m_angularGap = Context::getSolverProperty<MFloat>("injectorAngularGap", solverId(), AT_, &m_angularGap);
440 }
441
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.");
456 }
457 }
458
462
463
464 if(Context::propertyLength("sprayAngle", solverId()) == 1) {
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]);
467 } else {
468 m_sprayAngle[0] = Context::getSolverProperty<MFloat>("sprayAngle", solverId(), AT_, 0);
469 m_sprayAngle[1] = Context::getSolverProperty<MFloat>("sprayAngle", solverId(), AT_, 1);
470 }
471
472
473 m_sprayConeAngle = Context::getSolverProperty<MFloat>("sprayConeAngle", solverId(), AT_, &m_sprayConeAngle);
474
482 m_spraySeed = 5489U;
483 if(Context::propertyExists("spawnSeed", solverId())) {
484 m_spraySeed = (MLong)Context::getSolverProperty<MInt>("spawnSeed", solverId(), AT_);
485 }
486
487 if(m_activePrimaryBUp) m_PRNGPrimBreakUp.seed(m_spraySeed);
488 if(m_activeSecondaryBUp) m_PRNGSecBU.seed(m_spraySeed);
489
491 // secondary breakup properties
493 if(m_activeSecondaryBUp) {
501 m_secBUDisplayMessage = false;
502 m_secBUDisplayMessage =
503 Context::getSolverProperty<MBool>("spraysecBUPOutput", solverId(), AT_, &m_secBUDisplayMessage);
504
505 m_RTsecBreakUp = true;
506 m_RTsecBreakUp = Context::getSolverProperty<MBool>("RayleighTaylorBreakUp", solverId(), AT_, &m_RTsecBreakUp);
507
508 m_KHsecBreakUp = true;
509 m_KHsecBreakUp = Context::getSolverProperty<MBool>("KelvinHelmholtzBreakUp", solverId(), AT_, &m_KHsecBreakUp);
510
511 m_RTsecBreakUpLength = true;
512 m_RTsecBreakUpLength =
513 Context::getSolverProperty<MBool>("RayleighTaylorBreakUpLength", solverId(), AT_, &m_RTsecBreakUpLength);
514
515
527 m_predictivePRNG = true;
528 m_predictivePRNG = Context::getSolverProperty<MBool>("predictablePRNGSecBU", solverId(), AT_, &m_predictivePRNG);
529
530 // initialise PRNG
531 m_PRNGSecBU.discard(globalTimeStep);
532
544 m_sprayAngleKH = 2;
545 m_sprayAngleKH = Context::getSolverProperty<MFloat>("sprayAngleKH", solverId(), AT_, &m_sprayAngleKH);
546
547 // default values for secondary break-up are based on:
548 // LARGE EDDY SIMULATION OF HIGH-VELOCITY FUEL SPRAYS:
549 // STUDYING MESH RESOLUTION AND BREAKUP MODEL EFFECTS FOR SPRAY A
550 // A. Wehrfritz, V. Vuorinen, O. Kaario, & M. Larmi
551 // Atomization and Sprays, 23 (5): 419–442 (2013)
552
553
561 m_weberLimitRT = Context::getSolverProperty<MFloat>("sprayWeberLimitRT", solverId(), AT_, &m_weberLimitRT);
562
570 m_B0 = Context::getSolverProperty<MFloat>("sprayB0", solverId(), AT_, &m_B0);
571
579 // valid range of B1 is 10-60 according to Reitz1999
580 m_B1 = Context::getSolverProperty<MFloat>("sprayB1", solverId(), AT_, &m_B1);
581
589 // other values used by Sven Berger:
590 // 0.15 AIA-Nozzle MPa 10
591 // 0.20 AIA-Nozzle MPa 20
592 m_CRT = Context::getSolverProperty<MFloat>("sprayCRT", solverId(), AT_, &m_CRT);
593
601 m_CT = Context::getSolverProperty<MFloat>("sprayCT", solverId(), AT_, &m_CT);
602
610 m_weberLimitKH = Context::getSolverProperty<MFloat>("sprayWeberLimitKH", solverId(), AT_, &m_weberLimitKH);
611
620 // In literature sometimes denoted as B2
621 m_massLimit = Context::getSolverProperty<MFloat>("sprayMassLimit", solverId(), AT_, &m_massLimit);
622
630 m_maxRTChildParcels = Context::getSolverProperty<MInt>("sprayMaxRTChilds", solverId(), AT_, &m_maxRTChildParcels);
631
642 // NOTE: version 1 is following: "Implementation and validation of a Lagrangian spray model using
643 // experimental data of the EVN Spray G injector" Horacia J. Aguerre, Norberto M. Nigro
644 // Computers and Fluids 190(2019) 30-48
645 m_RTDiameterMode = Context::getSolverProperty<MInt>("sprayRTDiameterMode", solverId(), AT_, &m_RTDiameterMode);
646
655 // default based on Reitz et al.
656 //"Modeling Spray Atomization with the Kelvin-Helmholtz/Rayleigh-Taylor
657 // Hybrid Model" Atomization and Sprays 1999
658 m_Cbu = m_B1 / 2;
659 m_Cbu = Context::getSolverProperty<MFloat>("sprayCbu", solverId(), AT_, &m_Cbu);
660
661 MFloat breakUpLength = m_Cbu * m_injectorNozzleDiameter * sqrt(material().ambientDensityRatio());
662
663 cerr0 << "Liquid break-up length for KH-RT model is " << breakUpLength << endl;
664 }
665}
666
667
671template <MInt nDim>
673 // injection has not started yet, this also means that zero particles should be present!
674 if(m_injectionCrankAngle > maia::math::crankAngle(s_backPtr->m_time, m_Strouhal, m_initialCad, 0)) {
675 ASSERT(s_backPtr->a_noParticles() == 0, "");
676 if(domainId() == 0) {
677 cerr << "Before injection" << m_injectionCrankAngle << " "
678 << maia::math::crankAngle(s_backPtr->m_time, m_Strouhal, m_initialCad, 0) << endl;
679 }
680 for(MInt i = 0; i < m_injDataSize; i++) {
681 m_injData[i] = 0;
682 }
683 return;
684 }
685
686 m_timeSinceSOI += dt;
687
688
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;
693
694 // injection has already ended
695 if(endOfInjection) {
696 cerr0 << "After injection " << m_timeSinceSOI << " " << m_injDuration << endl;
697
698 for(MInt i = 0; i < m_injDataSize; i++) {
699 m_injData[i] = 0;
700 }
701 if(s_backPtr->m_spawnCellId > -1) {
702 m_injData[0] = m_timeSinceSOI;
703 }
704 return;
705 }
706
707
708 // injection location is on another domain
709 if(s_backPtr->m_spawnCellId < 0) {
710 ASSERT(s_backPtr->noDomains() > 1, "");
711
712 // receive injected particles instead
713 if(m_broadcastInjected) {
714 s_backPtr->recvInjected();
715 }
716 for(MInt i = 0; i < m_injDataSize; i++) {
717 m_injData[i] = 0;
718 }
719
720 } else {
721 // injection on this domain
722 ASSERT(domainId() == s_backPtr->m_spawnDomainId, "");
723
724 // update injection rate
725 updateInjectionRate();
726
727 // current number of particles
728 const MInt prevNo = s_backPtr->a_noParticles();
729
730 // inject new particles
731 primaryBreakUp(dt);
732
733 if(m_broadcastInjected) {
734 s_backPtr->broadcastInjected(prevNo);
735 }
736 }
737}
738
743template <MInt nDim>
745 TRACE();
746
747 ASSERT(domainId() == s_backPtr->m_spawnDomainId, "");
748 ASSERT(s_backPtr->m_spawnCellId > -1, "");
749
750 // invalidate injection data:
751 MInt noInjParticles = -1;
752 MInt noInjDroplets = -1;
753 MFloat injMass = NAN;
754 MFloat injParticleDiameter = NAN;
755 static MFloat lastD = -1;
756
757 // 1) ramp-up and other general injection parameters:
758
759 MFloat injectionProgress = m_timeSinceSOI / m_injDuration;
760
761
762 // the rampFactor is used to control/ramp the massflowrate and velocity
763 // during the injector opening and closing
764 // which linearly increases/decreases to the full injection velocity massflowrate
765 MFloat rampFactor = 1.0;
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];
770 }
771
772 // theoretically injected mass during this timestep
773 // current injection rate + left over mass from previous timestep
774 injMass = m_currentInjectionRate * rampFactor * dt + s_backPtr->m_particleResiduum;
775
776 // considering ramp-factor for parcel count
777 MInt noDroplets = mMax(1, (MInt)(m_primBrkUpParcelSize * rampFactor));
778
779 // 2) injector dependant injection parameters
780 // and injection of new particles
781 if(m_singleHoleInjector) { //(Spray A)
782
783 if(m_primRosinRammler) {
784 // droplet diameter distribution function
785
786 const MFloat diameterMin = m_injectorNozzleDiameter / m_RosinRammlerMin;
787 const MFloat diameterMean = m_injectorNozzleDiameter / m_RosinRammlerMean;
788 const MFloat diameterMax = m_injectorNozzleDiameter / m_RosinRammlerMax;
789
790 // mass of the smallest possible droplet
791 const MFloat minMass = sphericalMass(diameterMin, material().T());
792
793 noInjParticles = 0;
794 noInjDroplets = 0;
795
796 MFloat remainingMass = injMass;
797 injMass = 0;
798
799 while(remainingMass > minMass) {
800 if(lastD > 0) {
801 injParticleDiameter = lastD;
802 } else {
803 injParticleDiameter =
804 rosinRammler(diameterMin, diameterMean, diameterMax, m_RosinRammlerSpread, randomPrimBreakUp(1));
805 }
806
807 const MFloat particleMass = noDroplets * sphericalMass(injParticleDiameter, material().T());
808 const MFloat tempMass = remainingMass - particleMass;
809
810 // limit mass to be positive
811 if(tempMass < 0) {
812 lastD = injParticleDiameter;
813 s_backPtr->m_particleResiduum = remainingMass;
814 break;
815 }
816
817 // the initial cone diameter needs to be reduced by half of the droplet diameter
818 // since the droplet is supposed to be fully within the cone
819 const MFloat coneDiameter = m_injectorOrficeSize - 0.5 * injParticleDiameter;
820 lastD = -1;
821 injectParticles(1, s_backPtr->m_spawnCoord, injParticleDiameter, material().densityRatio(),
822 m_injectionSpeed * rampFactor, m_sprayConeAngle, coneDiameter, false, 0.0, noDroplets, 1);
823
824 noInjParticles++;
825 noInjDroplets += noDroplets;
826 injMass += particleMass;
827 remainingMass -= particleMass;
828 }
829
830 } else {
831 // blob-injection with all particles at the injector center
832
833 injParticleDiameter = m_injectorNozzleDiameter;
834 noInjDroplets = 1;
835
836 const MFloat initParticleMass = sphericalMass(injParticleDiameter, material().T());
837
838
839 MFloat noNewParticles = (m_currentInjectionRate * dt) / initParticleMass + s_backPtr->m_particleResiduum;
840 noInjParticles = floor(noNewParticles);
841
842 s_backPtr->m_particleResiduum = noNewParticles - noInjParticles;
843
844 injectParticles(noInjParticles, s_backPtr->m_spawnCoord, injParticleDiameter, material().densityRatio(),
845 m_injectionSpeed, m_sprayConeAngle, 0.0, false, 0.0, 1, 1);
846 }
847 } else if(m_hollowConeInjector) {
848 // hollow cone injector, such as the TMFB Injector
849 // This implementation of the hollow cone injector model of:
850 // Spray Modeling for Outwardly-Opening Hollow-Cone Injectors, Sim, Badra et al., SAE, 2016
851
852 // maximum number of parcels to be generated during this time step
853 const MInt parcelTargetNo = m_primParcelsPerInj / m_injDuration * dt;
854
855 // limit the number of parcels to be injected to maxNoPrimParcels
856 noInjParticles = std::min(parcelTargetNo, m_maxNoPrimParcels);
857
858 // injection velocity is scaled by this factor
859 const MFloat scaledInjectionV = rampFactor * m_injectionSpeed;
860
861 // thickness of the fluid flow within the angular gap
862 MFloat initialLiquidSheetThickness = NAN;
863
864 if(m_angularGap > -1) {
865 initialLiquidSheetThickness = max(m_angularGap * rampFactor, m_primMinDropletSize);
866 } else {
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);
873 }
874
875 // initial particle size when considered as string like coneshaped
876 // starting from the injector angular gap
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));
884 }
885
886 // mass of a particle to be added
887 const MFloat initialDropletMass = sphericalMass(injParticleDiameter, material().T());
888
889 const MFloat newDroplets = injMass / initialDropletMass;
890
891 noInjDroplets = floor(newDroplets / noInjParticles);
892
893 if(noInjDroplets > 0) {
894 const MFloat injectorGapDiameter = m_injectorNozzleDiameter - 0.5 * injParticleDiameter;
895 s_backPtr->m_particleResiduum = (newDroplets - noInjParticles * noInjDroplets) * initialDropletMass;
896
897 injectParticles(noInjParticles, s_backPtr->m_spawnCoord, injParticleDiameter, material().densityRatio(),
898 scaledInjectionV, m_sprayConeAngle, injectorGapDiameter * s_lengthFactor, true, m_angularGap,
899 noInjDroplets, 1.0);
900
901 injMass = noInjParticles * noInjDroplets * initialDropletMass;
902
903 } else {
904 injMass = 0;
905 noInjParticles = 0;
906 }
907
908 } else if(m_multiHoleInjector) {
909 // injectors with multiple injection holes
910 // MULTIHOLE : ECN Spray-G
911 // MULTIHOLE_Opt : optimized 7-hole Spray-G
912 // MULTIHOLE_TME : 5-hole Injection of the TME
913
914 MInt noHoles = -1;
915 switch(m_injectorType) {
916 case MULTIHOLE: {
917 noHoles = 8;
918 break;
919 }
920 case MULTIHOLE_OPT: {
921 noHoles = 7;
922 break;
923 }
924 case MULTIHOLE_TME: {
925 noHoles = 5;
926 break;
927 }
928 default:
929 mTerm(1, AT_, "Unknown multi-hole injector-type!");
930 }
931
932 static constexpr MFloat C_a = 0.65;
933 // NOTE: Wehrfritz et. al. propose to use the nominal injector nozzle diameter!
934 // otherwise the area-coefficient of ~0.65 can be used
935
936 const MFloat initialConeAngle = m_sprayConeAngle;
937 // as in:
938 //"Validation of a comprehensive computational fluid dynamics methodology to predict the direct
939 // injection process of gasoline sprays using Spray G experimental data"
940 // Davide Paredi , Tommaso Lucchini , Gianluca D’Errico, Angelo Onorati,
941 // Lyle Pickett and Joshua Lacey
942 // International J of Engine Research 2020, Vol. 21(1) 199–216
943
944 ASSERT(m_primRosinRammler, "");
945
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;
950
951 MFloat remainingMass = injMass;
952 const MFloat minMass = sphericalMass(diameterMin, material().T());
953
954 injMass = 0;
955 noInjParticles = 0;
956 noInjDroplets = 0;
957
958 std::function<MFloat(MInt)> holePosition = [&](MInt id) {
959 const MFloat angleDelta = 2 * M_PI / 8;
960 switch(m_injectorType) {
961 case MULTIHOLE: {
962 return (id * angleDelta);
963 }
964 case MULTIHOLE_OPT: {
965 // 0: x=0, z=-1
966 // 2: x=-1, z=0
967 // 4: x=0, z=1
968 // 6: x=1, z=0
969 if(id > 1) {
970 id = id + 1;
971 }
972 return id * angleDelta;
973 }
974 case MULTIHOLE_TME: {
975 return m_orificePositionAngle[id] / 180 * M_PI;
976 }
977 default:
978 mTerm(1, AT_, "Unknown multi-hole injector-type!");
979 }
980 };
981
982 std::function<MFloat(MInt)> injectorDiameter = [&](MInt id) {
983 switch(m_injectorType) {
984 case MULTIHOLE:
985 case MULTIHOLE_OPT: {
986 return m_injectorDiameter[0];
987 }
988 case MULTIHOLE_TME: {
989 return m_injectorDiameter[id + 1];
990 }
991 default:
992 mTerm(1, AT_, "Unknown multi-hole injector-type!");
993 }
994 };
995
996 // return angle [degree] of the hole orientation
997 std::function<MFloat(MInt)> holeAngle = [&](MInt id) {
998 switch(m_injectorType) {
999 case MULTIHOLE:
1000 case MULTIHOLE_OPT: {
1001 return 37.0;
1002 }
1003 case MULTIHOLE_TME: {
1004 return m_injectorOrificeAngle[id];
1005 }
1006 default:
1007 mTerm(1, AT_, "Unknown multi-hole injector-type!");
1008 }
1009 };
1010
1011 while(remainingMass > minMass) {
1012 if(lastD > 0) {
1013 injParticleDiameter = lastD;
1014 } else {
1015 injParticleDiameter =
1016 rosinRammler(diameterMin, diameterMean, diameterMax, m_RosinRammlerSpread, randomPrimBreakUp(1));
1017 }
1018
1019 const MFloat particleMass = noHoles * noDroplets * sphericalMass(injParticleDiameter, material().T());
1020
1021 const MFloat tempMass = remainingMass - particleMass;
1022
1023 if(tempMass < 0) {
1024 lastD = injParticleDiameter;
1025 s_backPtr->m_particleResiduum = remainingMass;
1026 break;
1027 }
1028
1029 lastD = -1;
1030 injectParticles(noHoles, s_backPtr->m_spawnCoord, injParticleDiameter, material().densityRatio(),
1031 m_injectionSpeed, initialConeAngle, m_injectorOrficeSize, holePosition, injectorDiameter,
1032 holeAngle, false, noDroplets, noHoles);
1033
1034 injMass += particleMass;
1035 noInjParticles += noHoles;
1036 noInjDroplets += noHoles * noDroplets;
1037 remainingMass -= particleMass;
1038 }
1039 }
1040
1041 // 3) store injection data for possible post-processing:
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;
1049
1050 array<MFloat, 3> injMomentum = {};
1051 for(MInt i = 0; i < 3; i++) {
1052 injMomentum[i] = 0;
1053 }
1054 MFloat injEnergy = 0;
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]);
1063 }
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;
1067 }
1068 m_injData[7] = injMomentum[0];
1069 m_injData[8] = injMomentum[1];
1070 m_injData[9] = injMomentum[2];
1071 m_injData[10] = injEnergy;
1072
1073 cerr << globalTimeStep << " " << m_timeSinceSOI << " injMass " << injMass << " injParticles " << noInjParticles
1074 << endl;
1075}
1076
1080template <MInt nDim>
1082 TRACE();
1083
1084 m_noRTsecBreakUp = 0;
1085 m_noKHsecBreakUp = 0;
1086
1087 if(!m_RTsecBreakUp && !m_KHsecBreakUp) return;
1088
1089 // use const to avoid breakup of particles which have just been created
1090 const MInt noPart = s_backPtr->a_noParticles();
1091
1092 // increase size of particle vector if necessary before joining next loop
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());
1095 }
1096
1097 // 0) Compute break-up length:
1098 // (usually used when no initial droplet size distribution (IDS) is available and
1099 // the blob method is used
1100 // Break-up length can be calculated by Levich theory according to Levich, 1962
1101
1102 const MFloat breakupLength =
1103 m_RTsecBreakUpLength ? m_Cbu * m_injectorNozzleDiameter * sqrt(material().ambientDensityRatio()) : 0;
1104
1105 // if(droplet.m_breakUpTime > 0.43839) {
1106 // reset for particles which have been affected significantly by tumble motion!
1107 // breakupLength = 0.0;
1108 //}
1109
1110
1111 ASSERT(m_injectorNozzleDiameter > 0, "ERROR: No valid nozzle diameter set.");
1112 ASSERT(breakupLength > -MFloatEps, "ERROR: Invalid breakup length.");
1113
1114 // 1) Loop over all particles and check for break-up
1115 for(MInt i = 0; i < noPart; i++) {
1116 auto& droplet = s_backPtr->m_partList.at(i);
1117
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;
1124
1125 ASSERT(droplet.m_cellId > 0, "ERROR: Invalid cellId");
1126
1127 // increase breakup time
1128 droplet.m_breakUpTime += dt;
1129
1130 ASSERT(!isnan(droplet.m_shedDiam) && droplet.m_shedDiam > -MFloatEps, "");
1131
1132 // 2) Computations for Kelvin-Helmholtz Break up Model mainly based on:
1133 // "Modeling Spray Atomization with the Kelvin-Helmholtz/Rayleigh-Taylor Hibrid Model"
1134 // by Beale and Reitz, Atomization and Sprays (1999)
1135 // NOTE: in this case the radius and not the diameter is used as reference length!
1136
1137 MFloat liquidLength = MFloatMax;
1138 if(m_RTsecBreakUpLength) {
1139 liquidLength = maia::math::distance(s_backPtr->m_spawnCoord, droplet.m_position);
1140 /*
1141 if(m_injectorType == MULTIHOLE) {
1142 const MFloat dist1 = maia::math::distance(s_backPtr->m_spawnCoord, droplet.m_position);
1143 const MFloat diff = POW2(dist1) - POW2(m_injectorDiameter);
1144 if(diff > 0) {
1145 liquidLength = sqrt(diff);
1146 } else {
1147 liquidLength = 0;
1148 }
1149 }
1150 */
1151 }
1152
1153 const MFloat dropRadius = 0.5 * droplet.m_diameter;
1154 const MFloat relV = droplet.magRelVel(&droplet.m_fluidVel[0], &droplet.m_velocity[0]);
1155
1156 //#ifdef LPT_DEBUG
1157 if(isnan(relV)) {
1158 cerr << droplet.m_fluidVel[0] << " " << droplet.m_fluidVel[1] << " " << droplet.m_fluidVel[nDim - 1] << endl;
1159 mTerm(1, AT_, "");
1160 }
1161 //#endif
1162
1163 // radius based particle Weber-number
1164 const MFloat We_l =
1165 droplet.WeberNumber(s_backPtr->m_material->density(droplet.m_temperature), pow(relV, 2), droplet.m_temperature)
1166 * droplet.s_We;
1167
1168 // if the weber number is to small the further calculated values are nan or inf
1169 if(We_l < m_weberLimitKH) continue;
1170
1171 // particle Reynolds number
1172 // factor 0.5 to convert from diameter to radius based
1173 const MFloat Re_l = 0.5
1174 * droplet.particleRe(material().density(droplet.m_temperature), relV,
1175 material().dynamicViscosity(droplet.m_temperature))
1176 * droplet.s_Re;
1177
1178 // gas Weber number
1179 const MFloat We_g = droplet.WeberNumber(droplet.m_fluidDensity, pow(relV, 2), droplet.m_temperature) * droplet.s_We;
1180 // Ohnesorge number
1181 const MFloat Z = sqrt(We_l) / Re_l;
1182 // Taylor number
1183 const MFloat T = Z * sqrt(We_g);
1184
1185 // Lamda_KH
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);
1188
1189
1190 // Kelvin-Helmholtz reference diameter (thus factor 2.0)
1191 const MFloat KH_diameter = 2.0 * m_B0 * KH_waveL;
1192
1193 // ASSERT(KH_diameter > MFloatEps,
1194 // to_string(KH_diameter) + " surfWL " + to_string(KH_waveL));
1195 if(KH_diameter < 2 * MFloatEps) {
1196 cerr << KH_diameter << " d_KH " << dropRadius << " " << Z << " " << T << " " << We_g << " " << We_l << endl;
1197 }
1198
1199 // Omega_KH
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);
1204
1205
1206 // tau_kH
1207 const MFloat KH_time = 3.726 * m_B1 * dropRadius / (KH_waveL * KH_growthR);
1208
1209 // const MFloat KH_speed = (droplet.m_diameter - KH_diameter) / KH_time;
1210
1211 // 3) prepare initialisation of new droplets
1212 // parent velocity magnitude
1213 const MFloat magVel = droplet.magVel();
1214 array<MFloat, nDim> parentTrajectory{};
1215 // get parent droplet direction
1216 for(MInt v = 0; v < nDim; v++) {
1217 parentTrajectory.at(v) = droplet.m_velocity.at(v) / magVel;
1218 }
1219
1220 // in case parent parcel is not moving
1221 if(magVel < MFloatEps) {
1222 // TODO labels:LPT,totest check if this exception is necessary!
1223 continue;
1224 /*
1225 // get direction from flow direction
1226 MFloat fluidVel = 0;
1227 for(MInt j = 0; j < nDim; j++) {
1228 fluidVel += POW2(s_backPtr->a_fluidVelocity(droplet.m_cellId, j));
1229 }
1230 fluidVel = sqrt(fluidVel);
1231
1232 for(MInt v = 0; v < nDim; v++) {
1233 parentTrajectory[v] = s_backPtr->a_fluidVelocity(droplet.m_cellId, v) / fluidVel;
1234 }
1235 */
1236 }
1237
1238 // magnitude of droplet acceleration
1239 MFloat drop_accel = 0.0;
1240 for(MInt dir = 0; dir < nDim; dir++) {
1241 drop_accel += POW2(droplet.m_accel.at(dir));
1242 }
1243 drop_accel = sqrt(drop_accel);
1244
1245 // storing some old droplet properties
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];
1250 }
1251 const MFloat oldDiameter = droplet.m_diameter;
1252
1253 // 4) Computations for Rayleigh-Taylor Break up Model mainly based on:
1254 // "Modeling Spray Atomization with the Kelvin-Helmholtz/Rayleigh-Taylor Hibrid Model"
1255 // by Beale and Reitz, Atomization and Sprays (1999)
1256 // or Baumgarten
1257
1258 // NOTE wavelength = 2 * PI * 1/ waveNumber(=K_RT)
1259 const MFloat RT_waveL =
1260 2.0 * PI
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);
1264
1265 // Omega_RT
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);
1271
1272 const MFloat RT_time = m_CT / RT_frequency;
1273
1274 // NOTE: in this case its the diameter!
1275 const MFloat RT_diameter = m_CRT * RT_waveL;
1276
1277 ASSERT(RT_diameter > 0, "");
1278
1279 // const MFloat RT_speed = (droplet.m_diameter - RT_diameter) / RT_time;
1280
1281
1282 // 6) Rayleigh-Taylor break-up
1283 if(m_RTsecBreakUp && liquidLength >= breakupLength && droplet.m_breakUpTime >= RT_time
1284 && droplet.m_diameter > RT_diameter && We_l < m_weberLimitRT) {
1285 // choose RT child diameter version
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) {
1294 // no breakup
1295 continue;
1296 }
1297 } else if(m_RTDiameterMode == 1) {
1298 childDiameter = pow(oldDiameter, 2 / 3) * pow(RT_diameter, 1 / 3);
1299 }
1300
1301 if(m_maxRTChildParcels == 0) {
1302 // no child parcles are added, only the diameter and no-particles is changed
1303 auto noDroplets = static_cast<MInt>(droplet.m_noParticles * POW3(oldDiameter / childDiameter));
1304 // NOTE: casting obove means a runding down for all positive values
1305 // version below applies additiona round up for values above 0.75
1306
1307 // NOTE: applying rounding to avoid mass-losses!
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;
1315 m_noRTsecBreakUp++;
1316
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;
1320 }
1321 }
1322
1323 } else {
1324 if(m_RTDiameterMode != 2) {
1325 initPRNG(globalTimeStep, s_backPtr->c_globalId(droplet.m_cellId));
1326 }
1327
1328 // calculate new number of droplets
1329 const MFloat dropletVolume = pow(droplet.m_diameter, 3.0);
1330 const MFloat dropletRTVolume = pow(childDiameter, 3.0);
1331
1332 MInt noBrokenUpDroplets = dropletVolume / dropletRTVolume;
1333
1334 // don't allow inconsequential breakup
1335 if(noBrokenUpDroplets == 1) {
1336 continue;
1337 }
1338 noBrokenUpDroplets *= droplet.m_noParticles;
1339
1340 // number of drops per new parcel
1341 auto dropsPerParcel = static_cast<MInt>(
1342 noBrokenUpDroplets > m_maxRTChildParcels ? floor(noBrokenUpDroplets / m_maxRTChildParcels) : 1);
1343
1344 // limit the number of new parcels
1345 MInt noNewParcels = noBrokenUpDroplets;
1346 if(dropsPerParcel > 1) {
1347 noNewParcels = m_maxRTChildParcels - 1;
1348 } else {
1349 --noNewParcels;
1350 }
1351
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);
1357
1358 ASSERT(parentMass > 0, "Invalid mass!");
1359 ASSERT(parentDiam > 0 && !isnan(parentDiam), "");
1360 ASSERT(dropsPerParcel > 0, "");
1361
1362 // Reset values
1363 droplet.m_diameter = parentDiam;
1364 droplet.m_shedDiam = parentDiam;
1365 droplet.m_noParticles = dropsPerParcel;
1366 droplet.m_breakUpTime = 0.0;
1367 m_noRTsecBreakUp++;
1368
1369#ifdef _OPENMP
1370#pragma omp critical
1371#endif
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;
1382 }
1383
1384 MFloat sprayAngle = m_sprayAngle[0];
1385 if(m_hollowConeInjector) {
1386 // hollow cone injectors have a large spray angle because of the hollow cone
1387 // geometry the actual relevant spray angle is much smaller
1388 static constexpr MFloat angularGapAngle = 45.0;
1389 sprayAngle = m_sprayAngle[0] - 2 * angularGapAngle;
1390 }
1391
1392#ifdef LPT_DEBUG
1393 vector<MFloat> mom(nDim);
1394#endif
1395
1396 // create new droplets
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);
1401
1402 const MInt childId = s_backPtr->addParticle(droplet.m_cellId, childDiameter, droplet.m_densityRatio, 0, 3,
1403 &childVelocity[0], &droplet.m_position[0], dropsPerParcel);
1404
1405 // set child temperature based on parent value
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;
1409
1410
1411#ifdef LPT_DEBUG
1412 // check kin. energy conservation
1413 MFloat magVelChild = 0;
1414 for(MInt n = 0; n < nDim; n++) {
1415 magVelChild += POW2(childVelocity[n]);
1416 }
1417 magVelChild = sqrt(magVelChild);
1418 if(fabs(magVel - magVelChild) > MFloatEps) {
1419 cerr << "RT kin. energy loss: " << magVel << " " << magVelChild << endl;
1420 }
1421 for(MInt n = 0; n < nDim; n++) {
1422 mom[n] += sphericalMass(childDiameter, droplet.m_temperature) * dropsPerParcel * childVelocity[n];
1423 }
1424#endif
1425 }
1426
1427#ifdef LPT_DEBUG
1428 // check momentum conservation:
1429 for(MInt n = 0; n < nDim; n++) {
1430 mom[n] +=
1431 droplet.m_velocity[n] * sphericalMass(droplet.m_diameter, droplet.m_temperature) * droplet.m_noParticles;
1432
1433 if(fabs(oldMom[n] - mom[n]) > MFloatEps) {
1434 cerr << "RT momentum change: " << oldMom[n] << " " << mom[n] << endl;
1435 }
1436 }
1437
1438 // check conservation:
1439 const MFloat newMass = sphericalMass(droplet.m_diameter, droplet.m_temperature) * droplet.m_noParticles;
1440 MFloat childMass = 0;
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);
1445 }
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;
1450 }
1451#endif
1452 }
1453 continue;
1454 }
1455
1456 // Kelvin-Helmholtz induced break-up
1457 // if((KH_speed >= RT_speed || liquidLength < breakupLength) && droplet.m_diameter >= KH_diameter && m_KHsecBreakUp)
1458 // {
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;
1462 // ASSERT(KH_mass > MFloatEps, to_string(KH_mass) + " KH_diameter:" + to_string(KH_diameter));
1463
1464 // Kelvin-Helmholtz induced mass stripping
1465 // NOTE: shedDiam is the diameter of the parent particle after undergoing KH-breakup!
1466 // i.e. the remaining mass in the parent!
1467 droplet.m_shedDiam -= (droplet.m_shedDiam - KH_diameter) / KH_time * dt;
1468
1469 // if the shed diameter is smaller than the KH diameter
1470 // no child parcel is added, the parent diameter is set to the KH-diameter
1471 // and the noParticles is set by mass conservation
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);
1477 // fully transversion to KH-diameter
1478 // after applying rounding to avoid mass-losses!
1479 // droplet.m_breakUpTime = -dt;
1480 droplet.m_diameter = diameter;
1481 droplet.m_shedDiam = diameter;
1482 droplet.m_noParticles = noDroplets;
1483 m_noKHsecBreakUp++;
1484
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;
1488 }
1489 continue;
1490 }
1491 }
1492
1493 // mass going into the particle = oldMass - mass remaining in the parant parcel
1494 const MFloat sheddedMass =
1495 oldMass - sphericalMass(droplet.m_shedDiam, droplet.m_temperature) * droplet.m_noParticles;
1496
1497 if(m_massLimit > sheddedMass / oldMass) {
1498 continue;
1499 }
1500
1501 MInt noSheddedOffDrops = floor(sheddedMass / KH_mass);
1502 if(noSheddedOffDrops <= 0) continue;
1503 const MFloat parentMass = oldMass - KH_mass * noSheddedOffDrops;
1504 ASSERT(parentMass > 0, "");
1505 // when reformulating the equations above, parentDiam is equal to the shedDiam!
1506 const MFloat parentDiam =
1507 pow(6.0 / PI * parentMass / (droplet.m_noParticles * material().density(droplet.m_temperature)), 1.0 / 3.0);
1508
1509 // adding additional child parcel
1510 // child parcel properties set as prescribed in
1511 //"Modeling Atomization Processes in High-Pressure Vaporizing Sprays"
1512 // R. Reitz , Atom.& Sprayz 3 (1987) 309-337
1513 initPRNG(globalTimeStep, s_backPtr->c_globalId(droplet.m_cellId));
1514
1515 // droplet.m_breakUpTime = -dt;
1516 droplet.m_diameter = parentDiam;
1517 droplet.m_shedDiam = parentDiam;
1518
1519 // different spray-angle versions
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;
1525 }
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];
1532 } else {
1533 angle = m_sprayAngle[1];
1534 }
1535 }
1536
1537 if(angle < 0 || angle > 90) {
1538 cerr << "Large KH-spray angle : " << angle << endl;
1539 }
1540
1541 array<MFloat, nDim> childVelocity{};
1542 randomVectorInCone(&childVelocity[0], &parentTrajectory[0], magVel, angle, string2enum("PART_EMITT_DIST_NONE"),
1543 randomSecBreakUp(), m_injectionDistSigmaCoeff);
1544
1545 // adding a single additional parcel with KH_diameter
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);
1548
1549 // set child temperature based on parent value
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;
1553
1554
1555 m_noKHsecBreakUp++;
1556
1557#ifdef _OPENMP
1558#pragma omp critical
1559#endif
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;
1571 }
1572
1573 // apply velocity change to parent particle
1574 // under momentum and kin. energy conservation!
1575 /*
1576 vector<MFloat> childMom(nDim);
1577 for(MInt n = 0; n < nDim; n++){
1578 childMom[n] = sphericalMass(KH_diameter) * noSheddedOffDrops * childVelocity[n];
1579 droplet.m_velocity[n] = (oldMom[n] - childMom[n])
1580 / (sphericalMass(droplet.m_diameter) * droplet.m_noParticles);
1581 }
1582
1583 maia::math::normalize(&droplet.m_velocity[0], 3);
1584 for(MInt n = 0; i < nDim; n++){
1585 droplet.m_velocity[i] *= magVel;
1586 }
1587 */
1588
1589#ifdef LPT_DEBUG
1590 // check momentum conservation
1591 vector<MFloat> mom(nDim);
1592 for(MInt n = 0; n < nDim; n++) {
1593 mom[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
1600 << endl;
1601 }
1602 }
1603
1604 // check kin. energy conservation
1605 MFloat magVelChild = 0;
1606 for(MInt n = 0; n < nDim; n++) {
1607 magVelChild += POW2(childVelocity[n]);
1608 }
1609 magVelChild = sqrt(magVelChild);
1610 if(fabs(magVel - magVelChild) > MFloatEps) {
1611 cerr << "KH kin. energy loss: " << magVel << " " << magVelChild << endl;
1612 }
1613 MFloat magVelNew = droplet.magVel();
1614 if(fabs(magVel - magVelNew) > MFloatEps) {
1615 cerr << "KH kin. energy loss: " << magVel << " " << magVelNew << endl;
1616 }
1617 // check mass conservation
1618 const MFloat newMass = sphericalMass(droplet.m_diameter, droplet.m_temperature) * droplet.m_noParticles;
1619 const MInt childId = s_backPtr->a_noParticles() - 1;
1620 const MFloat childMass =
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;
1627 }
1628#endif
1629 } /*else if( m_KHsecBreakUp && KH_diameter > droplet.m_diameter &&
1630 droplet.m_breakUpTime > KH_time ) {
1631 //version by Reitz (1987)
1632 MFloat d_KH = 2 * mMin(
1633 pow(3 * PI * POW2(droplet.m_diameter) * relV / ( 2 * 4 * KH_growthR), 0.33),
1634 pow(3 * POW2(droplet.m_diameter) * KH_waveL / 16, 0.33));
1635
1636 auto noDroplets = static_cast<MInt>(droplet.m_noParticles *
1637 POW3(droplet.m_diameter) / POW3(d_KH));
1638 if(noDroplets > 0 ) {
1639 const MFloat diameter = pow(6.0 / PI, 1.0 / 3.0) *
1640 pow(oldMass / (noDroplets * material().density(droplet.m_temperature)), 1.0/ 3.0);
1641 //change droplet diameter
1642 droplet.m_breakUpTime = 0.0;
1643 droplet.m_diameter = diameter;
1644 droplet.m_shedDiam = diameter;
1645 droplet.m_noParticles = noDroplets;
1646 m_noKHsecBreakUp++;
1647 }
1648 }
1649 */
1650 }
1651}
1652
1657template <MInt nDim>
1659 const MInt numInjections = m_injectionRateList.dim0();
1660 if(numInjections == 0) return;
1661
1662 if(m_injectionTimings[0] > m_timeSinceSOI) {
1663 m_currentInjectionRate = (m_timeSinceSOI + 1) / (m_injectionTimings[0] + 1) * m_injectionRateList[0];
1664 } else {
1665 for(MInt i = 0; i <= numInjections; i++) {
1666 if(m_injectionTimings[i] > m_timeSinceSOI) {
1667 const MFloat t =
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]);
1670 break;
1671 }
1672 }
1673 }
1674}
1675
1676
1695template <MInt nDim>
1696void SprayModel<nDim>::injectParticles(const MInt spawnParticlesCount, const MFloat* spawnCoord,
1697 const MFloat particleDiameter, const MFloat particleDensityRatio,
1698 const MFloat particleVelocity, const MFloat sprayConeAngle,
1699 const MFloat coneDiameter, std::function<MFloat(MInt)> holePositionAngle,
1700 std::function<MFloat(MInt)> injectorDiameter,
1701 std::function<MFloat(MInt)> holeNozzleAngle, const MBool hull,
1702 const MInt parcelSize, const MInt noInjectionHoles) {
1703 std::array<MFloat, nDim> spawnParticlesInitVelo;
1704 array<MFloat, 3> injectionLocation{};
1705 array<MFloat, 3> holeLocation{};
1706
1707 if(spawnParticlesCount > 0) {
1708 m_injStep++;
1709 }
1710
1711 if(particleDiameter < s_backPtr->m_sizeLimit) {
1712 cerr << "Inj. small particle " << particleDiameter << " " << s_backPtr->m_sizeLimit << endl;
1713 }
1714
1715 for(MInt i = 0; i < spawnParticlesCount; i++) {
1716 MInt randomShift = 1;
1717 if(m_partEmittDist != PART_EMITT_DIST_NONE) {
1718 randomShift = 0;
1719 m_PRNGPrimBreakUpCount +=
1720 randomVectorInCone(&spawnParticlesInitVelo[0], &m_injectorDir[0], particleVelocity, sprayConeAngle,
1721 m_partEmittDist, randomPrimBreakUp(0), m_injectionDistSigmaCoeff, 0);
1722
1723 if(coneDiameter > numeric_limits<MFloat>::epsilon() && !hull) {
1724 // determine a random point within the injector orfice
1725 randomPointInCircle(&holeLocation[0], &m_injectorDir[0], coneDiameter, randomPrimBreakUp(2));
1726
1727 if(noInjectionHoles > 1) {
1728 // determine position of the hole
1729 pointOnCircle(&injectionLocation[0], &m_injectorDir[0], injectorDiameter(i), holePositionAngle(i));
1730 }
1731 } else if(hull) {
1732 m_PRNGPrimBreakUpCount +=
1733 randomVectorInCone(&spawnParticlesInitVelo[0], &m_injectorDir[0], particleVelocity, sprayConeAngle,
1734 m_partEmittDist, randomPrimBreakUp(0), m_injectionDistSigmaCoeff, holeNozzleAngle(i));
1735
1736 // determine a point just on the hull of the cone (like for a hollow-cone injector)
1737 // TODO labels:LPT make settable
1738 static constexpr MBool randomDist = false;
1739 if(randomDist) {
1740 randomPointOnCircle(&injectionLocation[0], &m_injectorDir[0], coneDiameter, randomPrimBreakUp(1),
1741 spawnParticlesCount, i);
1742 } else {
1743 // generate one particle per 1 deg
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);
1747 }
1748 }
1749 }
1750
1751 // rotate injection by the angle of the nozzle
1752 if(fabs(holeNozzleAngle(i)) > numeric_limits<MFloat>::epsilon()) {
1753 const MFloat holeAngleRad = -holeNozzleAngle(i) / 180 * M_PI;
1754 const MFloat absDist = maia::math::norm(injectionLocation);
1755
1756 array<MFloat, nDim> crossP{};
1757 array<MFloat, nDim> revInjLoc{};
1758 for(MInt j = 0; j < nDim; j++) {
1759 revInjLoc[j] = -injectionLocation[j] / absDist;
1760 }
1761
1762 maia::math::cross(&m_injectorDir[0], &revInjLoc[0], &crossP[0]);
1763
1764 // boundary basis at the injection location
1765 MFloatScratchSpace injLocBasis(nDim, nDim, FUN_, "R");
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];
1770 }
1771
1772 // find inverse for reverse-tranformation
1773 MFloatScratchSpace inverse(nDim, nDim, AT_, "inverse");
1774 for(MInt j = 0; j < nDim; j++) {
1775 for(MInt k = 0; k < nDim; k++) {
1776 inverse(j, k) = injLocBasis(j, k);
1777 }
1778 }
1779
1780 maia::math::invert(&inverse(0, 0), 3, 3);
1781 // NOTE: inverse might not be normalized correctly, meaning that the magnitude is not conserved!
1782
1783 // rotation matrix around the z-axis in the injection location basis
1784 MFloatScratchSpace R(nDim, nDim, FUN_, "R");
1785 R.fill(0.0);
1786 R(2, 2) = 1.0;
1787 R(0, 0) = cos(holeAngleRad);
1788 R(1, 0) = sin(holeAngleRad);
1789 R(0, 1) = -sin(holeAngleRad);
1790 R(1, 1) = cos(holeAngleRad);
1791
1792 // global coordinate transformation is then found be B*R*B^-1
1793 MFloatScratchSpace result1(nDim, nDim, FUN_, "result1");
1794 MFloatScratchSpace result(nDim, nDim, FUN_, "result");
1795 MFloatScratchSpace tempV(nDim, FUN_, "tempV");
1796 for(MInt j = 0; j < nDim; j++) {
1797 tempV[j] = spawnParticlesInitVelo[j];
1798 }
1799
1800 maia::math::multiplyMatricesSq(injLocBasis, R, result1, nDim);
1801 maia::math::multiplyMatricesSq(result1, inverse, result, nDim);
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];
1806 }
1807 }
1808
1809 // reset vector length and direction
1810 maia::math::normalize(&spawnParticlesInitVelo[0], nDim);
1811 // reverse sign if velocity has opposite direction to the injector direction
1812 if(scalarProduct(&spawnParticlesInitVelo[0], &m_injectorDir[0], nDim) < 0.0) {
1813 for(MInt j = 0; j < nDim; j++) {
1814 spawnParticlesInitVelo[j] *= -1;
1815 }
1816 }
1817 for(MInt j = 0; j < nDim; j++) {
1818 spawnParticlesInitVelo[j] *= particleVelocity;
1819 }
1820 }
1821 // final injection position = hole position + orifice position + injector center
1822 for(MInt j = 0; j < nDim; j++) {
1823 injectionLocation[j] += spawnCoord[j] + holeLocation[j];
1824 }
1825
1826 s_backPtr->addParticle(s_backPtr->m_spawnCellId, particleDiameter, particleDensityRatio, randomShift, 2,
1827 &spawnParticlesInitVelo[0], &injectionLocation[0], parcelSize);
1828 }
1829}
1830
1836template <MInt nDim>
1838 if(domainId() != 0) return;
1839
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;
1844
1845 std::mt19937_64 randomNumberGenerator;
1846 randomNumberGenerator.seed(m_spraySeed);
1847
1848 const MInt noBinns = 1000;
1849 MFloat binWidth = (maxValue - minValue) / noBinns;
1850 MIntScratchSpace binValue(noBinns, AT_, "binValue");
1851 binValue.fill(0);
1852
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);
1856 binValue[binId]++;
1857 }
1858
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;
1866
1867 for(MInt i = 0; i < noBinns; i++) {
1868 m_log << minValue + i * binWidth << " " << binValue[i] << endl;
1869 }
1870
1871 minValue = 0.0;
1872 maxValue = 5 * meanValue;
1873 binWidth = (maxValue - minValue) / noBinns;
1874
1875 binValue.fill(0);
1876
1877 for(MInt i = 0; i < noDraws; i++) {
1878 const MFloat drawValue = NTDistribution(meanValue, randomNumberGenerator);
1879 const MInt binId = std::floor((drawValue - minValue) / binWidth);
1880 binValue[binId]++;
1881 }
1882
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;
1887 }
1888}
1889// Explicit instantiations for 2D and 3D
1890template class SprayModel<3>;
void mAlloc(T *&a, const MLong N, const MString &objectName, MString function)
allocates memory for one-dimensional array 'a' of size N
Definition: alloc.h:173
static MInt propertyLength(const MString &name, MInt solverId=m_noSolvers)
Returns the number of elements of a property.
Definition: context.cpp:538
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
Definition: context.cpp:494
Definition: lpt.h:82
std::vector< LPTSpherical< nDim > > m_partList
Definition: lpt.h:527
This class is a ScratchSpace.
Definition: scratch.h:758
void fill(T val)
fill the scratch with a given value
Definition: scratch.h:311
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.
Definition: lptspray.cpp:1696
void logDistributionFunction()
create log of distribution functions for plots and debugging
Definition: lptspray.cpp:1837
void injection(const MFloat dt)
Perform a spray injection.
Definition: lptspray.cpp:672
void primaryBreakUp(const MFloat dt)
Inject new particles from the injector exit.
Definition: lptspray.cpp:744
void updateInjectionRate()
Update injection rate.
Definition: lptspray.cpp:1658
void init(LPT< nDim > *_particle)
Spray model initialisation.
Definition: lptspray.cpp:29
void readSprayProperties()
void secondaryBreakUp(const MFloat dt)
Atomization of particles.
Definition: lptspray.cpp:1081
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
Definition: enums.cpp:20
@ PART_EMITT_DIST_NONE
Definition: enums.h:348
@ HOLLOWCONE
Definition: enums.h:361
@ FULLCONE
Definition: enums.h:361
@ MULTIHOLE_OPT
Definition: enums.h:361
@ MULTIHOLE
Definition: enums.h:361
@ MULTIHOLE_TME
Definition: enums.h:361
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
constexpr Real POW3(const Real x)
Definition: functions.h:123
constexpr Real POW2(const Real x)
Definition: functions.h:119
constexpr T mMax(const T &x, const T &y)
Definition: functions.h:94
MInt globalTimeStep
InfoOutFile m_log
std::ostream cerr0
int32_t MInt
Definition: maiatypes.h:62
std::basic_string< char > MString
Definition: maiatypes.h:55
double MFloat
Definition: maiatypes.h:52
int64_t MLong
Definition: maiatypes.h:64
bool MBool
Definition: maiatypes.h:58
MInt id
Definition: maiatypes.h:71
Definition: lpt.h:47
MFloat NTDistribution(const MFloat x_mean, std::mt19937_64 &PRNG)
Definition: lptlib.h:601
MFloat scalarProduct(const MFloat *a, const MFloat *b, const MInt length)
Definition: lptlib.h:149
void randomPointInCircle(MFloat *vec, const MFloat *normalDirection, const MFloat diameter, std::mt19937_64 &PRNG)
Definition: lptlib.h:383
void randomPointOnCircle(MFloat *vec, const MFloat *normalDirection, const MFloat diameter, std::mt19937_64 &PRNG, const MInt circleSplit=1, const MInt splitNo=0)
Definition: lptlib.h:459
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.
Definition: lptlib.h:252
MFloat rosinRammler(const MFloat min, const MFloat mean, const MFloat max, const MFloat spread, std::mt19937_64 &PRNG)
Definition: lptlib.h:584
void pointOnCircle(MFloat *vec, const MFloat *normalDirection, const MFloat diameter, MFloat phi)
Definition: lptlib.h:528
void cross(const T *const u, const T *const v, T *const c)
Definition: maiamath.h:101
MFloat distance(const MFloat *a, const MFloat *b)
Definition: maiamath.h:249
void invert(MFloat *A, const MInt m, const MInt n)
Definition: maiamath.cpp:171
void normalize(std::array< T, N > &u)
Definition: maiamath.h:191
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,...
Definition: maiamath.h:506
MFloat norm(const std::array< T, N > &u)
Definition: maiamath.h:148
void multiplyMatricesSq(MFloatScratchSpace &m1, MFloatScratchSpace &m2, MFloatScratchSpace &result, MInt dim)
Definition: maiamath.h:279
Namespace for auxiliary functions/classes.