ODT
Loading...
Searching...
No Matches
eddy.cc
Go to the documentation of this file.
1
6#include "eddy.h"
7#include "domain.h"
8#include <cmath>
9#include <numeric> //accumulate
10#include <algorithm> // min_element
11
13
19void eddy::init(domain *p_domn, domain *p_eddl) {
20 domn = p_domn;
21 eddl = p_eddl;
22
23 esdp1 = -2.0*domn->pram->Lp;
24 esdp2 = exp(-2.0*domn->pram->Lp/domn->pram->Lmax) - exp(-2.0*domn->pram->Lp/domn->pram->Lmin);
25 esdp3 = exp(-2.0*domn->pram->Lp/domn->pram->Lmin);
26 esdp4 = -esdp1 / esdp2;
27
28 cCoef = vector<double> (3,0.0);
29 bCoef = vector<double> (3,0.0);
30
31 LperiodicEddy = false;
32
33 cca = vector<double>(7);
34 ccb = vector<double>(5);
35 ccc = vector<double>(5);
36 ccd = vector<double>(5);
37
38 cca[0] = -6.58411371e+02;
39 cca[1] = 1.19518310e+03;
40 cca[2] = -8.21169900e+02;
41 cca[3] = 2.61583310e+02;
42 cca[4] = -3.38649746e+01;
43 cca[5] = -6.11440937e-01;
44 cca[6] = 1.01054351e+00;
45
46 ccb[0] = -0.48987445;
47 ccb[1] = 2.39482488;
48 ccb[2] = -4.48286508;
49 ccb[3] = 3.89924937;
50 ccb[4] = -0.40448995;
51
52 ccc[0] = -0.00169718;
53 ccc[1] = 0.02173516;
54 ccc[2] = -0.10620329;
55 ccc[3] = 0.2399716;
56 ccc[4] = 0.77720962;
57
58 ccd[0] = -1.98503860e-07;
59 ccd[1] = 1.31972798e-05;
60 ccd[2] = -3.21487789e-04;
61 ccd[3] = 3.44046508e-03;
62 ccd[4] = 9.85964172e-01;
63
64}
65
67
71void eddy::sampleEddySize() {
72
73 eddySize = esdp1 / log( domn->rand->getRand() * esdp2 + esdp3 );
74
75}
76
78
86void eddy::sampleEddyPosition() {
87
88 if(!domn->pram->Lperiodic) {
89 leftEdge = domn->rand->getRand() * (domn->Ldomain() - eddySize) + domn->posf->d.at(0);
90 rightEdge = leftEdge + eddySize;
91 if(leftEdge < domn->posf->d.at(0)) leftEdge = domn->posf->d.at(0);
92 if(rightEdge > domn->posf->d.at(domn->ngrd)) rightEdge = domn->posf->d.at(domn->ngrd);
93 }
94 else {
95 LperiodicEddy = false; // reset the value
96 leftEdge = domn->rand->getRand() * domn->Ldomain() + domn->posf->d.at(0);
97 rightEdge = leftEdge + eddySize;
98 if(rightEdge > (domn->Ldomain() + domn->posf->d.at(0)))
99 LperiodicEddy = true;
100 }
101
102}
103
105
116void eddy::tripMap(domain *line, const int iS, int iE, const double C, const bool LsplitAtEddy) {
117
118 int negrd = iE-iS+1;
119 int negrd2 = 2*negrd;
120 int negrd3 = 3*negrd;
121
122 double fracVleft;
123 double fracVmidl;
124 double fracVrght;
125 double pm1; // +1 or -1 factor to change sign
126 double invC = 1.0/C;
127
128 int k, i, ip, im, ie, iw;
129
130 //--------- Set fracVlft, fracVmid, fracVrgt
131
132 double r1 = (2*leftEdge+rightEdge)/3;
133 double r2 = (leftEdge+2*rightEdge)/3;
134
135 pm1 =leftEdge*r1 < 0 ? -1.0 : 1.0; // handles cells that split x=0
136 double Vleft = abs(pow(abs(leftEdge), C) - pm1*pow(abs(r1), C));
137
138 pm1 =r1*r2 < 0 ? -1.0 : 1.0; // handles cells that split x=0
139 double Vmidl = abs(pow(abs(r1), C) - pm1*pow(abs(r2), C));
140
141 pm1 =r2*rightEdge < 0 ? -1.0 : 1.0; // handles cells that split x=0
142 double Vrght = abs(pow(abs(r2), C) - pm1*pow(abs(rightEdge), C));
143
144 double Vtot = Vleft + Vmidl + Vrght;
145 fracVleft = Vleft / Vtot;
146 fracVmidl = Vmidl / Vtot;
147 fracVrght = Vrght / Vtot;
148
149 if(domn->pram->LTMA) { // use three equal volume segments instead of three equal length sements
150 fracVleft = 1.0/3; // no difference for planar, but affects cylindrical or spherical
151 fracVmidl = fracVleft;
152 fracVrght = fracVleft;
153 }
154
155 //---------
156
157 pos0 = line->pos->d; // for filling kernel
158 pos0.at(0) = 0.5*(line->posf->d.at(1)+leftEdge);
159 pos0.at(line->ngrd-1) = 0.5*(line->posf->d.at(line->ngrd-1)+rightEdge);
160
161 //----------- Grab cell "volume" array (delta(x^c)) in the eddy region
162
163 dxc.resize(negrd,0.0);
164
165 i=0; ip=1;
166 pm1 = line->posf->d.at(iS+ip)*leftEdge < 0 ? -1.0 : 1.0; // handles cells that split x=0
167 dxc.at(i) = abs(pow(abs(line->posf->d.at(iS+ip)), C) - pm1*pow(abs(leftEdge), C));
168
169 for(int i=1, ip=2; i<negrd-1; i++, ip++) {
170 pm1 = line->posf->d.at(iS+ip)*line->posf->d.at(iS+i) < 0 ? -1.0 : 1.0; // handles cells that split x=0
171 dxc.at(i) = abs(pow(abs(line->posf->d.at(iS+ip)), C) - pm1*pow(abs(line->posf->d.at(iS+i) ), C));
172 }
173
174 i=negrd-1;
175 pm1 = line->posf->d.at(iS+i)*rightEdge < 0 ? -1.0 : 1.0; // handles cells that split x=0
176 dxc.at(i) = abs(pow(abs(rightEdge), C) - pm1*pow(abs(line->posf->d.at(iS+i)), C));
177
178 //---------- insert cells to be filled with values
179
180 for(k=0; k<line->v.size(); k++)
181 line->v.at(k)->d.insert(line->v.at(k)->d.begin()+iE+1, negrd2, 0.0);
182
183 line->ngrd += negrd2;
184 line->ngrdf = line->ngrd+1;
185 iE += negrd2;
186
187 //---------- fill in variable data (except pos, posf)
188
189 for(k=0; k<line->v.size(); k++) {
190 if(line->v.at(k)->var_name=="pos" || line->v.at(k)->var_name=="posf") continue;
191 for(i=0; i<negrd; i++) {
192 line->v.at(k)->d.at(iS+negrd2+i) = line->v.at(k)->d.at(iS+i);
193 line->v.at(k)->d.at(iS+negrd2-1-i) = line->v.at(k)->d.at(iS+i);
194 }
195 }
196
197 //---------- update face positions
198
199 if(C==1 || leftEdge > 0.0) { // case: planar, or eddy on the right of the centerline
200
201 // posf.at(iS) is the same
202 line->posf->d.at(iS+1) = pow(pow(leftEdge,C) + fracVleft*dxc.at(0), invC);
203 for(ie=2, iw=1, i=1; ie<=negrd; ie++, iw++, i++)
204 line->posf->d.at(iS+ie) = pow(pow(line->posf->d.at(iS+iw),C) + fracVleft*dxc.at(i), invC);
205 for(ie=negrd+1,iw=negrd,i=negrd-1; ie<=negrd2; ie++, iw++, i--)
206 line->posf->d.at(iS+ie) = pow(pow(line->posf->d.at(iS+iw),C) + fracVmidl*dxc.at(i), invC);
207 for(ie=negrd2+1,iw=negrd2,i=0; ie<negrd3; ie++, iw++, i++)
208 line->posf->d.at(iS+ie) = pow(pow(line->posf->d.at(iS+iw),C) + fracVrght*dxc.at(i), invC);
209 // line->posf->d.at(iS+negrd3) is the same due to the inserted cells before.
210
211 }
212 else if(rightEdge < 0.0) { // case: eddy on the left of centerline
213
214 // posf.at(iS) is the same
215 line->posf->d.at(iS+1) = -pow(pow(abs(leftEdge),C) - fracVleft*dxc.at(0), invC);
216 for(ie=2, iw=1, i=1; ie<=negrd; ie++, iw++, i++)
217 line->posf->d.at(iS+ie) = -pow(pow(abs(line->posf->d.at(iS+iw)),C) - fracVleft*dxc.at(i), invC);
218 for(ie=negrd+1,iw=negrd,i=negrd-1; ie<=negrd2; ie++, iw++, i--)
219 line->posf->d.at(iS+ie) = -pow(pow(abs(line->posf->d.at(iS+iw)),C) - fracVmidl*dxc.at(i), invC);
220 for(ie=negrd2+1,iw=negrd2,i=0; ie<negrd3; ie++, iw++, i++)
221 line->posf->d.at(iS+ie) = -pow(pow(abs(line->posf->d.at(iS+iw)),C) - fracVrght*dxc.at(i), invC);
222 // line->posf->d.at(iS+negrd3) is the same due to the inserted cells before.
223
224 }
225 else { // case: eddy splits the centerline
226 // Note, this can directly replace the above two cases!
227
228 double dmb;
229
230 // posf.at(iS) is the same
231 pm1 = leftEdge < 0 ? -1 : 1;
232 dmb = pow(abs(leftEdge),C) + pm1*fracVleft*dxc.at(0);
233 line->posf->d.at(iS+1) = dmb < 0.0 ? pow(-dmb, invC) : pm1*pow(dmb, invC);
234 for(ie=2, iw=1, i=1; ie<=negrd; ie++, iw++, i++) {
235 pm1 = line->posf->d.at(iS+iw) < 0 ? -1 : 1;
236 dmb = pow(abs(line->posf->d.at(iS+iw)),C) + pm1*fracVleft*dxc.at(i);
237 line->posf->d.at(iS+ie) = dmb < 0.0 ? pow(-dmb, invC) : pm1*pow(dmb, invC);
238 }
239 for(ie=negrd+1,iw=negrd,i=negrd-1; ie<=negrd2; ie++, iw++, i--) {
240 pm1 = line->posf->d.at(iS+iw) < 0 ? -1 : 1;
241 dmb = pow(abs(line->posf->d.at(iS+iw)),C) + pm1*fracVmidl*dxc.at(i);
242 line->posf->d.at(iS+ie) = dmb < 0.0 ? pow(-dmb, invC) : pm1*pow(dmb, invC);
243 }
244 for(ie=negrd2+1,iw=negrd2,i=0; ie<negrd3; ie++, iw++, i++) {
245 pm1 = line->posf->d.at(iS+iw) < 0 ? -1 : 1;
246 dmb = pow(abs(line->posf->d.at(iS+iw)),C) + pm1*fracVrght*dxc.at(i);
247 line->posf->d.at(iS+ie) = dmb < 0.0 ? pow(-dmb, invC) : pm1*pow(dmb, invC);
248 }
249 // line->posf->d.at(iS+negrd3) is the same due to the inserted cells before.
250 }
251
252 //---------- update center positions
253
254 line->pos->setVar();
255
256 //---------- split at eddy faces if needed
257
258 if(LsplitAtEddy) {
259
260 vector<double> interPos(3); // used for cell split
261
262 if(abs(line->posf->d.at(iE+1) - rightEdge) > 1.0e-15) {
263 interPos.at(0) = line->posf->d.at(iE);
264 interPos.at(1) = rightEdge;
265 interPos.at(2) = line->posf->d.at(iE+1);
266 domn->mesher->splitCell(iE, 1, interPos); // don't interpolate here
267 }
268 if(abs(line->posf->d.at(iS) - leftEdge) > 1.0e-15) {
269 interPos.at(0) = line->posf->d.at(iS);
270 interPos.at(1) = leftEdge;
271 interPos.at(2) = line->posf->d.at(iS+1);
272 domn->mesher->splitCell(iS, 1, interPos); // again, don't interp
273 }
274 }
275
276}
277
279
285bool eddy::eddyTau(const double Z_value, const double C) {
286
287 double KK=0; // equivalent of the 4/27 fac
288 double rhoK=0, rhoJ=0, rhoKK=0, rhoJK=0;
289 vector<double> uRhoJ(3,0.0), uRhoK(3,0.0);
290
291 vector<double> P(3);
292 vector<double> Q(3);
293 double S;
294 double A;
295 double eKinEddy = 0.0;
296 double ePeEddy = 0.0;
297 double eViscPenalty = 0.0;
298 double eDL = 0.0;
299 double rho_ref = 0.0;
300
301 double rhoEddy=0.0, viscEddy=0.0;
302
303 vector<double> intRhoKi(eddl->ngrd);
304 vector<double> intRhoJi(eddl->ngrd);
305
306 int i;
307 double Etot;
308
309 //----------- set dxc
310
311 domn->mesher->setGridDxc(eddl, dxc, C);
312
313 if(leftEdge >= 0.0 || eddl->posf->d.at(1) <= 0.0)
314 dxc[0] = abs( pow(abs(eddl->posf->d.at(1)), C) - pow(abs(leftEdge),C) );
315 else
316 dxc[0] = abs( pow(abs(eddl->posf->d.at(1)), C) + pow(abs(leftEdge),C) );
317
318 if(eddl->posf->d.at(eddl->ngrd-1) >= 0.0 || rightEdge <= 0.0)
319 dxc[eddl->ngrd-1] = abs( pow(abs(eddl->posf->d.at(eddl->ngrd-1)), C) - pow(abs(rightEdge),C) );
320 else
321 dxc[eddl->ngrd-1] = abs( pow(abs(eddl->posf->d.at(eddl->ngrd-1)), C) + pow(abs(rightEdge),C) );
322
323 double VolE = accumulate(dxc.begin(), dxc.end(), 0.0);
324
325 //-----------
326
327 if(domn->pram->LplanarTau)
328 fillKernel_planarAnalytic();
329 else
330 fillKernel();
331
332 if(domn->pram->LdoDL) { // Darrieus Landau instability for variable density flows
333 vector<double> tempVec(eddl->ngrd);
334 transform(eddl->rho->d.begin(), eddl->rho->d.end(), dxc.begin(), tempVec.begin(), multiplies<double>());
335 rho_ref = accumulate(tempVec.begin(), tempVec.end(), 0.0) / VolE;
336 }
337
339
340 //---------- rhoK, rhoJ, UrhoK, UrhoJ
341
342 for(i=0; i<eddl->ngrd; i++) {
343 intRhoKi.at(i) = K.at(i)*eddl->rho->d.at(i)*dxc.at(i);
344 intRhoJi.at(i) = abs(intRhoKi.at(i));
345 KK += K.at(i)*K.at(i)*dxc.at(i);
346 if(!domn->pram->Lspatial) {
347 rhoK += intRhoKi.at(i);
348 rhoJ += abs(intRhoKi.at(i));
349 rhoKK += K.at(i)*intRhoKi.at(i);
350 rhoJK += K.at(i)*intRhoJi.at(i);
351 uRhoK.at(0) += intRhoKi.at(i)*eddl->uvel->d.at(i);
352 uRhoK.at(1) += intRhoKi.at(i)*eddl->vvel->d.at(i);
353 uRhoK.at(2) += intRhoKi.at(i)*eddl->wvel->d.at(i);
354 uRhoJ.at(0) += intRhoJi.at(i)*eddl->uvel->d.at(i);
355 uRhoJ.at(1) += intRhoJi.at(i)*eddl->vvel->d.at(i);
356 uRhoJ.at(2) += intRhoJi.at(i)*eddl->wvel->d.at(i);
357 if(domn->pram->LdoDL)
358 eDL += eddl->aDL->d.at(i)*K[i]*(eddl->rho->d[i] - rho_ref ) * dxc[i];
359 }
360 else { // mass flux, not mass
361 rhoK += eddl->uvel->d.at(i) * intRhoKi.at(i);
362 rhoJ += eddl->uvel->d.at(i) * abs(intRhoKi.at(i));
363 rhoKK += eddl->uvel->d.at(i) * K.at(i)*intRhoKi.at(i);
364 rhoJK += eddl->uvel->d.at(i) * K.at(i)*intRhoJi.at(i);
365 uRhoK.at(0) += eddl->uvel->d.at(i) * intRhoKi.at(i)*eddl->uvel->d.at(i);
366 uRhoK.at(1) += eddl->uvel->d.at(i) * intRhoKi.at(i)*eddl->vvel->d.at(i);
367 uRhoK.at(2) += eddl->uvel->d.at(i) * intRhoKi.at(i)*eddl->wvel->d.at(i);
368 uRhoJ.at(0) += eddl->uvel->d.at(i) * intRhoJi.at(i)*eddl->uvel->d.at(i);
369 uRhoJ.at(1) += eddl->uvel->d.at(i) * intRhoJi.at(i)*eddl->vvel->d.at(i);
370 uRhoJ.at(2) += eddl->uvel->d.at(i) * intRhoJi.at(i)*eddl->wvel->d.at(i);
371 if(domn->pram->LdoDL)
372 eDL += eddl->uvel->d[i] * eddl->aDL->d.at(i)*K[i]*(eddl->rho->d[i] - rho_ref ) * dxc[i];
373 }
374 }
375
377
378 A = rhoK/rhoJ;
379 S = 0.5*(A*A+1.0)*rhoKK - A*rhoJK;
380 for(i=0; i<3; i++) {
381 P.at(i) = uRhoK.at(i) - A * uRhoJ.at(i);
382 Q.at(i) = 0.25*P.at(i)*P.at(i)/S;
383 }
384
385 eKinEddy = KK/(eddySize*eddySize*VolE) * (Q.at(0)+Q.at(1)+Q.at(2));
386
387 ePeEddy = (domn->pram->LPeEddy ? domn->pram->g*rhoK : 0.0);
388
390
391 for(i=0; i<eddl->ngrd; i++) {
392 rhoEddy += eddl->rho->d.at(i) *dxc.at(i);
393 viscEddy += eddl->dvisc->d.at(i)*dxc.at(i);
394 }
395 rhoEddy /= VolE;
396 viscEddy /= VolE;
397
398 eViscPenalty = Z_value*0.5 * VolE/(eddySize*eddySize)*viscEddy*viscEddy/rhoEddy;
399
400 double Ufavre;
401 if(domn->pram->Lspatial) {
402 Ufavre = eddyFavreAvgVelocity(dxc);
403 eViscPenalty *= Ufavre;
404 }
405
406 invTauEddy = 0.0;
407
408 Etot = eKinEddy - eViscPenalty + eDL + ePeEddy;
409
410 if(Etot < 0.0) return false;
411
412 invTauEddy = sqrt(2.0*KK/(rhoKK*VolE*eddySize*eddySize) * Etot);
413
414 if(domn->pram->Lspatial)
415 invTauEddy = invTauEddy/Ufavre; // 1/s --> 1/m
416
417 return true;
418}
419
421
430void eddy::set_kernel_coefficients() {
431
432 double rhoK=0, rhoJ=0, rhoKK=0, rhoJK=0;
433 vector<double> uRhoJ(3,0.0), uRhoK(3,0.0);
434
435 vector<double> P(3);
436 vector<double> Q(3);
437 double S;
438 double A;
439
440 vector<double> intRhoKi(eddl->ngrd);
441 vector<double> intRhoJi(eddl->ngrd);
442
443 int i;
444
445 //----------- set dxc
446
447 double C = domn->pram->cCoord;
448 domn->mesher->setGridDxc(eddl, dxc, C);
449
450 if(leftEdge >= 0.0 || eddl->posf->d.at(1) <= 0.0)
451 dxc[0] = abs( pow(abs(eddl->posf->d.at(1)), C) - pow(abs(leftEdge),C) );
452 else
453 dxc[0] = abs( pow(abs(eddl->posf->d.at(1)), C) + pow(abs(leftEdge),C) );
454
455 if(eddl->posf->d.at(eddl->ngrd-1) >= 0.0 || rightEdge <= 0.0)
456 dxc[eddl->ngrd-1] = abs( pow(abs(eddl->posf->d.at(eddl->ngrd-1)), C) - pow(abs(rightEdge),C) );
457 else
458 dxc[eddl->ngrd-1] = abs( pow(abs(eddl->posf->d.at(eddl->ngrd-1)), C) + pow(abs(rightEdge),C) );
459
460 //-----------
461
462 if(domn->pram->LplanarTau) // todo: not strictly needed if using the same kernel as in eddyTau
463 fillKernel_planarAnalytic(); // todo: probably is needed if using testThirds large eddy suppression
464
466
467 //---------- rhoK, rhoJ, UrhoK, UrhoJ
468
469 for(i=0; i<eddl->ngrd; i++) {
470 intRhoKi.at(i) = K.at(i)*eddl->rho->d.at(i)*dxc.at(i);
471 intRhoJi.at(i) = abs(intRhoKi.at(i));
472 if(!domn->pram->Lspatial) {
473 rhoK += intRhoKi.at(i);
474 rhoJ += abs(intRhoKi.at(i));
475 rhoKK += K.at(i)*intRhoKi.at(i);
476 rhoJK += K.at(i)*intRhoJi.at(i);
477 uRhoK.at(0) += intRhoKi.at(i)*eddl->uvel->d.at(i);
478 uRhoK.at(1) += intRhoKi.at(i)*eddl->vvel->d.at(i);
479 uRhoK.at(2) += intRhoKi.at(i)*eddl->wvel->d.at(i);
480 uRhoJ.at(0) += intRhoJi.at(i)*eddl->uvel->d.at(i);
481 uRhoJ.at(1) += intRhoJi.at(i)*eddl->vvel->d.at(i);
482 uRhoJ.at(2) += intRhoJi.at(i)*eddl->wvel->d.at(i);
483 }
484 else { // mass flux, not mass
485 rhoK += eddl->uvel->d.at(i) * intRhoKi.at(i);
486 rhoJ += eddl->uvel->d.at(i) * abs(intRhoKi.at(i));
487 rhoKK += eddl->uvel->d.at(i) * K.at(i)*intRhoKi.at(i);
488 rhoJK += eddl->uvel->d.at(i) * K.at(i)*intRhoJi.at(i);
489 uRhoK.at(0) += eddl->uvel->d.at(i) * intRhoKi.at(i)*eddl->uvel->d.at(i);
490 uRhoK.at(1) += eddl->uvel->d.at(i) * intRhoKi.at(i)*eddl->vvel->d.at(i);
491 uRhoK.at(2) += eddl->uvel->d.at(i) * intRhoKi.at(i)*eddl->wvel->d.at(i);
492 uRhoJ.at(0) += eddl->uvel->d.at(i) * intRhoJi.at(i)*eddl->uvel->d.at(i);
493 uRhoJ.at(1) += eddl->uvel->d.at(i) * intRhoJi.at(i)*eddl->vvel->d.at(i);
494 uRhoJ.at(2) += eddl->uvel->d.at(i) * intRhoJi.at(i)*eddl->wvel->d.at(i);
495 }
496 }
497
499
500 A = rhoK/rhoJ;
501 S = 0.5*(A*A+1.0)*rhoKK - A*rhoJK;
502 for(i=0; i<3; i++) {
503 P.at(i) = uRhoK.at(i) - A * uRhoJ.at(i);
504 Q.at(i) = 0.25*P.at(i)*P.at(i)/S;
505 }
506
508
509 cCoef.at(0) = 0.5/S * (-P.at(0) + (P.at(0)>0 ? 1.0 : -1.0)
510 * sqrt( (1-domn->pram->A_param)*P.at(0)*P.at(0)
511 + 0.5*domn->pram->A_param*(P.at(1)*P.at(1)+P.at(2)*P.at(2)) ));
512 cCoef.at(1) = 0.5/S * (-P.at(1) + (P.at(1)>0 ? 1.0 : -1.0)
513 * sqrt( (1-domn->pram->A_param)*P.at(1)*P.at(1)
514 + 0.5*domn->pram->A_param*(P.at(0)*P.at(0)+P.at(2)*P.at(2)) ));
515 cCoef.at(2) = 0.5/S * (-P.at(2) + (P.at(2)>0 ? 1.0 : -1.0)
516 * sqrt( (1-domn->pram->A_param)*P.at(2)*P.at(2)
517 + 0.5*domn->pram->A_param*(P.at(0)*P.at(0)+P.at(1)*P.at(1)) ));
518
519 for(i=0; i<3; i++)
520 bCoef.at(i) = -A*cCoef.at(i);
521
522}
523
525
528void eddy::computeEddyAcceptanceProb(const double dtSample) {
529
530 double f, g;
531
532 // f = esdp4/(eddySize*eddySize) * exp(esdp1/eddySize);
533 // Pa = dtSample * invTauEddy * domn->pram->C_param / (eddySize * eddySize * f * g);
534 // note cancellation of eddySize * eddySize
535
536 f = esdp4 * exp(esdp1/eddySize);
537
538 if(!domn->pram->Lperiodic)
539 g = 1.0/(domn->Ldomain() - eddySize);
540 else
541 g = 1.0/domn->Ldomain();
542
543 Pa = dtSample * invTauEddy * domn->pram->C_param / (f * g);
544
545}
546
548
551void eddy::applyVelocityKernels(domain *line, const int iS, const int iE) {
552
554
555 int npts = iE-iS+1;
556
557 set_kernel_coefficients();
558
559 if(domn->pram->Lspatial) {
560
561 bool Lflag = false; // if kernel results in negative (or < umin) uvel for Lspatial, then don't apply kernels.
562 for(int i=0; i<npts; i++) // a better approach would be to find the limiting alpha and use that.
563 if ( line->uvel->d.at(i+iS) + cCoef.at(0)*K.at(i) + bCoef.at(0)*abs(K.at(i)) < domn->pram->umin_spatial ) {
564 Lflag = true;
565 break;
566 }
567 if(!Lflag) {
568 for(int i=0; i<npts; i++) {
569 line->uvel->d.at(i+iS) += cCoef.at(0)*K.at(i) + bCoef.at(0)*abs(K.at(i));
570 line->vvel->d.at(i+iS) += cCoef.at(1)*K.at(i) + bCoef.at(1)*abs(K.at(i));
571 line->wvel->d.at(i+iS) += cCoef.at(2)*K.at(i) + bCoef.at(2)*abs(K.at(i));
572 }
573 }
574 }
575 else {
576 for(int i=0; i<npts; i++) {
577 line->uvel->d.at(i+iS) += cCoef.at(0)*K.at(i) + bCoef.at(0)*abs(K.at(i));
578 line->vvel->d.at(i+iS) += cCoef.at(1)*K.at(i) + bCoef.at(1)*abs(K.at(i));
579 line->wvel->d.at(i+iS) += cCoef.at(2)*K.at(i) + bCoef.at(2)*abs(K.at(i));
580 }
581 }
582
583}
584
586
590void eddy::fillKernel_planarAnalytic() {
591
592 K.resize(eddl->ngrd);
593 double L = eddySize;
594 double y0 = leftEdge;
595 double y;
596
597 for(int i=0; i<eddl->ngrd; i++) {
598
599 y = eddl->pos->d.at(i);
600
601 if(y >= y0 && y <= y0 + L/3)
602 K.at(i) = y- (y0 + 3*(y-y0));
603 else if(y >= y0 + L/3 && y <= y0 + 2*L/3)
604 K.at(i) = y- (y0 + 2*L - 3*(y-y0));
605 else if(y >= y0 + 2*L/3 && y <= y0 + L)
606 K.at(i) = y- (y0 + 3*(y-y0)-2*L);
607 else
608 K.at(i) = 0.0;
609 }
610
611}
612
614
618void eddy::fillKernel() {
619
620 int nseg = eddl->ngrd / 3;
621
622 K.resize(eddl->ngrd);
623
624 //---------- 1st Segment
625
626 K.at(0) = 0.5*(eddl->posf->d.at(1) + leftEdge) - pos0.at(0); // pos is 0.5*(faces), but we want pos=0.5*(face+eddy_edge)
627 for(int i=1, j=1; i<nseg; i++, j++) // pos0 was already updated this way in tripmap
628 K.at(i) = eddl->pos->d.at(i) - pos0.at(j);
629
630 //---------- Second Segment
631
632 for(int i=nseg*2-1, j=0; i>=nseg; i--, j++)
633 K.at(i) = eddl->pos->d.at(i) - pos0.at(j);
634
635 //---------- Third Segment
636
637 for(int i=nseg*2, j=0; i<nseg*3-1; i++, j++)
638 K.at(i) = eddl->pos->d.at(i) - pos0.at(j);
639 int i=nseg*3-1, j=nseg-1;
640 K.at(i) = 0.5*(rightEdge + eddl->posf->d.at(i)) - pos0.at(j); // pos is 0.5*(faces), but we want pos=0.5*(face+eddy_edge)
641
642}
643
645
652double eddy::eddyFavreAvgVelocity(const vector<double> &dxc) {
653
654 double ufavg = 0.0;
655 double ravg = 0.0;
656 double dmb;
657 for(int i=0; i<eddl->ngrd; i++) {
658 dmb = eddl->rho->d.at(i) * dxc.at(i);
659 ufavg += dmb*eddl->uvel->d.at(i);
660 ravg += dmb;
661 }
662 return ufavg/ravg;
663}
664
domain::ngrd
int ngrd
number of grid cells
Definition domain.h:42
domain::uvel
dv * uvel
Definition domain.h:51
domain::Ldomain
double Ldomain()
Definition domain.cc:156
domain::posf
dv * posf
access as: posf->d[i], or posf->var_name, etc.
Definition domain.h:48
domain::vvel
dv * vvel
Definition domain.h:52
domain::wvel
dv * wvel
Definition domain.h:53
domain::aDL
dv * aDL
Definition domain.h:62
domain::pos
dv * pos
pointers to gas properties
Definition domain.h:47
domain::mesher
meshManager * mesher
pointer to mesh manager object
Definition domain.h:78
domain::rand
randomGenerator * rand
Definition domain.h:80
domain::dvisc
dv * dvisc
Definition domain.h:50
domain::pram
param * pram
pointer to the parameters object
Definition domain.h:73
domain::rho
dv * rho
Definition domain.h:49
dv::d
vector< double > d
the data
Definition dv.h:30
eddy::LperiodicEddy
bool LperiodicEddy
a wrap-around eddy
Definition eddy.h:35
eddy::Pa
double Pa
eddy acceptance probability
Definition eddy.h:34
eddy::esdp1
double esdp1
eddy size distribution parameters.
Definition eddy.h:42
eddy::fillKernel_planarAnalytic
void fillKernel_planarAnalytic()
Definition eddy.cc:590
eddy::esdp4
double esdp4
Definition eddy.h:45
eddy::rightEdge
double rightEdge
right edge location of eddy
Definition eddy.h:32
eddy::ccb
vector< double > ccb
Definition eddy.h:47
eddy::cCoef
vector< double > cCoef
coefficient of K kernel
Definition eddy.h:36
eddy::init
void init(domain *p_domn, domain *p_eddl)
Definition eddy.cc:19
eddy::applyVelocityKernels
void applyVelocityKernels(domain *line, const int iS, const int iE)
Definition eddy.cc:551
eddy::computeEddyAcceptanceProb
void computeEddyAcceptanceProb(const double dtSample)
Definition eddy.cc:528
eddy::eddl
domain * eddl
pointer to eddy line object
Definition eddy.h:28
eddy::invTauEddy
double invTauEddy
inverse eddy timescale
Definition eddy.h:33
eddy::bCoef
vector< double > bCoef
coefficient of J kernel
Definition eddy.h:37
eddy::ccd
vector< double > ccd
polynomial coefficient arrays for cylindricalAnomalyHack
Definition eddy.h:47
eddy::ccc
vector< double > ccc
Definition eddy.h:47
eddy::eddySize
double eddySize
size of eddy
Definition eddy.h:30
eddy::dxc
vector< double > dxc
\delta(x^cCoord) is prop. to cell "volume"
Definition eddy.h:39
eddy::leftEdge
double leftEdge
left edge location of eddy
Definition eddy.h:31
eddy::set_kernel_coefficients
void set_kernel_coefficients()
Definition eddy.cc:430
eddy::pos0
vector< double > pos0
initial eddy cell locations, for kernel
Definition eddy.h:40
eddy::domn
domain * domn
pointer to domain object
Definition eddy.h:27
eddy::eddyTau
bool eddyTau(const double Z_value, const double C)
Definition eddy.cc:285
eddy::tripMap
void tripMap(domain *line, const int iS, int iE, const double C, const bool LsplitAtEddy=false)
Definition eddy.cc:116
eddy::sampleEddySize
void sampleEddySize()
Definition eddy.cc:71
eddy::K
vector< double > K
eddy kernel K
Definition eddy.h:38
eddy::esdp2
double esdp2
Definition eddy.h:43
eddy::sampleEddyPosition
void sampleEddyPosition()
Definition eddy.cc:86
eddy::esdp3
double esdp3
Definition eddy.h:44
eddy::fillKernel
void fillKernel()
Definition eddy.cc:618
eddy::eddyFavreAvgVelocity
double eddyFavreAvgVelocity(const vector< double > &dxc)
Definition eddy.cc:652
eddy::cca
vector< double > cca
Definition eddy.h:47
meshManager::splitCell
void splitCell(const int isplt, const int nsplt, const vector< double > &cellFaces)
Definition meshManager.cc:971
meshManager::setGridDxc
void setGridDxc(const domain *line, vector< double > &dxc, double C)
Definition meshManager.cc:1432
param::cCoord
int cCoord
1 = planar, 2 = cylindrical, 3 = spherical
Definition param.h:68
param::LPeEddy
bool LPeEddy
flag to turn on potential energy for eddies (vertical domain)
Definition param.h:57
param::LTMA
bool LTMA
true for the triplet map TMA: 3 = vol segments; false for TMB: 3 equal length segments
Definition param.h:63
param::g
double g
gravity (default -9.81)
Definition param.h:59
param::Lspatial
bool Lspatial
spatial formulation if true
Definition param.h:62
param::Lperiodic
bool Lperiodic
periodic if true
Definition param.h:61
param::Lp
double Lp
Most probable eddy size frac of domainLength.
Definition param.h:83
param::Lmax
double Lmax
Max eddy size frac of domainLength.
Definition param.h:84
param::C_param
double C_param
Eddy frequency parameter.
Definition param.h:46
param::A_param
double A_param
Energy Distribution parameter alpha.
Definition param.h:45
param::LdoDL
bool LdoDL
flag to do the DL energy from the DL instability
Definition param.h:54
param::LplanarTau
bool LplanarTau
true for computing cylindrical/spherical tau_eddy using a planar formulation. If accepted,...
Definition param.h:64
param::umin_spatial
double umin_spatial
min u for spatial flows; used when kernels pull velocity
Definition param.h:110
param::Lmin
double Lmin
Min eddy size frac of domainLength.
Definition param.h:85
domain.h
Header file for class domain.
eddy.h
Header file for class eddy.
Generated by doxygen 1.10.0