conductor_tension/core.cs

using System;
namespace conductortension
{

    public partial class Core
    {

        //解状态方程的牛顿法
        public static double nt_equation(double A, double B, double sigma)
        {
            return Math.Pow(sigma, 3) + A * Math.Pow(sigma, 2) - B;
        }


        public static double d_eqution(double A, double sigma)
        {
            return 3 * Math.Pow(sigma, 2) + 2 * A * sigma;
        }



        public static double newton(double A, double B, double sigma,ref bool succeed)
        {
            succeed = false;
            double _sigma = sigma;
            int i;
            for (i=0; i < 20; i++)
            {
                double fx = nt_equation(A, B, _sigma);
                double d_fx = d_eqution(A, _sigma);
                double d_sigma = -fx / d_fx;
                if (Math.Abs(d_sigma) < 1e-5)
                {
                    succeed = true;
                    break;
                }

                _sigma = _sigma + d_sigma;
            }
            return _sigma;
        }



        //从状态方程接触张力
        //start_tension  牛顿法初值
        //alpha 是膨胀系数
        //E 弹性模量
        //area 导线截面
        //conrol_t 控制条件温度
        //conrol_load 控制条件荷载
        //conrol_tension 控制条件张力
        //load 计算工况荷载
        //t 计算工况温度
        //l 档距
        public static double tensionFrmoStateEquation(double start_tension, double E, double alpha, double area, double control_load, double control_tension, double control_t, double load, double t, double l)
        {
            double A, B;
            //A = E / 24 * Math.Pow(control_load * l / control_tension, 2) - control_tension / area + alpha * E * (t - control_t);
            A = Fx(E, control_load, control_tension, alpha, area, control_t, l) + alpha * E * t;
            B = E * Math.Pow(load/area * l, 2) / 24;
            double tension;//牛顿法初值
            bool succeed=false;
            tension = newton(A, B, start_tension / area, ref succeed) * area;
            if (succeed)
            {
                return tension;
            }
            return double.NaN;    
        }
        //alpha 是膨胀系数
        //E 弹性模量
        //area 导线截面
        //t 温度
        // l 档距
        public static double Fx(double E, double load, double tension, double alpha, double area, double t, double l)
        {
            return E * Math.Pow(load * l / tension, 2) / 24 - tension / area - alpha * E * t;
        }

        public static bool CalZhangLi(calStressLoadNew inPut, ref double ZhangLi)
        {
            //张力=T6
            //计算最大使用张力
            double ZuiDaShiYongZhangLi;	//最大使用张力=计算拉断力*保证破断张力取计算拉断力的/安全系数
            ZuiDaShiYongZhangLi = inPut.EDingLaDuanLi * inPut.BaoZhengPoDuanZhangLiXiShu / 100 / inPut.AnQuan;
            //ZuiDaShiYongZhangLi = inPut.LaDuanLi / inPut.AnQuan;
            //计算温度线膨胀系数
            double XianPengZhang = inPut.XianPengZhang;
            //计算平均运行张力
            double PingJunYunXingZhangLi = inPut.EDingLaDuanLi * inPut.BaoZhengPoDuanZhangLiXiShu / 100 * inPut.NianPingJunXishu/100;//年平均单位是%
            //double a, b, T1, T2, T3, T4, T5, T6;
            //计算待求工况比载
            double DaiHeZai = 0;	//待求工况比载
            if (CalHeZai(inPut, ref DaiHeZai) == false)
            {
                return false;
            }
            //临界档距表计算
            //一、初始计算数组,5项分别代表:1:气温，2:荷载，3:容许张力，4:P/T，5:Fo=-(T+A*E*a*t)(N)
            double[] ChuShiZuiDiQiWen = new double[5];		//最低温
            double[] ChuShiNianPingJun = new double[5];	//年平均气温
            double[] ChuShiDaFeng = new double[5];			//大风
            double[] ChuShiFuBing = new double[5];			//最大覆冰
            //1:存储各工况温度
            ChuShiZuiDiQiWen[0] = inPut.DiWen;
            ChuShiNianPingJun[0] = inPut.PingWen;
            ChuShiDaFeng[0] = inPut.FengWen;
            ChuShiFuBing[0] = inPut.BingWen;
            //定义各个工况结构体，将每个结构体中的待求工况设置为对应工况数据，因为计算荷载时前面写的计算荷载函数中调取的是待求工况的温度等信息，所以需要对应各个工况
            calStressLoadNew ZuiDiWen = inPut;
            calStressLoadNew NianPingJun = inPut;
            calStressLoadNew DaFeng = inPut;
            calStressLoadNew FuBing = inPut;
            ZuiDiWen.DaiBing = ZuiDiWen.DiBing;
            ZuiDiWen.DaiFeng = ZuiDiWen.DiFeng;
            ZuiDiWen.DaiZheSuanFengSu = false;
            ZuiDiWen.DaiWen = ZuiDiWen.DiWen;
            NianPingJun.DaiBing = NianPingJun.PingBing;
            NianPingJun.DaiFeng = NianPingJun.PingFeng;
            NianPingJun.DaiZheSuanFengSu = false;
            NianPingJun.DaiWen = NianPingJun.PingWen;
            DaFeng.DaiBing = DaFeng.FengBing;
            //double wind_velocity;
            //wind_velocity = Math.Round(DaFeng.FengFeng * Math.Pow(DaFeng.PingJunGaoDu / DaFeng.SheJiFengSuGaoDu, 0.16), 2);
            DaFeng.DaiFeng = DaFeng.FengFeng;
            DaFeng.DaiZheSuanFengSu = true;
            DaFeng.DaiWen = DaFeng.FengWen;
            FuBing.DaiBing = FuBing.BingBing;
            FuBing.DaiFeng = FuBing.BingFeng;
            FuBing.DaiZheSuanFengSu = false;
            FuBing.DaiWen = FuBing.BingWen;
            //2:计算并存储各个工况荷载
            CalHeZai(ZuiDiWen, ref ChuShiZuiDiQiWen[1]);
            CalHeZai(NianPingJun, ref ChuShiNianPingJun[1]);
            CalHeZai(DaFeng, ref ChuShiDaFeng[1]);
            CalHeZai(FuBing, ref ChuShiFuBing[1]);
            //3:存储各工况容许张力
            ChuShiZuiDiQiWen[2] = ZuiDaShiYongZhangLi;
            ChuShiNianPingJun[2] = PingJunYunXingZhangLi;
            ChuShiDaFeng[2] = ZuiDaShiYongZhangLi;
            ChuShiFuBing[2] = ZuiDaShiYongZhangLi;
            //4:计算并存储P/T  荷载/张力
            ChuShiZuiDiQiWen[3] = ChuShiZuiDiQiWen[1] / ChuShiZuiDiQiWen[2];
            ChuShiNianPingJun[3] = ChuShiNianPingJun[1] / ChuShiNianPingJun[2];
            ChuShiDaFeng[3] = ChuShiDaFeng[1] / ChuShiDaFeng[2];
            ChuShiFuBing[3] = ChuShiFuBing[1] / ChuShiFuBing[2];

            //计算4个工况下的Fx by 杜孟远2020.4.17
            //t_dangju 试探最大Fx时使用的档距
            int[] MaxFxConditionFlag = new int[inPut.DangJu + 1];//记录不同t_dangju下，哪种工况的Fx最大。1代表最低温，2代表年平均温，3代表大风，4代表覆冰。0位置元素无效。
            MaxFxConditionFlag[0] = 1;
            for (int t_dangju = inPut.DangJu; t_dangju <= inPut.DangJu; t_dangju++)
            {
                double Fx_dangjuZuiDiWen;
                double Fx_dangjuNianPingJun;
                double Fx_dangjuDaFeng;
                double Fx_dangjuFuBing;
                Fx_dangjuZuiDiWen = Fx(inPut.MoLiang, ChuShiZuiDiQiWen[1], ChuShiZuiDiQiWen[2], inPut.XianPengZhang, inPut.JieMianJi, ChuShiZuiDiQiWen[0], t_dangju);
                Fx_dangjuNianPingJun = Fx(inPut.MoLiang, ChuShiNianPingJun[1], ChuShiNianPingJun[2], inPut.XianPengZhang, inPut.JieMianJi, ChuShiNianPingJun[0], t_dangju);
                Fx_dangjuDaFeng = Fx(inPut.MoLiang, ChuShiDaFeng[1], ChuShiDaFeng[2], inPut.XianPengZhang, inPut.JieMianJi, ChuShiDaFeng[0], t_dangju);
                Fx_dangjuFuBing = Fx(inPut.MoLiang, ChuShiFuBing[1], ChuShiFuBing[2], inPut.XianPengZhang, inPut.JieMianJi, ChuShiFuBing[0], t_dangju);
                double maxFx;
                maxFx = Fx_dangjuZuiDiWen;
                MaxFxConditionFlag[t_dangju] = 1;
                if (maxFx < Fx_dangjuNianPingJun)
                {
                    maxFx = Fx_dangjuNianPingJun;
                    MaxFxConditionFlag[t_dangju] = 2;
                }
                if (maxFx < Fx_dangjuDaFeng)
                {
                    maxFx = Fx_dangjuDaFeng;
                    MaxFxConditionFlag[t_dangju] = 3;
                }
                if (maxFx < Fx_dangjuFuBing)
                {
                    maxFx = Fx_dangjuFuBing;
                    MaxFxConditionFlag[t_dangju] = 4;
                }
                double currentTension = 0;
                double currentLoad = 0;
                double currentTemperature = 0;
                switch (MaxFxConditionFlag[t_dangju])
                {
                    case 1:
                        currentLoad = ChuShiZuiDiQiWen[1];
                        currentTension = ChuShiZuiDiQiWen[2];
                        currentTemperature = ChuShiZuiDiQiWen[0];
                        break;
                    case 2:
                        currentLoad = ChuShiNianPingJun[1];
                        currentTension = ChuShiNianPingJun[2];
                        currentTemperature = ChuShiNianPingJun[0];
                        break;
                    case 3:
                        currentLoad = ChuShiDaFeng[1];
                        currentTension = ChuShiDaFeng[2];
                        currentTemperature = ChuShiDaFeng[0];
                        break;
                    case 4:
                        currentLoad = ChuShiFuBing[1];
                        currentTension = ChuShiFuBing[2];
                        currentTemperature = ChuShiFuBing[0];
                        break;
                    default:
                        break;
                }
                //TODO: 目前不进入这个判断 by 杜孟远 2020.4.18
                MaxFxConditionFlag[t_dangju] = MaxFxConditionFlag[t_dangju - 1];//TODO: 以后删掉这句 by 杜孟远 2020.4.18
                if (MaxFxConditionFlag[t_dangju] != MaxFxConditionFlag[t_dangju - 1])//有突变，找到了临界档距的范围 
                {
                    double criticalSpan;
                    double formerTension = 0;
                    double formerLoad = 0;
                    double formerTemperature = 0;
                    switch (MaxFxConditionFlag[t_dangju - 1])
                    {
                        case 1:
                            formerLoad = ChuShiZuiDiQiWen[1];
                            formerTension = ChuShiZuiDiQiWen[2];
                            formerTemperature = ChuShiZuiDiQiWen[0];
                            break;
                        case 2:
                            formerLoad = ChuShiNianPingJun[1];
                            formerTension = ChuShiNianPingJun[2];
                            formerTemperature = ChuShiNianPingJun[0];
                            break;
                        case 3:
                            formerLoad = ChuShiDaFeng[1];
                            formerTension = ChuShiDaFeng[2];
                            formerTemperature = ChuShiDaFeng[0];
                            break;
                        case 4:
                            formerLoad = ChuShiFuBing[1];
                            formerTension = ChuShiFuBing[2];
                            formerTemperature = ChuShiFuBing[0];
                            break;
                        default:
                            break;
                    }
                    double t1 = 24 / inPut.MoLiang * (formerTension - currentTension) / inPut.JieMianJi;
                    double t2 = 24 * inPut.XianPengZhang * (formerTemperature - currentTemperature);
                    double t3 = Math.Pow(formerLoad / formerTension, 2) - Math.Pow(currentLoad / currentTension, 2);
                    criticalSpan = Math.Sqrt((t1 + t2) / t3);//TODO:暂时不用
                    //Console.WriteLine("{0}", criticalSpan);

                }
                double start_tension = inPut.EDingLaDuanLi / inPut.AnQuan;
                ZhangLi = tensionFrmoStateEquation(start_tension,inPut.MoLiang, inPut.XianPengZhang, inPut.JieMianJi, currentLoad, currentTension, currentTemperature, DaiHeZai, inPut.DaiWen, inPut.DangJu);
            }

            ZhangLi = BaoLiuLiangWei(ZhangLi);
            return true;
        }

        public struct calStressLoadNew
        {
            public double WaiJing;                 //电线外径（mm）
            public double EDingLaDuanLi;                //电线计算拉断力（N）
            public double ZhongLiang;              //电线单位重量（kg/km）
            public double JieMianJi;               //电线截面积（mm2）
            public double MoLiang;                 //电线弹性模量（N/mm2）
            public double XianPengZhang;           //电线线性膨胀系数（10e-6/C）
            public double AnQuan;                  //安全系数
            public double NianPingJunXishu;             //年平均运行系数
            //public double XinXian;                 //新线系数
            //public double FuBingGuoZai;             //覆冰过载系数
            public double DiBing;                  //最低温工况覆冰（mm）
            public double DiFeng;                  //最低温工况风速（m/s）
            public double DiWen;                   //最低温工况气温（C）
            public double PingBing;                //年平均温工况覆冰（mm）
            public double PingFeng;                //年平均温工况风速（m/s）
            public double PingWen;                 //年平均温工况气温（C）
            public double BingBing;                //最大覆冰工况覆冰（mm）
            public double BingFeng;                //最大覆冰工况风速（m/s）
            public double BingWen;                 //最大覆冰工况气温（C）
            public double FengBing;                //最大风工况覆冰（mm）
            public double FengFeng;                //最大风工况风速（m/s）
            public double FengWen;                 //最大风工况气温（C）
            public double PingJunGaoDu;            //导线平均高度（m）	20
            public int CuCaoDu;                 //地面粗糙度类别（1-田野、乡村、山丘, 2-海岛, 3-城市）
            public double CuCaoDu_K;               //地面粗糙度系数K
            public double CuCaoDu_a;               //地面粗糙度指数a 
            //public double TiXingXiShu;             //电线体形系数			1.2
            //public double FengYaBuJunYun;          //风压不均匀系数（不需要可以删掉，计算得来）		0.75
            public double DaiBing;                 //待求工况覆冰（mm）
            public double DaiFeng;                 //待求工况风速（m/s）
            public double DaiWen;                  //待求工况气温（C）
            public int DangJu;                  //档距（m）
            public double SheJiFengSuGaoDu;		    //设计基准风速高度
            public double BaoZhengPoDuanZhangLiXiShu;	//保证破断张力取计算拉断力的%，即张力系数
            public double dxLoad;                   //线荷载
            public bool DaiZheSuanFengSu;//是否对待计算工况风速进行折算。
        }

        public static bool CalHeZai(calStressLoadNew inPut, ref double HeZai)
        {
            //  //综合荷载=Math.Pow(Math.Pow(垂直荷载,2)+Math.Pow(水平荷载,2),0.5)
            double ChuiZhiHeZai = 0;		//垂直荷载=自重荷载+9.80665*0.9*PI*待求工况覆冰厚度*(待求工况覆冰厚度+导线直径)/1000
            double ZiZhongHeZai = 0;		//自重荷载=导线单位重量*9.80665/1000/
            double bi = 1;//覆冰增大系数B
            //导线力学计算不考虑覆冰增大系数 by 杜孟远 2020.04.06
            double u;//体型系数

            if (inPut.DaiBing > 0 || inPut.WaiJing < 17)
                u = 1.2;
            else
                u = 1.1;

            ZiZhongHeZai = inPut.ZhongLiang * 9.80665 / 1000;//'导线自荷载

            ChuiZhiHeZai = 9.80665 * 0.9 * Math.PI * inPut.DaiBing * (inPut.DaiBing + inPut.WaiJing) / 1000; //'导线冰荷载

            double a = 1;//定义风压不均匀系数a
            double wind_v = inPut.DaiFeng;//计算荷载用风速
            if (inPut.DaiZheSuanFengSu)
            {
                wind_v = Math.Round(inPut.DaiFeng * Math.Pow(inPut.PingJunGaoDu / inPut.SheJiFengSuGaoDu, 0.16), 2);
            }
            double dx3 = 0;
            if (inPut.DaiFeng < 20)//要用没折算过的风速
                a = 1;
            else if (inPut.DaiFeng < 27)
                a = 0.85;
            else if (inPut.DaiFeng < 31.5)
                a = 0.75;
            else
                a = 0.7;
            if (inPut.DaiBing == 0) //'无冰工况
            {
                dx3 = Math.Pow(wind_v, 2) * inPut.WaiJing * u * a * bi / 1600;//  '导线无冰风荷载
            }
            else
            {
                dx3 = u * Math.Pow(wind_v, 2) * (inPut.WaiJing + 2 * inPut.DaiBing) * a * bi / 1600; //  '导线有冰风荷载
            }

            HeZai = Math.Sqrt(Math.Pow((ZiZhongHeZai + ChuiZhiHeZai), 2) + Math.Pow(dx3, 2));//  '存储导线荷载




            return true;
        }

    }


}