井筒多相管流的Beggs and Bril计算方法C++实现

#include<iostream>
#include<numbers>
#include<algorithm>
#include<vector>
using namespace std;

// 计算平均压力
double calculateAveragePressure(double currentPressure, double pressureIncrement, double Nx) {
    return (currentPressure + pressureIncrement) / Nx;
}

// 计算平均温度
double calculateAverageTemperature(double T) {
    // 根据具体情况编写平均温度的计算公式
    return T + 273;  // 替换为实际的平均温度计算公式
}

//计算原油密度
double calculateRhoL(double Go, double Rs, double Gg, double Bo) {
    return (Go + 1.293 * Rs * Gg) / Bo;
}

//计算天然气密度
double calculateRhog(double Gg, double avgP, double Z, double T) {
    return (Gg * 1.293 * avgP * 293) / (Z * T * 0.101);
}

//计算气的就地流量
double calculateQg(double avgT, double Z, double Rp, double Rs, double Qo, double avgP) {
    return 0.101325 * avgT * Z * (Rp - Rs) * Qo / (86400 * avgP * 293);
}

//计算液的就地流量
double calculateQL(double Qo, double Bo) {
    return Qo * Bo / 86400;
}
//计算就气液混合物的表观流速
double calculateVm(double InnerD, double Qg, double QL) {
    double Ap = numbers::pi * pow(InnerD / 2, 2);
    double VsL = QL / Ap;
    double Vsg = Qg / Ap;
    return VsL + Vsg;
}
//计算就液的表观流速
double calculateVsL(double InnerD, double QL) {
    double Ap = numbers::pi * pow(InnerD / 2, 2);
    double VsL = QL / Ap;
    return VsL;

}
//计算就气的表观流速
double calculateVsg(double InnerD, double Qg) {
    double Ap = numbers::pi * pow(InnerD / 2, 2);
    double Vsg = Qg / Ap;
    return Vsg;
}

//计算液体气体的总质量流量
double calculateG(double RhoL, double Rhog, double QL, double Qg) {
    double GL = RhoL * QL;
    double Gg = Rhog * Qg;
    return GL + Gg;
}
//计算液体的总质量流量
double calculateGL(double RhoL, double Rhog, double QL, double Qg) {
    double GL = RhoL * QL;
    return GL;
}
//计算气体的总质量流量
double calculateGg(double RhoL, double Rhog, double QL, double Qg) {
    double Gg = Rhog * Qg;
    return Gg;
}
//计算入口体积含液率（无滑脱持液率）
double calculateEL(double QL, double Qg) {
    return QL / (QL + Qg);
}
//计算弗鲁德数Nfr
double calculateNFr(double Vm, double InnerD) {
    return pow(Vm, 2) / (9.81 * InnerD);
}
//计算混合物粘度，MuL实际是原油粘度
double  calculateMum(double MuL, double EL, double Mug) {
    return MuL * EL + Mug * (1 - EL);
}
//计算液体速度准数
double  calculateNvL(double VsL, double RhoL, double e) {
    return VsL * pow(RhoL / (9.81 * e), 0.25);
}
//判断流型计算
double calculateHLtheta(double EL, double NFr, double NvL, double InCl) {
    double L1 = 316 * pow(EL, 0.302);//235.93
    double L2 = 92.52 * pow(10, -5) * pow(EL, -2.4684);//0.0101
    double L3 = 0.1 * pow(EL, -1.4516);//0.41
    double L4 = 0.5 * pow(EL, -6.738);//350.058
    if ((EL < 0.01 && NFr < L1) or (EL >= 0.01 && NFr < L2)) {//分离流
        double a = 0.98;
        double b = 0.4846;
        double c = 0.0868;
        double d = 0.011;
        double e = -3.768;
        double f = 3.539;
        double g = -1.614;
        double HLo = max(a * pow(EL, b) / pow(NFr, c), EL);
        double C = (1 - EL) * log(d * pow(EL, e) * pow(NvL, f) * pow(NFr, g));
        double thetaFai = 1 + C * (sin(1.8 * InCl) - 1.0 / 3 * pow(sin(1.8 * InCl), 3));
        double HLtheta = HLo * thetaFai;
        return HLtheta;
    }
    else if ((EL >= 0.01 && EL<0.4 && NFr > L2 && NFr <= L1) or (EL >= 0.4 && NFr > L3 && NFr < L4)) {//间歇流
        double a1 = 0.845;
        double b1 = 0.5351;
        double c1 = 0.0173;
        double d1 = 2.96;
        double e1 = 0.305;
        double f1 = -0.4473;
        double g1 = 0.0978;
        double HLo1 = max(a1 * pow(EL, b1) / pow(NFr, c1), EL);//0.508
        double C1 = (1 - EL) * log(d1 * pow(EL, e1) * pow(NvL, f1) * pow(NFr, g1));//0.301
        double thetaFai1 = 1 + C1 * (sin(1.8 * InCl * numbers::pi / 180) - 1.0 / 3 * pow(sin(1.8 * InCl * numbers::pi / 180), 3));//1.2
        double HLtheta1 = HLo1 * thetaFai1;//0.61
        return HLtheta1;
    }
    else if ((EL < 0.4 && NFr >= L1) or (EL >= 0.4 && NFr > L4)) {//分散流
        double a2 = 1.065;
        double b2 = 0.5929;
        double c2 = 0.0609;
        double d2 = 4.7;
        double e2 = -0.3692;
        double f2 = 0.1244;
        double g2 = -0.5056;
        double HLo2 = max(a2 * pow(EL, b2) / pow(NFr, c2), EL);
        double C2 = (1 - EL) * log(d2 * pow(EL, e2) * pow(NvL, f2) * pow(NFr, g2));
        double thetaFai2 = 1 + C2 * (sin(1.8 * InCl) - 1.0 / 3 * pow(sin(1.8 * InCl), 3));
        double HLtheta2 = HLo2 * thetaFai2;
        return HLtheta2;
    }
    else {//过渡流
        double a = 0.98;
        double b = 0.4846;
        double c = 0.0868;
        double d = 0.011;
        double e = -3.768;
        double f = 3.539;
        double g = -1.614;
        double HLo = max(a * pow(EL, b) / pow(NFr, c), EL);
        double C = (1 - EL) * log(d * pow(EL, e) * pow(NvL, f) * pow(NFr, g));
        double thetaFai = 1 + C * (sin(1.8 * InCl) - 1.0 / 3 * pow(sin(1.8 * InCl), 3));
        double HLtheta = HLo * thetaFai;

        double a1 = 0.845;
        double b1 = 0.5351;
        double c1 = 0.0173;
        double d1 = 2.96;
        double e1 = 0.035;
        double f1 = -0.4473;
        double g1 = 0.0978;
        double HLo1 = max(a1 * pow(EL, b1) / pow(NFr, c1), EL);
        double C1 = (1 - EL) * log(d1 * pow(EL, e1) * pow(NvL, f1) * pow(NFr, g1));
        double thetaFai1 = 1 + C1 * (sin(1.8 * InCl) - 1.0 / 3 * pow(sin(1.8 * InCl), 3));
        double HLtheta1 = HLo1 * thetaFai1;

        double A = (L3 - NFr) / (L3 - L2);
        double B = 1 - A;
        double HLtheta3 = A * HLtheta + B * HLtheta1;
    }
}
//计算两相流动的雷诺数
double calculateNRe(double InnerD, double Vm, double RhoL, double EL, double Rhog, double MuL, double Mug) {
    return (InnerD * Vm * (RhoL * EL + Rhog * (1 - EL))) / (MuL * EL + Mug * (1 - EL));
}
//计算阻力系数Lambda
double calculateLambda(double EL, double HLtheta, double NRe) {
    double y = EL / pow(HLtheta, 2);//1.01
    double S = 0;
    if (y > 1 && y < 1.2) {
        S = log(2.2 * y - 1.2);//0.028916088376603331
    }
    else {
        S = log(y) / (-0.0523 + 3.18 * log(y) - 0.8725 * pow(log(y), 2) + 0.01853 * pow(log(y), 4));
    }
    double LambdaPie = 0.0056 + 0.5 / pow(NRe, 0.32);//0.25797
    double Lambda = LambdaPie * exp(S);//0.2655
    return Lambda;
}
double Tr(double K, double L, double Mu, double Re, double Rw) {
    return 2 * numbers::pi * K * L / (Mu * log(pow(Re / Rw, 0.5))) * pow(10, 5);
}
double calculatechangeP(double RhoL, double HLtheta, double Rhog, double Lambda, double G, double Vm, double InnerD, double InCl, double P, double Vsg) {
    double Ap = numbers::pi * pow(InnerD / 2, 2);//0.003管子截断面积
    double a = sin(InCl * numbers::pi / 180);
    double changeP = ((RhoL * HLtheta + Rhog * (1 - HLtheta)) * 9.81 * sin(InCl * numbers::pi / 180) + (Lambda * G * Vm) / (2 * InnerD * Ap)) / (1 - ((RhoL * HLtheta + Rhog * (1 - HLtheta)) * Vm * Vsg) / (P * pow(10, 6))) * pow(10, -6);
    return changeP;
}
vector<double> BeggsAndBrillIteration(double Qo, double Rp, double Rs, double Go, double Gg, double T, double InnerD, double L, double Pbh, int NFor, double StartPressureDrop, double Bo, double Z, double e, double MuL, double Mug, double InCl) {
    vector<double> P;
    double avgT = calculateAverageTemperature(T);//313
    double RhoL = calculateRhoL(Go, Rs, Gg, Bo);//817.3
    double QL = calculateQL(Qo, Bo);//0.000617
    double VsL = calculateVsL(InnerD, QL);//0.20
    double NvL = calculateNvL(VsL, RhoL, e);//1.68
    double error = std::numeric_limits<double>::max(); // 初始化误差
    const double tolerance = 0.1; // 设置收敛阈值
    double currentPressure = 0;
    double newPressure = 0;
    for (int i = 0; i < NFor; i++) {
        currentPressure = Pbh + StartPressureDrop;//p21
        double avgP = calculateAveragePressure(currentPressure, Pbh, 2);//4.25
        double Rhog = calculateRhog(Gg, avgP, Z, avgT);//39.65
        double Qg = calculateQg(avgT, Z, Rp, Rs, Qo, avgP);//0.001
        double Vm = calculateVm(InnerD, Qg, QL);//0.54
        double Gm = calculateG(RhoL, Rhog, QL, Qg);//0.54
        double EL = calculateEL(QL, Qg);//0.38
        double NFr = calculateNFr(Vm, InnerD);//0.479
        double Mum = calculateMum(MuL, EL, Mug);//1.31
        double Vsg = calculateVsg(InnerD, Qg);//0.34
        double HLtheta = calculateHLtheta(EL, NFr, NvL, InCl);//0.61
        double NRe = calculateNRe(InnerD, Vm, RhoL, EL, Rhog, MuL, Mug);//8.47
        double Lambda = calculateLambda(EL, HLtheta, NRe);//0.2655
        double changeP = calculatechangeP(RhoL, HLtheta, Rhog, Lambda, Gm, Vm, InnerD, InCl, avgP, Vsg) * L;//0.37793190750543698
        newPressure = Pbh + changeP;//p22
        
        while (abs(newPressure - currentPressure) >= tolerance) {
            currentPressure = newPressure;
            avgP = calculateAveragePressure(currentPressure, Pbh, 2);//4.25
            Rhog = calculateRhog(Gg, avgP, Z, avgT);//39.65
            Qg = calculateQg(avgT, Z, Rp, Rs, Qo, avgP);//0.001
            Vm = calculateVm(InnerD, Qg, QL);//0.54
            Gm = calculateG(RhoL, Rhog, QL, Qg);//0.54
            EL = calculateEL(QL, Qg);//0.38
            NFr = calculateNFr(Vm, InnerD);//0.479
            Mum = calculateMum(MuL, EL, Mug);//1.31
            Vsg = calculateVsg(InnerD, Qg);//0.34
            HLtheta = calculateHLtheta(EL, NFr, NvL, InCl);//0.61
            NRe = calculateNRe(InnerD, Vm, RhoL, EL, Rhog, MuL, Mug);//8.47
            Lambda = calculateLambda(EL, HLtheta, NRe);//0.2655
            changeP = calculatechangeP(RhoL, HLtheta, Rhog, Lambda, Gm, Vm, InnerD, InCl, avgP, Vsg) * L;//0.37793190750543698
            newPressure = Pbh + changeP;//p22
        }
        cout <<"误差"<<abs(newPressure - currentPressure) << endl;
        P.emplace_back(newPressure);
        Pbh = newPressure;       
    }
    return P;
}


int main() {
   double Qo = 150;//不含水自喷井油井产量m3/day
   double Rp = 100;//生产气油比m3/m3
   double Rs = 54;//溶解气油比m3/m3
   double Go = 850;//地面脱气原油密度kg/m3
   double Gg = 0.7;//标准状况下气密度kg/m3
   double T = 57;//平均温度°C
   double InnerD = 0.062;//油管内径m
   double L = 200;//计算段长m
   double Pbh = 0.7;//起始点压力MPa
   double NFor = 100;//分段数
   double StartPressureDrop = 0.5;//初设压降Mpa
   double Bo = 1.066;//原油体积系数
   double Z = 0.899;//天然压缩因子
   double e = 0.0183;//原油表面张力N/m
   double MuL = 3.47;//原油粘度mPa·s
   double Mug = 0.0122;//天然气粘度mPa·s
   double InCl = 45;//角度
    auto BeggsAndBrillIterations = BeggsAndBrillIteration(Qo, Rp, Rs, Go, Gg, T, InnerD, L, Pbh, NFor, StartPressureDrop, Bo, Z, e, MuL, Mug, InCl);
    for (auto i : BeggsAndBrillIterations) {
        cout << "P=" << i<< "Mpa" << endl;
    }
    return 0;
}

参考书籍提取码xfsr
·流量卡·可以看看29元 235G~