12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394 |
- #include "physics.hpp"
- //------------------------
- // Brooks and Corey model
- //------------------------
- namespace Physics{
- double g=9.81;
- double rho=1000;
- double phi=0.3;
- double k0=3e-5;
- double nivrivsat=0.01;
- double (*s)(double)=&s_BC;
- double (*ds)(double)=&ds_BC;
- void (*s_ds)(double,double&,double&)=&s_ds_BC;
- double (*kr)(double)=&kr_BC;
- double (*dkr)(double)=&dkr_BC;
- void (*kr_dkr)(double,double&,double&)=&kr_dkr_BC;
- double model_data[6]={-2000,0,3,11,0,0};
-
- #define Psat model_data[0]
- #define sres model_data[1]
- #define lambda model_data[2]
- #define alpha model_data[3]
- void setModel(Model model){
- switch(model){
- case BrooksCorey:
- s=&s_BC;
- ds=&ds_BC;
- s_ds=&s_ds_BC;
- kr=&kr_BC;
- dkr=&dkr_BC;
- kr_dkr=&kr_dkr_BC;
- Psat=-2000;
- sres=0;
- lambda=3;
- alpha=11;
- break;
- default:
- assert(false);
- };
- }
-
-
- double
- s_BC(double P){
- if(P>=Psat) return 1;
- return sres+(1-sres)*pow(Psat/P,lambda);
- }
-
- double
- ds_BC(double P){
- if(P>=Psat) return 0;
- return ((sres-1)*lambda*pow(Psat/P,lambda))/P;
- }
-
- void
- s_ds_BC(double P,double& v,double& dv){
- if(P>=Psat){
- v=1;
- dv=0;
- }
- else{
- double t=(1-sres)*pow(Psat/P,lambda);
- v=sres+t;
- dv=-(lambda*t)/P;
- }
- }
-
- double
- kr_BC(double P){
- if(P>=Psat) return 1;
- return pow(Psat/P,alpha);
- }
- double
- dkr_BC(double P){
- if(P>=Psat) return 0;
- return -alpha*pow(Psat/P,alpha)/P;
- }
- void
- kr_dkr_BC(double P,double& v,double& dv){
- if(P>=Psat){
- v=1;
- dv=0;
- }
- else{
- double t=pow(Psat/P,alpha);
- v=t;
- dv=-(alpha*t)/P;
- }
- }
- }
|