Commit 28681ed8 authored by Philippe Delandmeter's avatar Philippe Delandmeter
Browse files

slim2D: bottom sediment flux

parent a62d979a
Pipeline #1776 passed with stage
in 34 minutes and 25 seconds
......@@ -1335,3 +1335,59 @@ void getHC::operator()(functorCache &cm, fullMatrix<double> &val) const {
}
}
sedimentBottomFlux2D::sedimentBottomFlux2D(const functor *sed, const functor *uvBot, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity) {
_sed = sed;
_uvBot = uvBot;
_wind = wind;
_bath = bath;
_sedGround = sedGround;
_gravity = gravity;
}
void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val) const {
const fullMatrix<double> &sed = cm.get(*_sed);
const fullMatrix<double> &uvBot = cm.get(*_uvBot);
const fullMatrix<double> &wind = cm.get(*_wind);
const fullMatrix<double> &bath = cm.get(*_bath);
const fullMatrix<double> &sedGround = cm.get(*_sedGround);
const fullMatrix<double> &gravity = cm.get(*_gravity);
double factorW = 1e-2;
double wMax = 1e-4;
double u0 = 0.3;
double w0 = 10;
double a01 = 2.8e-2;
double a02 = 1.44e-1;
double a1 = 1e-6;
double n = 4.;
double omega2 = 2.4538;
double eros0Fact = 0.245;
val.resize(cm.nEvaluationPoint(),1,true);
for (int i = 0; i < cm.nEvaluationPoint(); i++) {
if (sed(i,0) > 3) Msg::Warning("vertical Velocity for sediment parametrization not valid anymore if sediment concentration greater than 3kg m^-3");
double uBot = hypot(uvBot(i,0), uvBot(i,1));
double wS = std::max(0., std::min(wMax, factorW * sed(i,0)));
double uWind = hypot(wind(i,0), wind(i,1));
wS = (uWind < 10) ? wS : 0;
if (sed(i,0) < -1e-9) Msg::Warning("Wrong!! Negative sediment concentration");
double sed0 = std::max(0., sed(i,0));
double settling = (uBot < u0) ? - sed0 * wS * (1.- _pow2(uBot/u0)) : 0;
if (sedGround(i,0) < 1e-6)
val(i,0) = settling;
else{
double erosion0 = 0;
if (uWind > 1e-3){
double alpha = (wind(i,0)*uvBot(i,0) + wind(i,1)*uvBot(i,1))/(uWind*uBot);
double a0 = a01 + (a02 - a01) * (0.5 - atan(10*alpha)/(2*atan(10)));
double F = eros0Fact*exp(-0.905 * omega2 * bath(i,0) / gravity(i,0) - 0.0207);
erosion0 = a0 * pow(uWind/w0, 3)*pow(uBot/u0, n)*F;
}
double erosion1 = (uBot > u0) ? a1*(pow(uBot/u0, n) - 1) : 0;
double erosionTotal = std::min(0.01*sedGround(i,0), erosion0+erosion1); // supposing dt3d < 100, there will be no more erosion than available
if (erosionTotal < 0)
Msg::Fatal("Total erosion should always be positive! dt should be used for the max erosion term");
val(i,0) = (settling + erosionTotal)/bath(i,0);
}
}
}
......@@ -385,5 +385,14 @@ class getHC: public functor {
getHC(const functor *C, const functor *bath, const functor *hydro) ;
};
class sedimentBottomFlux2D : public functor {
const functor *_sed, *_uvBot, *_wind, *_bath, *_sedGround, *_gravity;
public :
void operator()(functorCache &cm, fullMatrix<double> &val) const ;
sedimentBottomFlux2D(const functor *sed, const functor *uvBot, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity);
};
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment