proteus  1.8.1
C/C++/Fortran libraries
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 
8 namespace proteus
9 {
10  template<int nSpace>
11  class cppHsuSedStress
12 {
13 public:
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 
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::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::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::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::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::grain_
double grain_
Definition: SedClosure.h:927
proteus::cppHsuSedStress::mu_sc
double mu_sc(double sedF, double rhoSolid, double theta)
Definition: SedClosure.h:363
proteus::cppHsuSedStress::mu_fr_limiter_
double mu_fr_limiter_
Definition: SedClosure.h:944
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::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::cppHsuSedStress2D
cppHsuSedStress< 2 > cppHsuSedStress2D
Definition: SedClosure.h:948
proteus::cppHsuSedStress::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::fContact_
double fContact_
Definition: SedClosure.h:936
pf
#define pf(x)
Definition: jf.h:23
proteus::cppHsuSedStress::l_sc
double l_sc(double sedF, double rhoSolid, double theta_n)
Definition: SedClosure.h:379
proteus::cppHsuSedStress::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
coeff
Double * coeff
Definition: Headers.h:42
proteus::cppHsuSedStress::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::psc
double psc(double sedF, double rhoSolid, double theta)
Definition: SedClosure.h:325
proteus::cppHsuSedStress::C3e_
double C3e_
Definition: SedClosure.h:933
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::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::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
proteus::cppHsuSedStress::small_
double small_
Definition: SedClosure.h:940
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::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_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::frFraction_
double frFraction_
Definition: SedClosure.h:930
proteus::cppHsuSedStress::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::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::p_friction
double p_friction(double sedF)
Definition: SedClosure.h:594
min
#define min(a, b)
Definition: jf.h:71
proteus::cppHsuSedStress::C4e_
double C4e_
Definition: SedClosure.h:934
proteus::cppHsuSedStress::maxFraction_
double maxFraction_
Definition: SedClosure.h:931
proteus::cppHsuSedStress::betaCoeff
double betaCoeff(double sedF, double rhoFluid, const double uFluid[nSpace], const double uSolid[nSpace], double nu)
Definition: SedClosure.h:71
proteus::cppHsuSedStress::nContact_
double nContact_
Definition: SedClosure.h:938
proteus::cppHsuSedStress::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::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_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::angFriction_
double angFriction_
Definition: SedClosure.h:939
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::gs0
double gs0(double sedF)
Definition: SedClosure.h:115
proteus::cppHsuSedStress::gradp_friction
double gradp_friction(double sedF)
Definition: SedClosure.h:609
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::aDarcy_
double aDarcy_
Definition: SedClosure.h:925
proteus::cppHsuSedStress::betaForch_
double betaForch_
Definition: SedClosure.h:926
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::sigmaC_
double sigmaC_
Definition: SedClosure.h:932
proteus::cppHsuSedStress::vos_limiter_
double vos_limiter_
Definition: SedClosure.h:943
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::packFraction_
double packFraction_
Definition: SedClosure.h:928
proteus::cppHsuSedStress::eR_
double eR_
Definition: SedClosure.h:935
proteus::cppHsuSedStress::mContact_
double mContact_
Definition: SedClosure.h:937
proteus
Definition: ADR.h:17
proteus::cppHsuSedStress::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::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
Definition: SedClosure.h:12
proteus::cppHsuSedStress::k_diff
double k_diff(double sedF, double rhoSolid, double theta)
Definition: SedClosure.h:578
proteus::cppHsuSedStress::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::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_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::packMargin_
double packMargin_
Definition: SedClosure.h:929
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::notSoLarge_
double notSoLarge_
Definition: SedClosure.h:941
proteus::cppHsuSedStress::large_
double large_
Definition: SedClosure.h:942
Generated on Fri Jul 1 2022 11:17:40 for proteus by doxygen 1.8.20