BUG: UListIO: byteSize overflowing on really big faceLists
[OpenFOAM-2.0.x.git] / src / lagrangian / dieselSpray / injector / swirlInjector / swirlInjector.C
blob7609535dc2778a1b89c23ca0d42bec03ca0d4b9b
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
11     OpenFOAM is free software: you can redistribute it and/or modify it
12     under the terms of the GNU General Public License as published by
13     the Free Software Foundation, either version 3 of the License, or
14     (at your option) any later version.
16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
19     for more details.
21     You should have received a copy of the GNU General Public License
22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
24 \*---------------------------------------------------------------------------*/
26 #include "swirlInjector.H"
27 #include "addToRunTimeSelectionTable.H"
28 #include "mathematicalConstants.H"
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 namespace Foam
34     defineTypeNameAndDebug(swirlInjector, 0);
36     addToRunTimeSelectionTable
37     (
38         injectorType,
39         swirlInjector,
40         dictionary
41     );
45 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
47 Foam::swirlInjector::swirlInjector
49     const Foam::Time& t,
50     const Foam::dictionary& dict
53     injectorType(t, dict),
54     propsDict_(dict.subDict(typeName + "Props")),
55     position_(propsDict_.lookup("position")),
56     direction_(propsDict_.lookup("direction")),
57     d_(readScalar(propsDict_.lookup("diameter"))),
58     mass_(readScalar(propsDict_.lookup("mass"))),
59     injectionPressure_(readScalar(propsDict_.lookup("injectionPressure"))),
60     T_(readScalar(propsDict_.lookup("temperature"))),
61     nParcels_(readLabel(propsDict_.lookup("nParcels"))),
62     X_(propsDict_.lookup("X")),
63     massFlowRateProfile_(propsDict_.lookup("massFlowRateProfile")),
64     injectionPressureProfile_(propsDict_.lookup("injectionPressureProfile")),
65     velocityProfile_(massFlowRateProfile_),
66     CdProfile_(massFlowRateProfile_),
67     TProfile_(massFlowRateProfile_),
68     averageParcelMass_(mass_/nParcels_),
69     pressureIndependentVelocity_(false)
71     // convert CA to real time
72     forAll(massFlowRateProfile_, i)
73     {
74         massFlowRateProfile_[i][0] =
75             t.userTimeToTime(massFlowRateProfile_[i][0]);
76     }
77     forAll(injectionPressureProfile_, i)
78     {
79         injectionPressureProfile_[i][0] =
80             t.userTimeToTime(injectionPressureProfile_[i][0]);
81     }
83     // check if time entries match
84     if
85     (
86         mag(massFlowRateProfile_[0][0] - injectionPressureProfile_[0][0])
87       > SMALL
88     )
89     {
90         FatalErrorIn
91         (
92             "swirlInjector::swirlInjector"
93             "(const time& t, const dictionary dict)"
94         )   << "Start-times do not match for "
95             << "injectionPressureProfile and massFlowRateProfile."
96             << abort(FatalError);
97     }
99     // check if time entries match
100     if
101     (
102         mag
103         (
104             massFlowRateProfile_.last()[0]
105           - injectionPressureProfile_.last()[0]
106         ) > SMALL
107     )
108     {
109         FatalErrorIn
110         (
111             "swirlInjector::swirlInjector"
112             "(const time& t, const dictionary dict)"
113         )   << "End-times do not match for "
114             << "injectionPressureProfile and massFlowRateProfile."
115             << abort(FatalError);
116     }
118     scalar integratedMFR = integrateTable(massFlowRateProfile_);
119     scalar integratedPressure =
120         integrateTable(injectionPressureProfile_)/(teoi()-tsoi());
122     forAll(massFlowRateProfile_, i)
123     {
124         // correct the massFlowRateProfile to match the injected mass
125         massFlowRateProfile_[i][1] *= mass_/integratedMFR;
127         velocityProfile_[i][0] = massFlowRateProfile_[i][0];
129         TProfile_[i][0] = massFlowRateProfile_[i][0];
130         TProfile_[i][1] = T_;
132         CdProfile_[i][0] = massFlowRateProfile_[i][0];
134     }
136     forAll(injectionPressureProfile_, i)
137     {
138         // correct the pressureProfile to match the injection pressure
139         injectionPressureProfile_[i][1] *=
140             injectionPressure_/integratedPressure;
141     }
143     // Normalize the direction vector
144     direction_ /= mag(direction_);
146     setTangentialVectors();
148     // check molar fractions
149     scalar Xsum = 0.0;
150     forAll(X_, i)
151     {
152         Xsum += X_[i];
153     }
155     if (mag(Xsum - 1.0) > SMALL)
156     {
157         WarningIn
158         (
159             "swirlInjector::swirlInjector"
160             "(const time& t, const dictionary dict)"
161         )   << "X does not add up to 1.0, correcting molar fractions." << endl;
163         forAll(X_, i)
164         {
165             X_[i] /= Xsum;
166         }
167     }
171 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
173 Foam::swirlInjector::~swirlInjector()
177 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
179 void Foam::swirlInjector::setTangentialVectors()
181     cachedRandom rndGen(label(0), -1);
182     scalar magV = 0.0;
183     vector tangent;
185     while (magV < SMALL)
186     {
187         vector testThis = rndGen.sample01<vector>();
189         tangent = testThis - (testThis & direction_)*direction_;
190         magV = mag(tangent);
191     }
193     tangentialInjectionVector1_ = tangent/magV;
194     tangentialInjectionVector2_ = direction_^tangentialInjectionVector1_;
198 Foam::label Foam::swirlInjector::nParcelsToInject
200     const scalar time0,
201     const scalar time1
202 ) const
205     scalar mInj =
206         mass_*(fractionOfInjection(time1) - fractionOfInjection(time0));
207     label nParcels = label(mInj/averageParcelMass_ + 0.49);
209     return nParcels;
213 const Foam::vector Foam::swirlInjector::position(const label n) const
215     return position_;
219 Foam::vector Foam::swirlInjector::position
221     const label n,
222     const scalar time,
223     const bool twoD,
224     const scalar angleOfWedge,
225     const vector& axisOfSymmetry,
226     const vector& axisOfWedge,
227     const vector& axisOfWedgeNormal,
228     cachedRandom& rndGen
229 ) const
231     if (twoD)
232     {
233         scalar is = position_ & axisOfSymmetry;
234         scalar magInj = mag(position_ - is*axisOfSymmetry);
236         vector halfWedge =
237             axisOfWedge*cos(0.5*angleOfWedge)
238           + axisOfWedgeNormal*sin(0.5*angleOfWedge);
239         halfWedge /= mag(halfWedge);
241         return (is*axisOfSymmetry + magInj*halfWedge);
242     }
243     else
244     {
245         // otherwise, disc injection
246         scalar iRadius = d_*rndGen.sample01<scalar>();
247         scalar iAngle = constant::mathematical::twoPi*rndGen.sample01<scalar>();
249         return
250         (
251             position_
252           + iRadius
253           * (
254               tangentialInjectionVector1_*cos(iAngle)
255             + tangentialInjectionVector2_*sin(iAngle)
256           )
257         );
259     }
261     return position_;
265 Foam::label Foam::swirlInjector::nHoles() const
267     return 1;
271 Foam::scalar Foam::swirlInjector::d() const
273     return d_;
277 const Foam::vector& Foam::swirlInjector::direction
279     const label i,
280     const scalar time
281 ) const
283     return direction_;
287 Foam::scalar Foam::swirlInjector::mass
289     const scalar time0,
290     const scalar time1,
291     const bool twoD,
292     const scalar angleOfWedge
293 ) const
295     scalar mInj =
296         mass_*(fractionOfInjection(time1) - fractionOfInjection(time0));
298     // correct mass if calculation is 2D
299     if (twoD)
300     {
301         mInj *= 0.5*angleOfWedge/constant::mathematical::pi;
302     }
304     return mInj;
308 Foam::scalar Foam::swirlInjector::mass() const
310     return mass_;
314 Foam::List<Foam::swirlInjector::pair>
315 Foam::swirlInjector::massFlowRateProfile() const
317     return massFlowRateProfile_;
321 Foam::scalar Foam::swirlInjector::massFlowRate(const scalar time) const
323     return getTableValue(massFlowRateProfile_, time);
327 Foam::List<Foam::swirlInjector::pair>
328 Foam::swirlInjector::injectionPressureProfile() const
330     return injectionPressureProfile_;
334 Foam::scalar Foam::swirlInjector::injectionPressure(const scalar time) const
336     return getTableValue(injectionPressureProfile_, time);
340 Foam::List<Foam::swirlInjector::pair>
341 Foam::swirlInjector::velocityProfile() const
343     return velocityProfile_;
347 Foam::scalar Foam::swirlInjector::velocity(const scalar time) const
349     return getTableValue(velocityProfile_, time);
353 Foam::List<Foam::swirlInjector::pair> Foam::swirlInjector::CdProfile() const
355     return CdProfile_;
359 Foam::scalar Foam::swirlInjector::Cd(const scalar time) const
361     return getTableValue(CdProfile_, time);
365 const Foam::scalarField& Foam::swirlInjector::X() const
367     return X_;
371 Foam::List<Foam::swirlInjector::pair> Foam::swirlInjector::T() const
373     return TProfile_;
377 Foam::scalar Foam::swirlInjector::T(const scalar time) const
379     return T_;
383 Foam::scalar Foam::swirlInjector::tsoi() const
385     return massFlowRateProfile_.first()[0];
389 Foam::scalar Foam::swirlInjector::teoi() const
391     return massFlowRateProfile_.last()[0];
395 Foam::scalar Foam::swirlInjector::fractionOfInjection(const scalar time) const
397     return integrateTable(massFlowRateProfile_, time)/mass_;
401 Foam::scalar Foam::swirlInjector::injectedMass
403     const scalar t
404 ) const
406     return mass_*fractionOfInjection(t);
410 void Foam::swirlInjector::correctProfiles
412     const liquidMixtureProperties& fuel,
413     const scalar referencePressure
416     scalar A = 0.25*constant::mathematical::pi*sqr(d_);
417     scalar pDummy = 1.0e+5;
418     scalar rho = fuel.rho(pDummy, T_, X_);
420     forAll(velocityProfile_, i)
421     {
422         scalar Pinj = getTableValue
423         (
424             injectionPressureProfile_,
425             massFlowRateProfile_[i][0]
426         );
427         scalar mfr = massFlowRateProfile_[i][1]/(rho*A);
428         scalar v = sqrt(2.0*(Pinj - referencePressure)/rho);
429         velocityProfile_[i][1] = v;
430         CdProfile_[i][1] = mfr/v;
431     }
435 Foam::vector Foam::swirlInjector::tan1(const label n) const
437     return tangentialInjectionVector1_;
441 Foam::vector Foam::swirlInjector::tan2(const label n) const
443     return tangentialInjectionVector2_;
447 // ************************************************************************* //