Efficient resolution of the Colebrook-White equation

개요

비압축성 유체가 파이프를 흐를 때 마찰에 의한 수두손실 계산을 위한 Darcy friction factor를 구하는 VBA Code이다. Reynolds Number 2300 미만은 층류유동으로 friction factor f=64/Re로 계산이 되고, 천이 및 난류유동에서는 Colebrook의 공식 의해 Darcy friction factor를 계산할 수 있다.

 

아래 VBA Code는 D. Clamond의 논문에 언급된 Matlab Code를 참고하여 수정한 것이다. 

 

VBA CODE

'Colebrook Equation
'
'This VBA code was coded by referring to matlab code by D. Clamond
'Clamond D. Efficient resolution of the colebrook equation. Industrial & Engineering Chemistry Research 2009; 48: p. 3665-3671.
'https://nl.mathworks.com/matlabcentral/fileexchange/21990-colebrook-m?focused=5105324&tab=function

' ======================================= MATLAP CODE =======================================
'
' F = COLEBROOK(R,K) fast, accurate and robust computation of the
'     Darcy-Weisbach friction factor F according to the Colebrook equation:
'                             -                       -
'      1                     |    K        2.51        |
'  ---------  =  -2 * Log_10 |  ----- + -------------  |
'   sqrt(F)                  |   3.7     R * sqrt(F)   |
'                             -                       -
' INPUT:
'   R : Reynolds' number (should be >= 2300).
'   K : Equivalent sand roughness height divided by the hydraulic diameter (default K=0).
'
' OUTPUT:
'   F : Friction factor.
'
' FORMAT:
'   R, K and F are either scalars or compatible arrays.
'
' ACCURACY:
'   Around machine precision forall R > 3 and forall 0 <= K,
'   i.e. forall values of physical interest.
'
' EXAMPLE: F = colebrook([3e3,7e5,1e100],0.01)
'
' Edit the m-file for more details.
' Method: Quartic iterations.
' Reference: http://arxiv.org/abs/0810.5564
' Read this reference to understand the method and to modify the code.
' Author: D. Clamond, 2008-09-16.
' Check for errors.
'if any(R(:)<2300) == 1,
'   warning('The Colebrook equation is valid for Reynolds'' numbers >= 2300.');
'end,
'if nargin == 1 || isempty(K) == 1,
'   K = 0;
'end,
'if any(K(:)<0) == 1,
'   warning('The relative sand roughness must be non-negative.');
'end,
'% Initialization.
'X1 = K .* R * 0.123968186335417556;              % X1 <- K * R * log(10) / 18.574.
'X2 = log(R) - 0.779397488455682028;              % X2 <- log( R * log(10) / 5.02 );
'% Initial guess.
'F = X2 - 0.2;
'% First iteration.
'E = ( log(X1+F) - 0.2 ) ./ ( 1 + X1 + F );
'F = F - (1+X1+F+0.5*E) .* E .*(X1+F) ./ (1+X1+F+E.*(1+E/3));
'% Second iteration (remove the next two lines for moderate accuracy).
'E = ( log(X1+F) + F - X2 ) ./ ( 1 + X1 + F );
'F = F - (1+X1+F+0.5*E) .* E .*(X1+F) ./ (1+X1+F+E.*(1+E/3));
'% Finalized solution.
'F = 1.151292546497022842 ./ F;                   % F <- 0.5 * log(10) / F;
'F = F .* F;                                      % F <- Friction factor.
' ===========================================================================================
'
Function Colebrook(R As Double, K As Double) As Double
    Dim X1 As Double, X2 As Double, F As Double, E As Double

'Laminar flow
    If R < 2300 And K < 0 Then
    Colebrook = CVErr(xlErrValue)
    Exit Function
    End If
    
    If R < 2300 And K >= 0 Then
    Colebrook = 64 / R
    End If

'Transient&Turbulent flow
    If R >= 2300 And K < 0 Then
    Colebrook = CVErr(xlErrValue)
    Exit Function
    End If

    If R >= 2300 And K >= 0 Then
    'Initialization.
    X1 = K * R * 0.123968186335418  'X1 <- K * R * log(10) / 18.574
    X2 = Log(R) - 0.779397488455682 'X2 <- log( R * log(10) / 5.02 )
    
    'Initial guess.
    F = X2 - 0.2
    
    'First iteration.
    E = (Log(X1 + F) + F - X2) / (1 + X1 + F)
    F = F - (1 + X1 + F + 0.5 * E) * E * (X1 + F) / (1 + X1 + F + E * (1 + E / 3))
    
    'Second iteration (remove the next two lines for moderate accuracy).
    E = (Log(X1 + F) + F - X2) / (1 + X1 + F)
    F = F - (1 + X1 + F + 0.5 * E) * E * (X1 + F) / (1 + X1 + F + E * (1 + E / 3))
    
    'Finalized solution.
    F = 1.15129254649702 / F 'F <- 0.5 * log(10) / F
    Colebrook = F * F        'F <- Friction factor.
    End If
End Function