FreeFOAM The Cross-Platform CFD Toolkit
ODEChemistryModel.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
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.
15 
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.
20 
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/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "ODEChemistryModel.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 template<class CompType, class ThermoType>
34 (
35  const fvMesh& mesh,
36  const word& compTypeName,
37  const word& thermoTypeName
38 )
39 :
40  CompType(mesh, thermoTypeName),
41 
42  ODE(),
43 
44  Y_(this->thermo().composition().Y()),
45 
46  reactions_
47  (
48  dynamic_cast<const reactingMixture<ThermoType>&>(this->thermo())
49  ),
50  specieThermo_
51  (
52  dynamic_cast<const reactingMixture<ThermoType>&>
53  (this->thermo()).speciesData()
54  ),
55 
56  nSpecie_(Y_.size()),
57  nReaction_(reactions_.size()),
58 
59  solver_
60  (
62  (
63  *this,
64  compTypeName,
65  thermoTypeName
66  )
67  ),
68 
69  RR_(nSpecie_)
70 {
71  // create the fields for the chemistry sources
72  forAll(RR_, fieldI)
73  {
74  RR_.set
75  (
76  fieldI,
77  new scalarField(mesh.nCells(), 0.0)
78  );
79  }
80 
81  Info<< "ODEChemistryModel: Number of species = " << nSpecie_
82  << " and reactions = " << nReaction_ << endl;
83 }
84 
85 
86 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
87 
88 template<class CompType, class ThermoType>
90 {}
91 
92 
93 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
94 
95 template<class CompType, class ThermoType>
97 (
98  const scalarField& c,
99  const scalar T,
100  const scalar p
101 ) const
102 {
103  scalar pf, cf, pr, cr;
104  label lRef, rRef;
105 
106  scalarField om(nEqns(), 0.0);
107 
108  forAll(reactions_, i)
109  {
110  const Reaction<ThermoType>& R = reactions_[i];
111 
112  scalar omegai = omega
113  (
114  R, c, T, p, pf, cf, lRef, pr, cr, rRef
115  );
116 
117  forAll(R.lhs(), s)
118  {
119  label si = R.lhs()[s].index;
120  scalar sl = R.lhs()[s].stoichCoeff;
121  om[si] -= sl*omegai;
122  }
123 
124  forAll(R.rhs(), s)
125  {
126  label si = R.rhs()[s].index;
127  scalar sr = R.rhs()[s].stoichCoeff;
128  om[si] += sr*omegai;
129  }
130  }
131 
132  return om;
133 }
134 
135 
136 template<class CompType, class ThermoType>
138 (
139  const Reaction<ThermoType>& R,
140  const scalarField& c,
141  const scalar T,
142  const scalar p,
143  scalar& pf,
144  scalar& cf,
145  label& lRef,
146  scalar& pr,
147  scalar& cr,
148  label& rRef
149 ) const
150 {
151  scalarField c2(nSpecie_, 0.0);
152  for (label i=0; i<nSpecie_; i++)
153  {
154  c2[i] = max(0.0, c[i]);
155  }
156 
157  scalar kf = R.kf(T, p, c2);
158  scalar kr = R.kr(kf, T, p, c2);
159 
160  pf = 1.0;
161  pr = 1.0;
162 
163  label Nl = R.lhs().size();
164  label Nr = R.rhs().size();
165 
166  label slRef = 0;
167  lRef = R.lhs()[slRef].index;
168 
169  pf = kf;
170  for (label s=1; s<Nl; s++)
171  {
172  label si = R.lhs()[s].index;
173 
174  if (c[si] < c[lRef])
175  {
176  scalar exp = R.lhs()[slRef].exponent;
177  pf *= pow(max(0.0, c[lRef]), exp);
178  lRef = si;
179  slRef = s;
180  }
181  else
182  {
183  scalar exp = R.lhs()[s].exponent;
184  pf *= pow(max(0.0, c[si]), exp);
185  }
186  }
187  cf = max(0.0, c[lRef]);
188 
189  {
190  scalar exp = R.lhs()[slRef].exponent;
191  if (exp<1.0)
192  {
193  if (cf > SMALL)
194  {
195  pf *= pow(cf, exp - 1.0);
196  }
197  else
198  {
199  pf = 0.0;
200  }
201  }
202  else
203  {
204  pf *= pow(cf, exp - 1.0);
205  }
206  }
207 
208  label srRef = 0;
209  rRef = R.rhs()[srRef].index;
210 
211  // find the matrix element and element position for the rhs
212  pr = kr;
213  for (label s=1; s<Nr; s++)
214  {
215  label si = R.rhs()[s].index;
216  if (c[si] < c[rRef])
217  {
218  scalar exp = R.rhs()[srRef].exponent;
219  pr *= pow(max(0.0, c[rRef]), exp);
220  rRef = si;
221  srRef = s;
222  }
223  else
224  {
225  scalar exp = R.rhs()[s].exponent;
226  pr *= pow(max(0.0, c[si]), exp);
227  }
228  }
229  cr = max(0.0, c[rRef]);
230 
231  {
232  scalar exp = R.rhs()[srRef].exponent;
233  if (exp<1.0)
234  {
235  if (cr>SMALL)
236  {
237  pr *= pow(cr, exp - 1.0);
238  }
239  else
240  {
241  pr = 0.0;
242  }
243  }
244  else
245  {
246  pr *= pow(cr, exp - 1.0);
247  }
248  }
249 
250  return pf*cf - pr*cr;
251 }
252 
253 
254 template<class CompType, class ThermoType>
256 (
257  const scalar time,
258  const scalarField &c,
259  scalarField& dcdt
260 ) const
261 {
262  scalar T = c[nSpecie_];
263  scalar p = c[nSpecie_ + 1];
264 
265  dcdt = omega(c, T, p);
266 
267  // constant pressure
268  // dT/dt = ...
269  scalar rho = 0.0;
270  scalar cSum = 0.0;
271  for (label i=0; i<nSpecie_; i++)
272  {
273  scalar W = specieThermo_[i].W();
274  cSum += c[i];
275  rho += W*c[i];
276  }
277  scalar mw = rho/cSum;
278  scalar cp = 0.0;
279  for (label i=0; i<nSpecie_; i++)
280  {
281  scalar cpi = specieThermo_[i].cp(T);
282  scalar Xi = c[i]/rho;
283  cp += Xi*cpi;
284  }
285  cp /= mw;
286 
287  scalar dT = 0.0;
288  for (label i=0; i<nSpecie_; i++)
289  {
290  scalar hi = specieThermo_[i].h(T);
291  dT += hi*dcdt[i];
292  }
293  dT /= rho*cp;
294 
295  // limit the time-derivative, this is more stable for the ODE
296  // solver when calculating the allowed time step
297  scalar dtMag = min(500.0, mag(dT));
298  dcdt[nSpecie_] = -dT*dtMag/(mag(dT) + 1.0e-10);
299 
300  // dp/dt = ...
301  dcdt[nSpecie_+1] = 0.0;
302 }
303 
304 
305 template<class CompType, class ThermoType>
307 (
308  const scalar t,
309  const scalarField& c,
310  scalarField& dcdt,
311  scalarSquareMatrix& dfdc
312 ) const
313 {
314  scalar T = c[nSpecie_];
315  scalar p = c[nSpecie_ + 1];
316 
317  scalarField c2(nSpecie_, 0.0);
318  for (label i=0; i<nSpecie_; i++)
319  {
320  c2[i] = max(c[i], 0.0);
321  }
322 
323  for (label i=0; i<nEqns(); i++)
324  {
325  for (label j=0; j<nEqns(); j++)
326  {
327  dfdc[i][j] = 0.0;
328  }
329  }
330 
331  // length of the first argument must be nSpecie()
332  dcdt = omega(c2, T, p);
333 
334  for (label ri=0; ri<reactions_.size(); ri++)
335  {
336  const Reaction<ThermoType>& R = reactions_[ri];
337 
338  scalar kf0 = R.kf(T, p, c2);
339  scalar kr0 = R.kr(T, p, c2);
340 
341  forAll(R.lhs(), j)
342  {
343  label sj = R.lhs()[j].index;
344  scalar kf = kf0;
345  forAll(R.lhs(), i)
346  {
347  label si = R.lhs()[i].index;
348  scalar el = R.lhs()[i].exponent;
349  if (i == j)
350  {
351  if (el < 1.0)
352  {
353  if (c2[si]>SMALL)
354  {
355  kf *= el*pow(c2[si] + VSMALL, el - 1.0);
356  }
357  else
358  {
359  kf = 0.0;
360  }
361  }
362  else
363  {
364  kf *= el*pow(c2[si], el - 1.0);
365  }
366  }
367  else
368  {
369  kf *= pow(c2[si], el);
370  }
371  }
372 
373  forAll(R.lhs(), i)
374  {
375  label si = R.lhs()[i].index;
376  scalar sl = R.lhs()[i].stoichCoeff;
377  dfdc[si][sj] -= sl*kf;
378  }
379  forAll(R.rhs(), i)
380  {
381  label si = R.rhs()[i].index;
382  scalar sr = R.rhs()[i].stoichCoeff;
383  dfdc[si][sj] += sr*kf;
384  }
385  }
386 
387  forAll(R.rhs(), j)
388  {
389  label sj = R.rhs()[j].index;
390  scalar kr = kr0;
391  forAll(R.rhs(), i)
392  {
393  label si = R.rhs()[i].index;
394  scalar er = R.rhs()[i].exponent;
395  if (i==j)
396  {
397  if (er<1.0)
398  {
399  if (c2[si]>SMALL)
400  {
401  kr *= er*pow(c2[si] + VSMALL, er - 1.0);
402  }
403  else
404  {
405  kr = 0.0;
406  }
407  }
408  else
409  {
410  kr *= er*pow(c2[si], er - 1.0);
411  }
412  }
413  else
414  {
415  kr *= pow(c2[si], er);
416  }
417  }
418 
419  forAll(R.lhs(), i)
420  {
421  label si = R.lhs()[i].index;
422  scalar sl = R.lhs()[i].stoichCoeff;
423  dfdc[si][sj] += sl*kr;
424  }
425  forAll(R.rhs(), i)
426  {
427  label si = R.rhs()[i].index;
428  scalar sr = R.rhs()[i].stoichCoeff;
429  dfdc[si][sj] -= sr*kr;
430  }
431  }
432  }
433 
434  // calculate the dcdT elements numerically
435  scalar delta = 1.0e-8;
436  scalarField dcdT0 = omega(c2, T - delta, p);
437  scalarField dcdT1 = omega(c2, T + delta, p);
438 
439  for (label i=0; i<nEqns(); i++)
440  {
441  dfdc[i][nSpecie()] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
442  }
443 
444 }
445 
446 
447 template<class CompType, class ThermoType>
450 {
451  scalar pf, cf, pr, cr;
452  label lRef, rRef;
453 
454  const volScalarField rho
455  (
456  IOobject
457  (
458  "rho",
459  this->time().timeName(),
460  this->mesh(),
461  IOobject::NO_READ,
462  IOobject::NO_WRITE,
463  false
464  ),
465  this->thermo().rho()
466  );
467 
469  (
470  new volScalarField
471  (
472  IOobject
473  (
474  "tc",
475  this->time().timeName(),
476  this->mesh(),
477  IOobject::NO_READ,
478  IOobject::NO_WRITE
479  ),
480  this->mesh(),
481  dimensionedScalar("zero", dimTime, SMALL),
482  zeroGradientFvPatchScalarField::typeName
483  )
484  );
485 
486  scalarField& tc = ttc();
487 
488  label nReaction = reactions_.size();
489 
490  if (this->chemistry_)
491  {
492  forAll(rho, celli)
493  {
494  scalar rhoi = rho[celli];
495  scalar Ti = this->thermo().T()[celli];
496  scalar pi = this->thermo().p()[celli];
497  scalarField c(nSpecie_);
498  scalar cSum = 0.0;
499 
500  for (label i=0; i<nSpecie_; i++)
501  {
502  scalar Yi = Y_[i][celli];
503  c[i] = rhoi*Yi/specieThermo_[i].W();
504  cSum += c[i];
505  }
506 
507  forAll(reactions_, i)
508  {
509  const Reaction<ThermoType>& R = reactions_[i];
510 
511  omega
512  (
513  R, c, Ti, pi, pf, cf, lRef, pr, cr, rRef
514  );
515 
516  forAll(R.rhs(), s)
517  {
518  scalar sr = R.rhs()[s].stoichCoeff;
519  tc[celli] += sr*pf*cf;
520  }
521  }
522  tc[celli] = nReaction*cSum/tc[celli];
523  }
524  }
525 
526 
527  ttc().correctBoundaryConditions();
528 
529  return ttc;
530 }
531 
532 
533 template<class CompType, class ThermoType>
536 {
538  (
539  new volScalarField
540  (
541  IOobject
542  (
543  "Sh",
544  this->mesh_.time().timeName(),
545  this->mesh_,
546  IOobject::NO_READ,
547  IOobject::NO_WRITE,
548  false
549  ),
550  this->mesh_,
552  zeroGradientFvPatchScalarField::typeName
553  )
554  );
555 
556  if (this->chemistry_)
557  {
558  scalarField& Sh = tSh();
559 
560  forAll(Y_, i)
561  {
562  forAll(Sh, cellI)
563  {
564  scalar hi = specieThermo_[i].Hc();
565  Sh[cellI] -= hi*RR_[i][cellI];
566  }
567  }
568  }
569 
570  return tSh;
571 }
572 
573 
574 template<class CompType, class ThermoType>
577 {
579  (
580  new volScalarField
581  (
582  IOobject
583  (
584  "dQ",
585  this->mesh_.time().timeName(),
586  this->mesh_,
587  IOobject::NO_READ,
588  IOobject::NO_WRITE,
589  false
590  ),
591  this->mesh_,
592  dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
593  zeroGradientFvPatchScalarField::typeName
594  )
595  );
596 
597  if (this->chemistry_)
598  {
599  volScalarField& dQ = tdQ();
600  dQ.dimensionedInternalField() = this->mesh_.V()*Sh()();
601  }
602 
603  return tdQ;
604 }
605 
606 
607 template<class CompType, class ThermoType>
609 {
610  // nEqns = number of species + temperature + pressure
611  return nSpecie_ + 2;
612 }
613 
614 
615 template<class CompType, class ThermoType>
617 {
618  const volScalarField rho
619  (
620  IOobject
621  (
622  "rho",
623  this->time().timeName(),
624  this->mesh(),
625  IOobject::NO_READ,
626  IOobject::NO_WRITE,
627  false
628  ),
629  this->thermo().rho()
630  );
631 
632  for (label i=0; i<nSpecie_; i++)
633  {
634  RR_[i].setSize(rho.size());
635  }
636 
637  if (this->chemistry_)
638  {
639  forAll(rho, celli)
640  {
641  for (label i=0; i<nSpecie_; i++)
642  {
643  RR_[i][celli] = 0.0;
644  }
645 
646  scalar rhoi = rho[celli];
647  scalar Ti = this->thermo().T()[celli];
648  scalar pi = this->thermo().p()[celli];
649 
650  scalarField c(nSpecie_);
651  scalarField dcdt(nEqns(), 0.0);
652 
653  for (label i=0; i<nSpecie_; i++)
654  {
655  scalar Yi = Y_[i][celli];
656  c[i] = rhoi*Yi/specieThermo_[i].W();
657  }
658 
659  dcdt = omega(c, Ti, pi);
660 
661  for (label i=0; i<nSpecie_; i++)
662  {
663  RR_[i][celli] = dcdt[i]*specieThermo_[i].W();
664  }
665  }
666  }
667 }
668 
669 
670 template<class CompType, class ThermoType>
672 (
673  const scalar t0,
674  const scalar deltaT
675 )
676 {
677  const volScalarField rho
678  (
679  IOobject
680  (
681  "rho",
682  this->time().timeName(),
683  this->mesh(),
684  IOobject::NO_READ,
685  IOobject::NO_WRITE,
686  false
687  ),
688  this->thermo().rho()
689  );
690 
691  for (label i=0; i<nSpecie_; i++)
692  {
693  RR_[i].setSize(rho.size());
694  }
695 
696  if (!this->chemistry_)
697  {
698  return GREAT;
699  }
700 
701  scalar deltaTMin = GREAT;
702 
703  tmp<volScalarField> thc = this->thermo().hc();
704  const scalarField& hc = thc();
705 
706  forAll(rho, celli)
707  {
708  for (label i=0; i<nSpecie_; i++)
709  {
710  RR_[i][celli] = 0.0;
711  }
712 
713  scalar rhoi = rho[celli];
714  scalar Ti = this->thermo().T()[celli];
715  scalar hi = this->thermo().hs()[celli] + hc[celli];
716  scalar pi = this->thermo().p()[celli];
717 
718  scalarField c(nSpecie_);
719  scalarField c0(nSpecie_);
720  scalarField dc(nSpecie_, 0.0);
721 
722  for (label i=0; i<nSpecie_; i++)
723  {
724  c[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
725  }
726  c0 = c;
727 
728  scalar t = t0;
729  scalar tauC = this->deltaTChem_[celli];
730  scalar dt = min(deltaT, tauC);
731  scalar timeLeft = deltaT;
732 
733  // calculate the chemical source terms
734  scalar cTot = 0.0;
735 
736  while (timeLeft > SMALL)
737  {
738  tauC = solver().solve(c, Ti, pi, t, dt);
739  t += dt;
740 
741  // update the temperature
742  cTot = sum(c);
743  ThermoType mixture(0.0*specieThermo_[0]);
744  for (label i=0; i<nSpecie_; i++)
745  {
746  mixture += (c[i]/cTot)*specieThermo_[i];
747  }
748  Ti = mixture.TH(hi, Ti);
749 
750  timeLeft -= dt;
751  this->deltaTChem_[celli] = tauC;
752  dt = min(timeLeft, tauC);
753  dt = max(dt, SMALL);
754  }
755  deltaTMin = min(tauC, deltaTMin);
756 
757  dc = c - c0;
758  scalar WTot = 0.0;
759  for (label i=0; i<nSpecie_; i++)
760  {
761  WTot += c[i]*specieThermo_[i].W();
762  }
763  WTot /= cTot;
764 
765  for (label i=0; i<nSpecie_; i++)
766  {
767  RR_[i][celli] = dc[i]*specieThermo_[i].W()/deltaT;
768  }
769  }
770 
771  // Don't allow the time-step to change more than a factor of 2
772  deltaTMin = min(deltaTMin, 2*deltaT);
773 
774  return deltaTMin;
775 }
776 
777 
778 // ************************ vim: set sw=4 sts=4 et: ************************ //