Category Archives: Scilab

コールブルック式による管摩擦係数(セカント法による数値解析)、Calculation of Friction Factor by Colebrook-White

Equation
\frac{1}{\sqrt{\lambda}}=-2\log_{10}(\frac{\varepsilon/d}{3.71}+\frac{2.51}{Re\sqrt{\lambda}})

Scilab
1980119802
1980319804

 -->function Result=func(re,ed)
-->    limit=10^-8;loop=10000;n=2;
-->    deff('[y]=f(lambda)','y=1/lambda^0.5+2*log10(ed/3.71+2.51/re/lambda^0.5)');//λ(lambda):Friction Factor
-->    lambda(1)=0.005;lambda(2)=0.0051;
-->    while n <= loop
-->        lambda(n+1) = lambda(n) - f(lambda(n))*(lambda(n) - lambda(n-1))/(f(lambda(n)) - f(lambda(n-1)));
-->        if abs(f(lambda(n+1))) < limit then
-->            break;
-->        end
-->        n = n +1;
-->    end
-->    Result=lambda(n);
-->endfunction
-->for www=1:2;
-->    funcprot(0);
-->    NNN1=50;
-->    NNN2=50;
-->    equation='Result=func(re,ed)';
-->    xset("window",www);clf();xgrid();
-->    /////////////////////////////////////
-->    ed=linspace(10^-6,0.05,NNN1);//Surface Roughness(ε)/Diameter of Pipe(d)
-->    re=linspace(10^4,10^5,NNN2);//Reynolds Number
-->    graphtitle='Friction Factor:λ';
-->    if www==1 | www==2 then
-->        xtitle(graphtitle,'Re','ε/d','FF',boxed=1);
-->        fnt='[Result]=func1(re,ed)';
-->        deff(fnt,equation);
-->        if www==1 then 
-->            R=245;G=255;B=250;color(R,G,B);idcolor=color(R,G,B);
-->            fplot3d(re,ed,func1,flag=[idcolor , 2 , 4]);
-->            re=5000;ed=0.05;
-->            printf("\nResult:ε/d= %e,Re=%e,λ=%.8f \n",ed,re,func1(re,ed));
-->            re=1*10^8;ed=0.00001;
-->            printf("\nResult:ε/d= %e,Re=%e,λ=%.8f \n",ed,re,func1(re,ed));
-->        elseif www==2 then
-->            nz=30;
-->            contour(re,ed,func1,nz,flag=[2,2,4]);
-->            axe=gca();
-->            for iii=1:nz;
-->                cnt=axe.children(iii).children;
-->                eee=length(cnt);
-->                if eee==2 then
-->                    cnt(2).font_size=4;
-->                end
-->            end
-->        end
-->    else
-->    end 
-->    currentaxes=gca();
-->    currentaxes.font_size=6;
-->    currentaxes.title.font_size=6;
-->    currentaxes.x_label.font_size=6;
-->    currentaxes.y_label.font_size=6;
-->    currentaxes.z_label.font_size=6;
-->    currentaxes.tight_limits="on";
-->    currentaxes.cube_scaling="on";
-->    currentaxes.auto_clear="on";
-->    currentaxes.auto_scale="on";
-->    currentaxes.isoview="off";
-->end

Result:ε/d= 5.000000e-02,Re=5.000000e+03,λ=0.07586375 

Result:ε/d= 1.000000e-05,Re=1.000000e+08,λ=0.00818443 
-->function Result=func(re,ed)
-->    limit=10^-8;loop=10000;n=2;
-->    deff('[y]=f(lambda)','y=1/lambda^0.5+2*log10(ed/3.71+2.51/re/lambda^0.5)');//λ(lambda):Friction Factor
-->    lambda(1)=0.005;lambda(2)=0.0051;
-->    while n <= loop
-->        lambda(n+1) = lambda(n) - f(lambda(n))*(lambda(n) - lambda(n-1))/(f(lambda(n)) - f(lambda(n-1)));
-->        if abs(f(lambda(n+1))) < limit then
-->            break;
-->        end
-->        n = n +1;
-->    end
-->    Result=lambda(n);
-->endfunction
-->for www=1:2;
-->    funcprot(0);
-->    NNN1=50;
-->    NNN2=50;
-->    equation='Result=func(re,ed)';
-->    xset("window",www);clf();xgrid();
-->    /////////////////////////////////////
-->    ed=linspace(10^-6,0.05,NNN1);//Surface Roughness(ε)/Diameter of Pipe(d)
-->    re=linspace(5000,10^4,NNN2);//Reynolds Number
-->    graphtitle='Friction Factor:λ';
-->    if www==1 | www==2 then
-->        xtitle(graphtitle,'Re','ε/d','FF',boxed=1);
-->        fnt='[Result]=func1(re,ed)';
-->        deff(fnt,equation);
-->        if www==1 then 
-->            R=245;G=255;B=250;color(R,G,B);idcolor=color(R,G,B);
-->            fplot3d(re,ed,func1,flag=[idcolor , 2 , 4]);
-->            re=5000;ed=0.05;
-->            printf("\nResult:ε/d= %e,Re=%e,λ=%.8f \n",ed,re,func1(re,ed));
-->            re=1*10^8;ed=0.00001;
-->            printf("\nResult:ε/d= %e,Re=%e,λ=%.8f \n",ed,re,func1(re,ed));
-->        elseif www==2 then
-->            nz=30;
-->            contour(re,ed,func1,nz,flag=[2,2,4]);
-->            axe=gca();
-->            for iii=1:nz;
-->                cnt=axe.children(iii).children;
-->                eee=length(cnt);
-->                if eee==2 then
-->                    cnt(2).font_size=4;
-->                end
-->            end
-->        end
-->    else
-->    end 
-->    currentaxes=gca();
-->    currentaxes.font_size=6;
-->    currentaxes.title.font_size=6;
-->    currentaxes.x_label.font_size=6;
-->    currentaxes.y_label.font_size=6;
-->    currentaxes.z_label.font_size=6;
-->    currentaxes.tight_limits="on";
-->    currentaxes.cube_scaling="on";
-->    currentaxes.auto_clear="on";
-->    currentaxes.auto_scale="on";
-->    currentaxes.isoview="off";
-->end

Result:ε/d= 5.000000e-02,Re=5.000000e+03,λ=0.07586375 

Result:ε/d= 1.000000e-05,Re=1.000000e+08,λ=0.00818443 
 -->limit=10^-8;
-->loop=10000;
-->ed=8*10^-5;//Surface Roughness(ε)/Diameter of Pipe(d)
-->re=6*10^3;//Reynolds Number
-->deff('[y]=f(lambda)','y=1/lambda^0.5+2*log10(ed/3.71+2.51/re/lambda^0.5)');//λ(lambda):Friction Factor
-->n=2;
-->lambda(1)=0.005;
-->lambda(2)=0.0051;
-->while n <= loop
-->    lambda(n+1) = lambda(n) - f(lambda(n))*(lambda(n) - lambda(n-1))/(f(lambda(n)) - f(lambda(n-1)));
-->    if abs(f(lambda(n+1))) < limit then
-->        break;
-->    end
-->    n = n +1;
-->end
-->printf("\nResult:ε/d= %e,Re=%e,λ=%.8f,y=%e",ed,re,lambda(n),f(lambda(n)));

Result:ε/d= 8.000000e-05,Re=6.000000e+03,λ=0.03559990,y=1.873973e-07

SciNotes

function Result=func(re,ed)
    limit=10^-8;loop=10000;n=2;
    deff('[y]=f(lambda)','y=1/lambda^0.5+2*log10(ed/3.71+2.51/re/lambda^0.5)');//λ(lambda):Friction Factor
    lambda(1)=0.005;lambda(2)=0.0051;
    while n <= loop
        lambda(n+1) = lambda(n) - f(lambda(n))*(lambda(n) - lambda(n-1))/(f(lambda(n)) - f(lambda(n-1)));
        if abs(f(lambda(n+1))) < limit then
            break;
        end
        n = n +1;
    end
    Result=lambda(n);
endfunction
for www=1:2;
    funcprot(0);
    NNN1=50;
    NNN2=50;
    equation='Result=func(re,ed)';
    xset("window",www);clf();xgrid();
    /////////////////////////////////////
    ed=linspace(10^-6,0.05,NNN1);//Surface Roughness(ε)/Diameter of Pipe(d)
    re=linspace(5000,10^4,NNN2);//Reynolds Number
    graphtitle='Friction Factor:λ';
    if www==1 | www==2 then
        xtitle(graphtitle,'Re','ε/d','FF',boxed=1);
        fnt='[Result]=func1(re,ed)';
        deff(fnt,equation);
        if www==1 then 
            R=245;G=255;B=250;color(R,G,B);idcolor=color(R,G,B);
            fplot3d(re,ed,func1,flag=[idcolor , 2 , 4]);
            re=5000;ed=0.05;
            printf("\nResult:ε/d= %e,Re=%e,λ=%.8f \n",ed,re,func1(re,ed));
            re=1*10^8;ed=0.00001;
            printf("\nResult:ε/d= %e,Re=%e,λ=%.8f \n",ed,re,func1(re,ed));
        elseif www==2 then
            nz=30;
            contour(re,ed,func1,nz,flag=[2,2,4]);
            axe=gca();
            for iii=1:nz;
                cnt=axe.children(iii).children;
                eee=length(cnt);
                if eee==2 then
                    cnt(2).font_size=4;
                end
            end
        end
    else
    end 
    currentaxes=gca();
    currentaxes.font_size=6;
    currentaxes.title.font_size=6;
    currentaxes.x_label.font_size=6;
    currentaxes.y_label.font_size=6;
    currentaxes.z_label.font_size=6;
    currentaxes.tight_limits="on";
    currentaxes.cube_scaling="on";
    currentaxes.auto_clear="on";
    currentaxes.auto_scale="on";
    currentaxes.isoview="off";
end
limit=10^-8;
loop=10000;
ed=8*10^-5;//Surface Roughness(ε)/Diameter of Pipe(d)
re=6*10^3;//Reynolds Number
deff('[y]=f(lambda)','y=1/lambda^0.5+2*log10(ed/3.71+2.51/re/lambda^0.5)');//λ(lambda):Friction Factor
n=2;
lambda(1)=0.005;
lambda(2)=0.0051;
while n <= loop
    lambda(n+1) = lambda(n) - f(lambda(n))*(lambda(n) - lambda(n-1))/(f(lambda(n)) - f(lambda(n-1)));
    if abs(f(lambda(n+1))) < limit then
        break;
    end
    n = n +1;
end
printf("\nResult:ε/d= %e,Re=%e,λ=%.8f,y=%e",ed,re,lambda(n),f(lambda(n)));

参考文献
Graham Woan(2007).『ケンブリッジ物理公式ハンドブック』.共立出版.286pp.
板東修(2008).『Execlで解く配管とポンプの流れ』.工業調査会.157pp.

アプリケーション
URL Scilab http://www.scilab.org/