proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
SedClosure.h
Go to the documentation of this file.
1#ifndef SEDCLOSURE_H
2#define SEDCLOSURE_H
3
4#include <cmath>
5#include <iostream>
6#include "xtensor-python/pyarray.hpp"
7
8namespace proteus
9{
10 template<int nSpace>
11 class cppHsuSedStress
12{
13public:
14 cppHsuSedStress(
15 double aDarcy, // darcy parameter for drag term. Default value from Ergun (1952) is 150
16 double betaForch, // forchheimer parameter for drag term. Default value from Ergun (1952) is 1.75
17 double grain, // Grain size, default assumed as d50
18 double packFraction, //Critical volume fraction for switching the drag relation 0.2 by default, see Chen and Hsu 2014
19 double packMargin, // For packFraction +- packmargin, the drag coefficient is calculated by taking a weighted combination of the two relation (for packed and non-packed sediment
20 double maxFraction,
21 double frFraction,
22 double sigmaC,
23 double C3e,
24 double C4e,
25 double eR,
26 double fContact,
27 double mContact,
28 double nContact,
29 double angFriction,
30 double vos_limiter,
31 double mu_fr_limiter
32):
33
34 aDarcy_(aDarcy),
35 betaForch_(betaForch),
36 grain_(grain),
37 packFraction_(packFraction),
38 packMargin_(packMargin),
39 frFraction_(frFraction),
40 maxFraction_(maxFraction),
41 sigmaC_(sigmaC),
42 C3e_(C3e),
43 C4e_(C4e),
44 eR_(eR),
45 fContact_(fContact),
46 mContact_(mContact),
47 nContact_(nContact),
48 angFriction_(angFriction),
49 small_(1e-100),
50 notSoLarge_(1e+6),
51 large_(1e+100),
52 vos_limiter_(vos_limiter),
53 mu_fr_limiter_(mu_fr_limiter)
54
55
56
57 {}
58
59 inline double xt_betaCoeff(
60 double sedF, // Sediment fraction
61 double rhoFluid,
62 const xt::pyarray<double>& uFluid, //Fluid velocity
63 const xt::pyarray<double>& uSolid, //Sediment velocity
64 double nu //Kinematic viscosity
65 )
66 {
67 return betaCoeff(sedF, rhoFluid, uFluid.data(), uSolid.data(), nu);
68 }
69
70
71 inline double betaCoeff(
72 double sedF, // Sediment fraction
73 double rhoFluid,
74 const double uFluid[nSpace], //Fluid velocity
75 const double uSolid[nSpace], //Sediment velocity
76 double nu //Kinematic viscosity
77 )
78 {
79 double du2 = 0.;
80 for (int ii=0; ii<nSpace; ii++)
81 {
82 du2+= (uFluid[ii] - uSolid[ii])*(uFluid[ii] - uSolid[ii]);
83 }
84 double du = sqrt(du2);
85 double weight = 1.;
86
87
88 double gDrag1 = (aDarcy_*sedF*nu/((1. - sedF)*grain_*grain_) + betaForch_*du/grain_); // Darcy forchheimer term
89 double Cd = 0.44;
90 double Re = (1-sedF)*du*grain_/nu;
91 double Re_min = 1.0e-3;
92 if (Re < 1000)
93 {
94 Cd = 24*(1. + 0.15*pow(fmax(Re_min,Re), 0.687))/fmax(Re_min,Re); //Particle cloud resistance Wen and Yu 1966
95 }
96
97 double gDrag2 = 0.75 * Cd * du * pow(1. - sedF, -1.65)/grain_;
98
99 //cek debug -- this makes the drag term nonlinear
100 if(sedF < packFraction_ + packMargin_)
101 {
102 if (sedF > packFraction_ - packMargin_)
103 {
104 weight = 0.5 + 0.5* (sedF - packFraction_) /packMargin_;
105 }
106 else
107 {
108 weight = 0.;
109 }
110 }
111
112 return (weight*gDrag1 + (1.-weight)*gDrag2)*rhoFluid;
113 }
114
115 inline double gs0(double sedF)
116 {
117 double g0(0.);
118 if(sedF< 0.635)
119 {
120 if(sedF< 0.49)
121 {
122 g0 = 0.5*(2.-sedF)/((1-sedF)*(1-sedF)*(1-sedF));
123 }
124 else
125 {
126 g0= 0.853744035/(0.64-sedF);
127 }
128 }
129 else g0 = 170.74880702;
130 return g0;
131 }
132
133 inline double xt_deps_sed_deps(
134 double sedF, // Sediment fraction
135 double rhoFluid,
136 double rhoSolid,
137 const xt::pyarray<double>& uFluid,
138 const xt::pyarray<double>& uSolid,
139 const xt::pyarray<double>& gradC, //Sediment velocity
140 double nu, //Kinematic viscosity
141 double theta_n,
142 double kappa_n,
143 double epsilon_n,
144 double nuT_n,
145 const xt::pyarray<double>& g)
146 {
147 return deps_sed_deps(sedF,
148 rhoFluid,
149 rhoSolid,
150 uFluid.data(),
151 uSolid.data(),
152 gradC.data(),
153 nu,
154 theta_n,
155 kappa_n,
156 epsilon_n,
157 nuT_n,
158 g.data());
159
160 }
161
162 inline double deps_sed_deps(
163 double sedF, // Sediment fraction
164 double rhoFluid,
165 double rhoSolid,
166 const double uFluid[nSpace],
167 const double uSolid[nSpace],
168 const double gradC[nSpace], //Sediment velocity
169 double nu, //Kinematic viscosity
170 double theta_n,
171 double kappa_n,
172 double epsilon_n,
173 double nuT_n,
174 const double g[nSpace])
175
176 {
177
178 double beta = betaCoeff(sedF,rhoFluid,uFluid,uSolid,nu)+small_;
179 double gs = gs0(sedF)+small_;
180 double l_c = sqrt(M_PI)*grain_/(24.*(sedF+small_)*gs);
181 double t_p = rhoSolid/beta;
182 double t_c = l_c/(sqrt(theta_n) + small_);
183 double t_l = 0.165*kappa_n/(epsilon_n + small_);
184 double t_cl = std::min(t_c,t_l);
185 double alpha= t_cl/(t_cl + t_p);
186 double term = beta/(rhoFluid*(1.-sedF));
187 double es_1 = 2.*term*(1-alpha)*sedF*kappa_n;
188 double g_gradC = 0.;
189 for (int ii=0; ii<nSpace; ii++)
190 {
191 g_gradC+= g[ii]*gradC[ii];
192 }
193
194 double ss = rhoSolid/rhoFluid - 1.;
195 double es_2 = nuT_n*ss*g_gradC/sigmaC_/(1.-sedF) ;
196
197 return C3e_ * es_1 / kappa_n +C4e_ * es_2 / kappa_n;
198
199 }
200
201
202 inline double xt_kappa_sed1(
203 double sedF, // Sediment fraction
204 double rhoFluid,
205 double rhoSolid,
206 const xt::pyarray<double>& uFluid,
207 const xt::pyarray<double>& uSolid,
208 const xt::pyarray<double>& gradC,
209 double nu,
210 double theta_n,
211 double kappa_n,
212 double epsilon_n,
213 double nuT_n,
214 const xt::pyarray<double>& g)
215
216 {
217 return kappa_sed1(sedF,
218 rhoFluid,
219 rhoSolid,
220 uFluid.data(),
221 uSolid.data(),
222 gradC.data(),
223 nu,
224 theta_n,
225 kappa_n,
226 epsilon_n,
227 nuT_n,
228 g.data());
229 }
230
231 inline double kappa_sed1(
232 double sedF, // Sediment fraction
233 double rhoFluid,
234 double rhoSolid,
235 const double uFluid[nSpace],
236 const double uSolid[nSpace],
237 const double gradC[nSpace],
238 double nu,
239 double theta_n,
240 double kappa_n,
241 double epsilon_n,
242 double nuT_n,
243 const double g[nSpace])
244 {
245
246 double beta = betaCoeff(sedF,rhoFluid,uFluid,uSolid,nu)+small_;
247 double gs = gs0(sedF)+small_;
248 double l_c = sqrt(M_PI)*grain_/(24.*(sedF+small_)*gs);
249 double t_p = rhoSolid/beta;
250 double t_c = l_c/(sqrt(theta_n) + small_);
251 double t_l = 0.165*kappa_n/(epsilon_n + small_);
252 double t_cl = std::min(t_c,t_l);
253 double alpha= t_cl/(t_cl + t_p);
254 double term = beta/(rhoFluid*(1.-sedF));
255 double es_1 = 2.*term*(1-alpha)*sedF*kappa_n;
256 double g_gradC = 0.;
257 for (int ii=0; ii<nSpace; ii++)
258 {
259 g_gradC+= g[ii]*gradC[ii];
260 }
261
262 double ss = rhoSolid/rhoFluid - 1.;
263 double es_2 = nuT_n*ss*g_gradC/sigmaC_/(1.-sedF) ;
264
265 return es_1 + es_2;
266
267 }
268
269 inline double xt_dkappa_sed1_dk(
270 double sedF, // Sediment fraction
271 double rhoFluid,
272 double rhoSolid,
273 const xt::pyarray<double>& uFluid, //Fluid velocity
274 const xt::pyarray<double>& uSolid, //Sediment velocity
275 const xt::pyarray<double>& gradC, //Sediment velocity
276 double nu, //Kinematic viscosity
277 double theta_n,
278 double kappa_n,
279 double epsilon_n,
280 double nuT_n)
281 {
282 return dkappa_sed1_dk(sedF,
283 rhoFluid,
284 rhoSolid,
285 uFluid.data(),
286 uSolid.data(),
287 gradC.data(),
288 nu,
289 theta_n,
290 kappa_n,
291 epsilon_n,
292 nuT_n);
293
294 }
295
296 inline double dkappa_sed1_dk(
297 double sedF, // Sediment fraction
298 double rhoFluid,
299 double rhoSolid,
300 const double uFluid[nSpace], //Fluid velocity
301 const double uSolid[nSpace], //Sediment velocity
302 const double gradC[nSpace], //Sediment velocity
303 double nu, //Kinematic viscosity
304 double theta_n,
305 double kappa_n,
306 double epsilon_n,
307 double nuT_n)
308 {
309 double beta = betaCoeff(sedF,rhoFluid,uFluid,uSolid,nu)+small_;
310 double gs = gs0(sedF)+small_;
311 double l_c = sqrt(M_PI)*grain_/(24.*(sedF+small_)*gs);
312 double t_p = rhoSolid/beta;
313 double t_c = l_c/(sqrt(theta_n) + small_);
314 double t_l = 0.165*kappa_n/(epsilon_n + small_);
315 double t_cl = std::min(t_c,t_l);
316 double alpha= t_cl/(t_cl + t_p);
317 double term = beta/(rhoFluid*(1.-sedF));
318 double es_1 = 2.*term*rhoSolid*(1-alpha)*sedF;
319 return es_1;
320
321 }
322
323
324
325 inline double psc( double sedF,
326 double rhoSolid,
327 double theta )
328
329 {
330 double gs_0 = gs0(sedF);
331 double eRp1 = 1.+eR_;
332 double psc = rhoSolid * sedF * (1.+2.*eRp1*sedF*gs_0)*theta;
333 return psc;
334 }
335
336 inline double psc_term( double sedF,
337 double rhoSolid,
338 double theta_np1,
339 double du_dx,
340 double dv_dy,
341 double dw_dz)
342
343 {
344
345 return -2.*psc(sedF,rhoSolid,theta_np1)*(du_dx + dv_dy + dw_dz)/(3.*rhoSolid * sedF);
346 }
347
348 inline double dpsc_term_dtheta( double sedF,
349 double rhoSolid,
350 double du_dx,
351 double dv_dy,
352 double dw_dz)
353
354 {
355
356 double gs_0 = gs0(sedF);
357 double eRp1 = 1.+eR_;
358 double dpsc = rhoSolid * sedF * (1.+2.*eRp1*sedF*gs_0);
359 return -2.*dpsc*(du_dx + dv_dy + dw_dz)/(3.*rhoSolid * sedF);
360 }
361
362
363 inline double mu_sc( double sedF,
364 double rhoSolid,
365 double theta )
366
367 {
368 double gs_0 = gs0(sedF);
369 double sq_pi = (sqrt(M_PI));
370 double sedF2 = sedF * sedF;
371 double eRp1 = 1.+eR_;
372 double sq_theta = sqrt(theta);
373
374 double mu_sc = rhoSolid*grain_*sq_theta*( 0.8*sedF2*gs_0*eRp1/(sq_pi) + (1./15.)*sq_pi*sedF2*gs_0*eRp1 + (1./6.)*sq_pi*sedF + (5./48)*sq_pi/(gs_0*eRp1) );
375
376 return mu_sc;
377 }
378
379 inline double l_sc( double sedF,
380 double rhoSolid,
381 double theta_n )
382
383 {
384 double gs_0 = gs0(sedF);
385 double sq_pi = (sqrt(M_PI));
386 double sedF2 = sedF * sedF;
387 double eRp1 = 1.+eR_;
388 double sq_theta = sqrt(theta_n);
389 double lambda = (4./3.)*sedF2*rhoSolid*grain_*gs_0*eRp1*sq_theta/sq_pi;
390
391 return lambda;
392 }
393
394
395
396
397
398
399
400
401 inline double tausc_term_theta(
402 double sedF, // IN
403 double rhoSolid, //IN
404 double theta_n, //IN
405 double du_dx,
406 double du_dy,
407 double du_dz,
408 double dv_dx,
409 double dv_dy,
410 double dv_dz,
411 double dw_dx,
412 double dw_dy,
413 double dw_dz)
414
415
416 {
417
418 double la = l_sc(sedF, rhoSolid, theta_n);
419 double mu = mu_sc(sedF, rhoSolid, theta_n);
420
421
422 double divU = du_dx + dv_dy + dw_dz;
423 double t11 = 2*mu*du_dx + (la - (2./3.)*mu)*divU ;
424 double t22 = 2*mu*dv_dy + (la - (2./3.)*mu)*divU ;
425 double t33 = 2*mu*dw_dz + (la - (2./3.)*mu)*divU ;
426 double t13 = mu*(du_dz + dw_dx);
427 double t23 = mu*(dv_dz + dw_dy);
428 double t12 = mu*(du_dy + dv_dx);
429 // No need to define t31, t32 and t21 as tensor is symmetric
430
431 double term = t11 * du_dx + t22 * dv_dy + t33 * dw_dz;
432 term += t12*(du_dy +dv_dx);
433 term += t13*(du_dz +dw_dx);
434 term += t23*(dv_dz +dw_dy);
435
436 return term* (2./3.) / (rhoSolid*sedF);
437
438 }
439
440 inline double gamma_s( double sedF,
441 double rhoSolid,
442 double theta_n,
443 double theta_np1,
444 double du_dx,
445 double dv_dy,
446 double dw_dz)
447
448
449 {
450 double sq_pi = (sqrt(M_PI));
451 double sq_theta = sqrt(theta_n);
452 double divU = du_dx + dv_dy + dw_dz;
453 double gamma_s = - 2*(1-eR_*eR_)*sedF*gs0(sedF)*(4.*sq_theta/(sq_pi*grain_) - divU)*theta_np1;
454 return gamma_s;
455
456 }
457 inline double dgamma_s_dtheta( double sedF,
458 double rhoSolid,
459 double theta_n,
460 double du_dx,
461 double dv_dy,
462 double dw_dz)
463
464
465 {
466 double sq_pi = (sqrt(M_PI));
467 double sq_theta = sqrt(theta_n);
468 double divU = du_dx + dv_dy + dw_dz;
469 double gamma_s = - 2*(1-eR_*eR_)*sedF*gs0(sedF)*(4.*sq_theta/(sq_pi*grain_) - divU);
470 return gamma_s;
471
472 }
473 inline double xt_jint1( double sedF,
474 double rhoFluid,
475 double rhoSolid,
476 const xt::pyarray<double>& uFluid,
477 const xt::pyarray<double>& uSolid,
478 double kappa,
479 double epsilon,
480 double theta,
481 double nu)
482 {
483 return jint1(sedF,
484 rhoFluid,
485 rhoSolid,
486 uFluid.data(),
487 uSolid.data(),
488 kappa,
489 epsilon,
490 theta,
491 nu);
492 }
493
494 inline double jint1( double sedF,
495 double rhoFluid,
496 double rhoSolid,
497 const double uFluid[nSpace],
498 const double uSolid[nSpace],
499 double kappa,
500 double epsilon,
501 double theta,
502 double nu)
503 {
504
505 double beta = betaCoeff(sedF,rhoFluid,uFluid,uSolid,nu)+small_;
506 double gs = gs0(sedF)+small_;
507 double l_c = sqrt(M_PI)*grain_/(24.*(sedF+small_)*gs);
508 double t_p = rhoSolid/beta;
509 double t_c = l_c/(sqrt(theta) + small_);
510 double t_l = 0.165*kappa/(epsilon + small_);
511 double t_cl = std::min(t_c,t_l);
512 double alpha= t_cl/(t_cl + t_p);
513
514 double Jint1 = 4. * alpha * beta * kappa /( 3.* rhoSolid) ;
515 return Jint1;
516
517 }
518 inline double xt_jint2( double sedF,
519 double rhoFluid,
520 double rhoSolid,
521 const xt::pyarray<double>& uFluid,
522 const xt::pyarray<double>& uSolid,
523 double theta,
524 double nu)
525
526 {
527 return jint2(sedF,
528 rhoFluid,
529 rhoSolid,
530 uFluid.data(),
531 uSolid.data(),
532 theta,
533 nu);
534 }
535
536 inline double jint2( double sedF,
537 double rhoFluid,
538 double rhoSolid,
539 const double uFluid[nSpace],
540 const double uSolid[nSpace],
541 double theta,
542 double nu)
543 {
544
545 double beta = betaCoeff(sedF,rhoFluid,uFluid,uSolid,nu)+small_;
546 return - 2*beta*theta/rhoSolid;
547
548 }
549
550 inline double xt_djint2_dtheta( double sedF,
551 double rhoFluid,
552 double rhoSolid,
553 const xt::pyarray<double>& uFluid,
554 const xt::pyarray<double>& uSolid,
555 double nu)
556 {
557 return djint2_dtheta(sedF,
558 rhoFluid,
559 rhoSolid,
560 uFluid.data(),
561 uSolid.data(),
562 nu);
563 }
564
565 inline double djint2_dtheta( double sedF,
566 double rhoFluid,
567 double rhoSolid,
568 const double uFluid[nSpace],
569 const double uSolid[nSpace],
570 double nu)
571 {
572
573 double beta = betaCoeff(sedF,rhoFluid,uFluid,uSolid,nu)+small_;
574 return - 2*beta/rhoSolid;
575
576 }
577
578 inline double k_diff( double sedF,
579 double rhoSolid,
580 double theta )
581
582 {
583 double gs_0 = gs0(sedF);
584 double sq_pi = (sqrt(M_PI));
585 double sedF2 = sedF * sedF;
586 double eRp1 = 1.+eR_;
587 double sq_theta = sqrt(theta);
588 double k_diff = rhoSolid*grain_*sq_theta*( 2.*sedF2*gs_0*eRp1/(sq_pi) + (0.5625)*sq_pi*sedF2*gs_0*eRp1 + (0.9375)*sq_pi*sedF + (0.390625)*sq_pi/(gs_0*eRp1) );
589
590
591 return k_diff;
592 }
593
594 inline double p_friction(double sedF)
595
596 {
597 double pf = 0.;
598 if (sedF > frFraction_)
599 {
600 double sedLim =std::min(sedF,vos_limiter_);
601 pf =fContact_*pow(sedLim-frFraction_,mContact_) / ( pow(maxFraction_ - sedLim,nContact_) + small_);
602 }
603
604 return pf;
605
606 }
607
608
609 inline double gradp_friction(double sedF)
610
611 {
612 double coeff = 0;
613 double pf = p_friction(sedF);
614 if (sedF > frFraction_ )
615 {
616 double sedLim =std::min(sedF,vos_limiter_);
617
618 double den1 = (sedLim - frFraction_);
619 double den2 = (maxFraction_ - sedLim);
620
621 coeff = pf *( (mContact_/den1) + (nContact_/den2) );
622 }
623
624 return coeff;
625 }
626 inline double mu_fr(double sedF,
627 double du_dx,
628 double du_dy,
629 double du_dz,
630 double dv_dx,
631 double dv_dy,
632 double dv_dz,
633 double dw_dx,
634 double dw_dy,
635 double dw_dz)
636 {
637 double divU = du_dx + dv_dy + dw_dz;
638 double pf = p_friction(sedF);
639 double s11 = du_dx - (1./nSpace)*divU;
640 double s22 = dv_dy - (1./nSpace)*divU;
641 double s33 = dw_dz - (1./nSpace)*divU;
642 double s12 = 0.5*(du_dy + dv_dx);
643 double s13 = 0.5*(du_dz + dw_dx);
644 double s23 = 0.5*(dv_dz + dw_dy);
645 double sumS = s11*s11 + s22*s22 + s33*s33 + 2.*s12*s12 + 2.*s13*s13 + 2.*s23*s23;
646 double mu_sf = 0.0;
647 double mu_fr = 0.0;
648 if (sedF > frFraction_)
649 {
650 mu_sf = pf * sqrt(2.) * sin(angFriction_) / (2 * sqrt(sumS) + small_);
651 mu_fr = std::min(mu_sf,mu_fr_limiter_);
652 }
653// if (sedF > 0.0 )
654// {
655// printf("sedF --> %2.20f", sedF);
656// }
657// if (mu_sf > 0.0 )
658// {
659// printf("mu_sf --> %2.20f", mu_sf);
660// }
661// if (mu_fr > 0.0 )
662// {
663// printf("mu_fr --> %2.20f", mu_fr);
664// }
665// if (pf > 0.0 )
666// {
667// printf("pf --> %2.20f", pf);
668// }
669
670 return mu_fr;
671 }
672
673
674
675
676
677
678 /*
679 inline double diffusion_theta_rhs( double sedF,
680 double rhoSolid,
681 double theta_n,
682 double dtheta_dx,
683 double dtheta_dy,
684 double dtheta_dz
685 )
686 {
687
688 double gs_0 = gs0(sedF);
689 double sq_pi = (sqrt(M_PI));
690 double sedF2 = sedF * sedF;
691 double eRp1 = 1.+eR_;
692 double sq_theta = sqrt(theta_n);
693
694 double divTheta = dtheta_dx + dtheta_dy + dtheta_dz;
695
696 double k_diff = k_diff(rhoSolid);
697
698 return kappa*divTheta;
699
700 }
701
702
703
704
705 }*/
706
707
708
709
710
711 inline xt::pyarray<double> xt_mIntFluid(double sedF,
712 double rhoFluid,
713 const xt::pyarray<double>& uFluid_n, //Fluid velocity
714 const xt::pyarray<double>& uSolid_n, //Sediment velocity
715 const xt::pyarray<double>& uFluid_np1, //Fluid velocity
716 double nu, //Kinematic viscosity
717 double nuT, //Turbulent viscosity
718 const xt::pyarray<double>& gradc
719 )
720 {
721 auto mint2 = xt::pyarray<double>::from_shape({2});
722 mIntFluid(mint2.data(),
723 sedF,
724 rhoFluid,
725 uFluid_n.data(),
726 uSolid_n.data(),
727 uFluid_np1.data(),
728 nu,
729 nuT,
730 gradc.data());
731 return mint2;
732 }
733
734 inline void mIntFluid( double * mint2, double sedF,
735 double rhoFluid,
736 const double uFluid_n[nSpace], //Fluid velocity
737 const double uSolid_n[nSpace], //Sediment velocity
738 const double uFluid_np1[nSpace], //Fluid velocity
739 double nu, //Kinematic viscosity
740 double nuT, //Turbulent viscosity
741 const double gradc[nSpace]
742 )
743 {
744
745 double beta = betaCoeff(sedF,rhoFluid,uFluid_n,uSolid_n,nu);
746 for (int ii=0; ii<nSpace; ii++)
747 {
748 mint2[ii] = -sedF*beta*(uFluid_np1[ii])/(rhoFluid * (1. - sedF)) ;
749 }
750
751 }
752
753 inline xt::pyarray<double> xt_mIntSolid(double sedF,
754 double rhoFluid,
755 const xt::pyarray<double>& uFluid_n, //Fluid velocity
756 const xt::pyarray<double>& uSolid_n, //Sediment velocity
757 const xt::pyarray<double>& uSolid_np1, //Sediment velocity
758 double nu, //Kinematic viscosity
759 double nuT, //Turbulent viscosity
760 const xt::pyarray<double>& gradc
761 )
762 {
763 auto mint2 = xt::pyarray<double>::from_shape({2});
764 mIntSolid(mint2.data(),
765 sedF,
766 rhoFluid,
767 uFluid_n.data(),
768 uSolid_n.data(),
769 uSolid_np1.data(),
770 nu,
771 nuT,
772 gradc.data());
773 return mint2;
774 }
775
776 inline void mIntSolid( double * mint2, double sedF,
777 double rhoFluid,
778 const double uFluid_n[nSpace], //Fluid velocity
779 const double uSolid_n[nSpace], //Sediment velocity
780 const double uSolid_np1[nSpace], //Sediment velocity
781 double nu, //Kinematic viscosity
782 double nuT, //Turbulent viscosity
783 const double gradc[nSpace]
784 )
785 {
786
787 double beta = betaCoeff(sedF,rhoFluid,uFluid_n,uSolid_n,nu);
788 for (int ii=0; ii<nSpace; ii++)
789 {
790 mint2[ii] = -sedF*beta*(-uSolid_np1[ii])/(rhoFluid * (1. - sedF)) ;
791 }
792
793 }
794 inline xt::pyarray<double> xt_mIntgradC(double sedF,
795 double rhoFluid,
796 const xt::pyarray<double>& uFluid_n, //Fluid velocity
797 const xt::pyarray<double>& uSolid_n, //Sediment velocity
798 double nu, //Kinematic viscosity
799 double nuT, //Turbulent viscosity
800 const xt::pyarray<double>& gradc
801 )
802 {
803
804 auto mint2 = xt::pyarray<double>::from_shape({2});
805 mIntgradC(mint2.data(),
806 sedF,
807 rhoFluid,
808 uFluid_n.data(),
809 uSolid_n.data(),
810 nu,
811 nuT,
812 gradc.data());
813 return mint2;
814 }
815
816 inline void mIntgradC(double * mint2, double sedF,
817 double rhoFluid,
818 const double uFluid_n[nSpace], //Fluid velocity
819 const double uSolid_n[nSpace], //Sediment velocity
820 double nu, //Kinematic viscosity
821 double nuT, //Turbulent viscosity
822 const double gradc[nSpace]
823 )
824
825 {
826
827 double beta = betaCoeff(sedF,rhoFluid,uFluid_n,uSolid_n,nu);
828 for (int ii=0; ii<nSpace; ii++)
829 {
830 mint2[ii] = - sedF*beta*nuT*gradc[ii]/sigmaC_/(rhoFluid * (1. - sedF));
831 }
832
833
834 }
835
836
837
838 inline double xt_dmInt_duFluid
839 ( double sedF,
840 double rhoFluid,
841 const xt::pyarray<double>& uFluid_n, //Fluid velocity
842 const xt::pyarray<double>& uSolid_n, //Sediment velocity
843 double nu //Kinematic viscosity
844
845 )
846 {
847 return dmInt_duFluid(sedF,
848 rhoFluid,
849 uFluid_n.data(),
850 uSolid_n.data(),
851 nu);
852 }
853
854 inline double dmInt_duFluid
855 ( double sedF,
856 double rhoFluid,
857 const double uFluid_n[nSpace], //Fluid velocity
858 const double uSolid_n[nSpace], //Sediment velocity
859 double nu //Kinematic viscosity
860
861 )
862 {
863 return -sedF*betaCoeff(sedF,rhoFluid,uFluid_n,uSolid_n,nu)/(rhoFluid * (1. - sedF));
864
865 }
866
867
868 inline double xt_dmInt_duSolid
869 ( double sedF,
870 double rhoFluid,
871 const xt::pyarray<double>& uFluid_n, //Fluid velocity
872 const xt::pyarray<double>& uSolid_n, //Sediment velocity
873 double nu //Kinematic viscosity
874
875 )
876 {
877 return dmInt_duSolid(sedF,
878 rhoFluid,
879 uFluid_n.data(),
880 uSolid_n.data(),
881 nu);
882 }
883
884 inline double dmInt_duSolid
885 ( double sedF,
886 double rhoFluid,
887 const double uFluid_n[nSpace], //Fluid velocity
888 const double uSolid_n[nSpace], //Sediment velocity
889 double nu //Kinematic viscosity
890
891 )
892 {
893
894 return +sedF*betaCoeff(sedF,rhoFluid,uFluid_n,uSolid_n,nu)/(rhoFluid * (1. - sedF));
895 }
896
897 inline double p_s( double sedF,
898 double rhoSolid,
899 double theta,
900 double du_dx,
901 double du_dy,
902 double du_dz,
903 double dv_dx,
904 double dv_dy,
905 double dv_dz,
906 double dw_dx,
907 double dw_dy,
908 double dw_dz)
909
910
911 {
912 double divU = du_dx + dv_dy + dw_dz;
913 double mf = mu_fr(sedF, du_dx, du_dy, du_dz, dv_dx, dv_dy, dv_dz, dw_dx, dw_dy, dw_dz);
914 double msc = mu_sc(sedF, rhoSolid, theta );
915 double lam = l_sc(sedF, rhoSolid, theta );
916 double pcorr = ( (2./3.)*(msc + mf) - lam ) * divU;
917 return p_friction(sedF) + psc(sedF, rhoSolid, theta) + pcorr ;
918 }
919
920
921
922
923
924
925 double aDarcy_;
926 double betaForch_;
927 double grain_;
928 double packFraction_;
929 double packMargin_;
930 double frFraction_;
931 double maxFraction_;
932 double sigmaC_;
933 double C3e_;
934 double C4e_;
935 double eR_;
936 double fContact_;
937 double mContact_;
938 double nContact_;
939 double angFriction_;
940 double small_;
941 double notSoLarge_;
942 double large_;
943 double vos_limiter_;
944 double mu_fr_limiter_;
945
946};
947
948 typedef cppHsuSedStress<2> cppHsuSedStress2D;
949}
950
951#endif
952
coeff
Double * coeff
Definition Headers.h:42
proteus::cppHsuSedStress::p_s
double p_s(double sedF, double rhoSolid, double theta, double du_dx, double du_dy, double du_dz, double dv_dx, double dv_dy, double dv_dz, double dw_dx, double dw_dy, double dw_dz)
Definition SedClosure.h:897
proteus::cppHsuSedStress::psc
double psc(double sedF, double rhoSolid, double theta)
Definition SedClosure.h:325
proteus::cppHsuSedStress::psc_term
double psc_term(double sedF, double rhoSolid, double theta_np1, double du_dx, double dv_dy, double dw_dz)
Definition SedClosure.h:336
proteus::cppHsuSedStress::xt_mIntSolid
xt::pyarray< double > xt_mIntSolid(double sedF, double rhoFluid, const xt::pyarray< double > &uFluid_n, const xt::pyarray< double > &uSolid_n, const xt::pyarray< double > &uSolid_np1, double nu, double nuT, const xt::pyarray< double > &gradc)
Definition SedClosure.h:753
proteus::cppHsuSedStress::xt_mIntgradC
xt::pyarray< double > xt_mIntgradC(double sedF, double rhoFluid, const xt::pyarray< double > &uFluid_n, const xt::pyarray< double > &uSolid_n, double nu, double nuT, const xt::pyarray< double > &gradc)
Definition SedClosure.h:794
proteus::cppHsuSedStress::tausc_term_theta
double tausc_term_theta(double sedF, double rhoSolid, double theta_n, double du_dx, double du_dy, double du_dz, double dv_dx, double dv_dy, double dv_dz, double dw_dx, double dw_dy, double dw_dz)
Definition SedClosure.h:401
proteus::cppHsuSedStress::k_diff
double k_diff(double sedF, double rhoSolid, double theta)
Definition SedClosure.h:578
proteus::cppHsuSedStress::xt_djint2_dtheta
double xt_djint2_dtheta(double sedF, double rhoFluid, double rhoSolid, const xt::pyarray< double > &uFluid, const xt::pyarray< double > &uSolid, double nu)
Definition SedClosure.h:550
proteus::cppHsuSedStress< 2 >::dmInt_duFluid
double dmInt_duFluid(double sedF, double rhoFluid, const double uFluid_n[nSpace], const double uSolid_n[nSpace], double nu)
Definition SedClosure.h:855
proteus::cppHsuSedStress< 2 >::jint1
double jint1(double sedF, double rhoFluid, double rhoSolid, const double uFluid[nSpace], const double uSolid[nSpace], double kappa, double epsilon, double theta, double nu)
Definition SedClosure.h:494
proteus::cppHsuSedStress< 2 >::deps_sed_deps
double deps_sed_deps(double sedF, double rhoFluid, double rhoSolid, const double uFluid[nSpace], const double uSolid[nSpace], const double gradC[nSpace], double nu, double theta_n, double kappa_n, double epsilon_n, double nuT_n, const double g[nSpace])
Definition SedClosure.h:162
proteus::cppHsuSedStress::xt_dkappa_sed1_dk
double xt_dkappa_sed1_dk(double sedF, double rhoFluid, double rhoSolid, const xt::pyarray< double > &uFluid, const xt::pyarray< double > &uSolid, const xt::pyarray< double > &gradC, double nu, double theta_n, double kappa_n, double epsilon_n, double nuT_n)
Definition SedClosure.h:269
proteus::cppHsuSedStress< 2 >::betaCoeff
double betaCoeff(double sedF, double rhoFluid, const double uFluid[nSpace], const double uSolid[nSpace], double nu)
Definition SedClosure.h:71
proteus::cppHsuSedStress< 2 >::mIntFluid
void mIntFluid(double *mint2, double sedF, double rhoFluid, const double uFluid_n[nSpace], const double uSolid_n[nSpace], const double uFluid_np1[nSpace], double nu, double nuT, const double gradc[nSpace])
Definition SedClosure.h:734
proteus::cppHsuSedStress::xt_betaCoeff
double xt_betaCoeff(double sedF, double rhoFluid, const xt::pyarray< double > &uFluid, const xt::pyarray< double > &uSolid, double nu)
Definition SedClosure.h:59
proteus::cppHsuSedStress::dpsc_term_dtheta
double dpsc_term_dtheta(double sedF, double rhoSolid, double du_dx, double dv_dy, double dw_dz)
Definition SedClosure.h:348
proteus::cppHsuSedStress::xt_deps_sed_deps
double xt_deps_sed_deps(double sedF, double rhoFluid, double rhoSolid, const xt::pyarray< double > &uFluid, const xt::pyarray< double > &uSolid, const xt::pyarray< double > &gradC, double nu, double theta_n, double kappa_n, double epsilon_n, double nuT_n, const xt::pyarray< double > &g)
Definition SedClosure.h:133
proteus::cppHsuSedStress< 2 >::kappa_sed1
double kappa_sed1(double sedF, double rhoFluid, double rhoSolid, const double uFluid[nSpace], const double uSolid[nSpace], const double gradC[nSpace], double nu, double theta_n, double kappa_n, double epsilon_n, double nuT_n, const double g[nSpace])
Definition SedClosure.h:231
proteus::cppHsuSedStress< 2 >::dkappa_sed1_dk
double dkappa_sed1_dk(double sedF, double rhoFluid, double rhoSolid, const double uFluid[nSpace], const double uSolid[nSpace], const double gradC[nSpace], double nu, double theta_n, double kappa_n, double epsilon_n, double nuT_n)
Definition SedClosure.h:296
proteus::cppHsuSedStress::xt_jint1
double xt_jint1(double sedF, double rhoFluid, double rhoSolid, const xt::pyarray< double > &uFluid, const xt::pyarray< double > &uSolid, double kappa, double epsilon, double theta, double nu)
Definition SedClosure.h:473
proteus::cppHsuSedStress< 2 >::mIntSolid
void mIntSolid(double *mint2, double sedF, double rhoFluid, const double uFluid_n[nSpace], const double uSolid_n[nSpace], const double uSolid_np1[nSpace], double nu, double nuT, const double gradc[nSpace])
Definition SedClosure.h:776
proteus::cppHsuSedStress::xt_dmInt_duSolid
double xt_dmInt_duSolid(double sedF, double rhoFluid, const xt::pyarray< double > &uFluid_n, const xt::pyarray< double > &uSolid_n, double nu)
Definition SedClosure.h:869
proteus::cppHsuSedStress::xt_jint2
double xt_jint2(double sedF, double rhoFluid, double rhoSolid, const xt::pyarray< double > &uFluid, const xt::pyarray< double > &uSolid, double theta, double nu)
Definition SedClosure.h:518
proteus::cppHsuSedStress::p_friction
double p_friction(double sedF)
Definition SedClosure.h:594
proteus::cppHsuSedStress< 2 >::mIntgradC
void mIntgradC(double *mint2, double sedF, double rhoFluid, const double uFluid_n[nSpace], const double uSolid_n[nSpace], double nu, double nuT, const double gradc[nSpace])
Definition SedClosure.h:816
proteus::cppHsuSedStress< 2 >::dmInt_duSolid
double dmInt_duSolid(double sedF, double rhoFluid, const double uFluid_n[nSpace], const double uSolid_n[nSpace], double nu)
Definition SedClosure.h:885
proteus::cppHsuSedStress< 2 >::jint2
double jint2(double sedF, double rhoFluid, double rhoSolid, const double uFluid[nSpace], const double uSolid[nSpace], double theta, double nu)
Definition SedClosure.h:536
proteus::cppHsuSedStress::gradp_friction
double gradp_friction(double sedF)
Definition SedClosure.h:609
proteus::cppHsuSedStress::l_sc
double l_sc(double sedF, double rhoSolid, double theta_n)
Definition SedClosure.h:379
proteus::cppHsuSedStress::mu_sc
double mu_sc(double sedF, double rhoSolid, double theta)
Definition SedClosure.h:363
proteus::cppHsuSedStress::xt_mIntFluid
xt::pyarray< double > xt_mIntFluid(double sedF, double rhoFluid, const xt::pyarray< double > &uFluid_n, const xt::pyarray< double > &uSolid_n, const xt::pyarray< double > &uFluid_np1, double nu, double nuT, const xt::pyarray< double > &gradc)
Definition SedClosure.h:711
proteus::cppHsuSedStress< 2 >::djint2_dtheta
double djint2_dtheta(double sedF, double rhoFluid, double rhoSolid, const double uFluid[nSpace], const double uSolid[nSpace], double nu)
Definition SedClosure.h:565
proteus::cppHsuSedStress::dgamma_s_dtheta
double dgamma_s_dtheta(double sedF, double rhoSolid, double theta_n, double du_dx, double dv_dy, double dw_dz)
Definition SedClosure.h:457
proteus::cppHsuSedStress::gs0
double gs0(double sedF)
Definition SedClosure.h:115
proteus::cppHsuSedStress::xt_dmInt_duFluid
double xt_dmInt_duFluid(double sedF, double rhoFluid, const xt::pyarray< double > &uFluid_n, const xt::pyarray< double > &uSolid_n, double nu)
Definition SedClosure.h:839
proteus::cppHsuSedStress::mu_fr
double mu_fr(double sedF, double du_dx, double du_dy, double du_dz, double dv_dx, double dv_dy, double dv_dz, double dw_dx, double dw_dy, double dw_dz)
Definition SedClosure.h:626
proteus::cppHsuSedStress::cppHsuSedStress
cppHsuSedStress(double aDarcy, double betaForch, double grain, double packFraction, double packMargin, double maxFraction, double frFraction, double sigmaC, double C3e, double C4e, double eR, double fContact, double mContact, double nContact, double angFriction, double vos_limiter, double mu_fr_limiter)
Definition SedClosure.h:14
proteus::cppHsuSedStress::xt_kappa_sed1
double xt_kappa_sed1(double sedF, double rhoFluid, double rhoSolid, const xt::pyarray< double > &uFluid, const xt::pyarray< double > &uSolid, const xt::pyarray< double > &gradC, double nu, double theta_n, double kappa_n, double epsilon_n, double nuT_n, const xt::pyarray< double > &g)
Definition SedClosure.h:202
proteus::cppHsuSedStress::gamma_s
double gamma_s(double sedF, double rhoSolid, double theta_n, double theta_np1, double du_dx, double dv_dy, double dw_dz)
Definition SedClosure.h:440
pf
#define pf(x)
Definition jf.h:23
proteus
Definition ADR.h:19
proteus::cppHsuSedStress2D
cppHsuSedStress< 2 > cppHsuSedStress2D
Definition SedClosure.h:948
Generated on for proteus by doxygen 1.18.0