Instruction FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 30, N o 1, March 2017, pp. 81 - 91 DOI: 10.2298/FUEE1701081V ON THE NUMERICAL COMPUTATION OF CYLINDRICAL CONDUCTOR INTERNAL IMPEDANCE FOR COMPLEX ARGUMENTS OF LARGE MAGNITUDE * Slavko Vujević, Dino Lovrić University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, Split, Croatia Abstract. In this paper a numerical algorithm for computation of per-unit-length internal impedance of cylindrical conductors under complex arguments of large magnitude is presented. The presented algorithm either numerically solves the scaled exact formula for internal impedance or employs asymptotic approximations of modified Bessel functions when applicable. The formulas presented can be used for computation of per-unit-length internal impedance of solid cylindrical conductors as well as tubular cylindrical conductors. Key words: internal impedance, modified Bessel functions, large function arguments, scaling. 1. INTRODUCTION Internal impedance per-unit-length (pul) or surface impedance of cylindrical conductors is required in analysis of numerous electromagnetic problems [1-5]. This pul internal impedance can be computed using various formulas which contain special functions such as Bessel functions and modified Bessel functions [6]. Whatever formula is employed the results are valid only for smaller function arguments whereas for larger function arguments stability issues often occur. These issues are directly connected with computing special functions (Bessel functions and modified Bessel functions) under large parameters which in some cases yield extremely large values and in some cases extremely low values. In addition, these extreme values are multiplied, divided, subtracted and added which considerably makes thing worse. In this paper an algorithm is presented which circumvents the mentioned issues by first scaling the employed formulas to avoid overflow/underflow issues and then solving the expressions for modified Bessel functions in two ways - either by numerical integration or by using asymptotic approximations when applicable [7].  Received February 25, 2016; received in revised form April 7, 2016 Corresponding author: Slavko Vujević University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, Split, Croatia (E-mail: vujevic@fesb.hr) * An earlier version of this paper was presented at the 12 th International Conference on Applied Electromagnetics (ПЕС 2015), August 31 - September 2, 2015, in Niš, Serbia [1]. 82 S. VUJEVIĆ, D. LOVRIĆ The formulas presented in the paper are applicable to solid and tubular cylindrical conductors. All presented formulas are for a tubular cylindrical conductor, but by introducing the value zero for internal radius of the tubular cylindrical conductor, the pul internal impedance of a solid cylindrical conductor can be obtained. This model for computing pul internal impedance of single-layer tubular conductors represent a basis for a more general model which will be able to compute pul internal impedance of a multi- layered tubular conductor which is currently in development. 2. FORMULA FOR COMPUTATION OF TUBULAR CYLINDRICAL CONDUCTOR INTERNAL IMPEDANCE Computation of pul internal impedance of tubular cylindrical conductors (Fig. 1), which takes the skin effect into account but ignores the proximity effect, can be performed using various formulas based on different special functions. It has been concluded in the previous work of the authors of this paper that, from the numerical stability standpoint, the most suitable formula for computation of pul internal impedance of tubular conductors is based on modified Bessel functions of the first and second kind [7]: 1 0 0 1 1 1 1 1 ( ) ( ) ( ) ( ) 2 ( ) ( ) ( ) ( ) i e e i e i e e i K r I r K r I r Z r K r I r K r I r                              (1) exp (1 ) 4 j j                    (2) where σ is the electrical conductivity of the conductor material, re is the external radius of the conductor, ri is the internal radius of the conductor, 0I and 1I are complex-valued modified Bessel function of the first kind of order zero and one, 0K and 1K are complex-valued modified Bessel function of the second kind of order zero and one (also called Kelvin functions),  is the complex wave propagation constant, α is the attenuation constant, µ is the permeability of the conductor material, ω is the circular frequency and j is the imaginary unit. Fig. 1 Cross-section of a tubular cylindrical conductor As it has been shown in [7] by rearranging formula (1) and scaling it by an appropriate factor, the following formula for pul internal impedance of tubular conductors can be obtained: Numerical Computation of Cylindrical Conductor Internal Impedance 83 0 0 0 1 1 1 1 1 1 1 ( ) ( ) exp[ 2 ( )] ( ) ( ) ( ) 2 ( ) ( ) ( ) exp[ 2 ( )] ( ) ( ) s s i e e is s s e i e s s s e e i e e is s i e K r K r r r I r I r I r Z r I r K r K r r r I r I r                                           (3) where the scaled modified Bessel functions are: ( ) exp( ) ( ) s n n I r r I r        (4) ( ) exp( ) ( ) s n n K r r K r       (5) Modified Bessel functions of the first kind are scaled down exp( )r  times whereas modified Bessel functions of the second kind are scaled up exp( )r  times. In such a way quantities of similar magnitudes are obtained which consequently enables more stable computation. The computation of internal impedance Z can be further simplified depending on the magnitude of ( ) e i r r   . Numerical analysis has shown that for ( ) 19 e i r r    computation of Z must be performed using (3) in order to maintain high accuracy. However, for larger magnitudes of ( ) e i r r   simplifications of formula (3) can be performed without loss of accuracy. The following relation presents these simplifications and their interval of applicability: 0 1 15 15 ( ) 19 ( ) 10 2 ( ) ( ) 10 2 s e e is e e e i e I r ; r r r I r Z ; r r r                               (6) As can be seen from (3) and (6), it is imperative to compute scaled modified Bessel functions of the first and second kind as accurately as possible. The proposed numerical procedure for achieving this is addressed in the following section of the paper. 3. COMPUTATION OF SCALED MODIFIED BESSEL FUNCTIONS In the developed algorithm for function parameters α∙r ≤ 25 integral representation of scaled modified Bessel functions of the first and second kind is used. Integral representation of modified Bessel functions is more suitable than the infinite sum representation because the scaling factors given in (4-5) can be easily included in the integral representation of modified Bessel functions. This is not the case when using the infinite sum representation. Integrals that occur in modified Bessel functions of the first and second kind are solved numerically using adaptive Simpson rule. On the other hand, for function parameters α∙r > 25 computation of scaled modified Bessel functions of the first and second kind is performed using asymptotic approximations. Through extensive numerical analysis it has been found that for function parameter values larger than 25, asymptotic approximations of modified Bessel functions produce results of equal accuracy as the numerical solution of integral representation of modified Bessel functions but in less computation time. 84 S. VUJEVIĆ, D. LOVRIĆ 3.1. Computation of scaled modified Bessel functions of the first kind for α∙r ≤ 25 Modified Bessel function of the first kind of order zero in its integral form can be expressed by the following equation [8]: / 2 0 0 0 2 ( ) cos( sin ) exp( ) ( ) s I r j r d r I r                      (7) Further simplification of the previous expression and separation of real and imaginary parts yields the following relation for scaled modified Bessel function of the first kind of order zero: / 2 0 0 / 2 0 1 ( ) [exp( ) cos exp( ) cos ] [exp( ) sin exp( ) sin ] s I r a a b b d j a a b b d                         (8) where a and b are given by: (sin 1)a r     (9) (sin 1)b r     (10) The separation of the real and imaginary parts is performed because these integrals are solved separately using adaptive Simpson numerical integration. Numerical integration yields highly accurate results because the separated functions are simple to integrate as can be seen from Fig. 2 and Fig. 3 which depict how the real and imaginary parts of equation (8) behave on the integration interval for various values of parameter α∙r. Fig. 2 Real part of scaled modified Bessel function of the first kind of order zero for various values of parameter α∙r Numerical Computation of Cylindrical Conductor Internal Impedance 85 Fig. 3 Imaginary part of scaled modified Bessel function of the first kind of order zero for various values of parameter α∙r Integral representation of modified Bessel function of the first kind of order one can be expressed by the following equation [8]: / 2 1 0 1 2 ( ) sin( sin ) sin exp( ) ( ) s I r j j r d r I r                         (11) As before, by simplification of expression (11) and separation of real and imaginary parts, the following relation for scaled modified Bessel function of the first kind of order one can be obtained: / 2 1 0 / 2 0 1 ( ) [exp( ) cos exp( ) cos ] sin [exp( ) sin exp( ) sin ] sin s I r a a b b d j a a b b d                             (12) Two integrals present in equation (12) are again solved numerically using adaptive Simpson rule. Fig. 4 and Fig. 5 depict how the real and imaginary parts of equation (12) behave on the integration interval for various values of parameter α∙r. 86 S. VUJEVIĆ, D. LOVRIĆ Fig. 4 Real part of scaled modified Bessel function of the first kind of order one for various values of parameter α∙r Fig. 5 Imaginary part of scaled modified Bessel function of the first kind of order one for various values of parameter α∙r 3.2. Computation of scaled modified Bessel functions of the first kind for α∙r > 25 Asymptotic approximation of scaled modified Bessel function of the first kind can be expressed by [8]: 2 2 1 1 [4 (2 1) ] 1 ( ) ~ 1 ( 1) ; 0, 1 ! (8 )2 m s m t n m m n t I r n m rr                                 (13) Numerical Computation of Cylindrical Conductor Internal Impedance 87 From the previous expression asymptotic approximations of scaled modified functions of the first kind of orders zero and one can easily be deduced: 0 1 1 ( ) ~ 1 ( )2 NA s m m m c I r rr                 (14) 1 1 1 ( ) ~ 1 ( )2 NA s m m m d I r rr                 (15) where:                r r r r r NA 10000for3 10000300for5 300100for7 10050for9 5025for12 (16) 2 1 ( 1) [ (2 1) ] 8 ! m m m m t c t m          (17) 1 2 1 ( 1) [4 (2 1) ] 8 ! m m m m t d t m           (18) The expressions for cm and dm are deduced from (13) and are also used for asymptotic approximations of modified Bessel functions of the second kind. Values of NA have been determined through numerical analysis. 3.3. Computation of scaled modified Bessel functions of the second kind for α∙r ≤ 25 Integral present in the expression for the modified Bessel function of the second kind of order zero has an upper integral limit that tends to infinity [8]. Fortunately, the integral function rapidly tends to zero as the function argument increases so the infinite limit can be substituted with a finite limit tm0 without loss of accuracy: 0 0 0 0 ( ) exp( cosh ) exp( cosh ) mt K r r t dt r t dt                (19)          r tm 65 1cosh 1 0 (20) Now the scaled modified Bessel function of the second kind of order zero can be deduced from (19): 0 0 0 0 0 ( ) exp( ) cos exp( ) sin m mt t s K r d d dt j d d dt            (21) (cosh 1)d r t    (22) 88 S. VUJEVIĆ, D. LOVRIĆ The two integrals present in equation (21) are again solved numerically using adaptive Simpson rule with high accuracy. Fig. 6 and Fig. 7 depict how the real and imaginary parts of equation (21) behave on the integration interval for various values of parameter α∙r. Fig. 6 Real part of scaled modified Bessel function of the second kind of order zero for various values of parameter α∙r Fig. 7 Imaginary part of scaled modified Bessel function of the second kind of order zero for various values of parameter α∙r Similarly as for the modified Bessel function of second kind of order zero, the integral present in the expression for modified Bessel function of second kind of order one [8] can be replaced with a finite limit tm1: 1 1 0 0 ( ) exp( cosh ) cosh exp( cosh ) cosh mt K r r t t dt r t t dt                  (23) 25.001  mm tt (24) Numerical Computation of Cylindrical Conductor Internal Impedance 89 Simplification of expression (23) yields the following expression for scaled modified Bessel function of the second kind of order one: 1 1 1 0 0 ( ) exp( ) cos cosh exp( ) sin cosh m mt t s K r d d t dt j d d t dt              (25) As before the two integrals present in equation (25) are solved numerically using adaptive Simpson rule with high accuracy. Fig. 8 and Fig. 9 depict how the real and imaginary parts of equation (25) behave on the integration interval for various values of parameter α∙r. Fig. 8 Real part of scaled modified Bessel function of the second kind of order one for various values of parameter α∙r Fig. 9 Imaginary part of scaled modified Bessel function of the second kind of order one for various values of parameter α∙r 90 S. VUJEVIĆ, D. LOVRIĆ 3.4. Computation of scaled modified Bessel functions of the second kind for α∙r > 25 Asymptotic approximation of scaled modified Bessel functions of the second kind is given by the following expression [8]: 2 2 1 1 4 (2 1) ( ) ~ 1 ; 0, 1 2 ! (8 ) m s t n m m n t K r n r m r                                  (26) From the previous expression asymptotic approximations of scaled modified functions of the second kind of orders zero and one can be deduced: 0 1 ( 1) ( ) ~ 1 2 ( ) mNA s m m m c K r r r                 (27) 1 1 ( 1) ( ) ~ 1 2 ( ) mNA s m m m d K r r r                 (28) where NA is given by (16) whereas the coefficients cm and dm are computed from (17) and (18). 4. NUMERICAL EXAMPLES The presented model for computation of pul internal impedance of tubular conductors was implemented into a FORTRAN program. In order to ascertain the accuracy of obtained results and numerical stability of the model itself, a comparison is made with MATLAB which is used to compute pul internal impedance using the initial formula (1). Both FORTRAN and MATLAB employ double precision computing. It is important to note here that by using a program package which can employ more decimal places higher robustness of results would be achieved but at the expense of execution time. In the numerical example magnitudes and phase angles of Z for a thin tubular copper conductor (internal radius ri = 3.8 mm and external radius re = 4 mm) are computed. The results of the comparison are presented in Table 1 and Table 2. Table 1 Comparison of magnitudes of tubular cylindrical conductor internal impedance. α∙re Z (Ω) Proposed MATLAB 10 -2 0.003643657122067 0.003643657122067 10 -1 0.003643657122745 0.003643657122745 10 0 0.003643663902873 0.003643663902873 10 1 0.003710702668820 0.003710702668820 10 2 0.025181394368712 0.025181394368712 10 3 0.251267138203603 NaN 10 5 25.12049572965153 NaN 10 10 2512043.292911872 NaN 10 15 251204329284.9072 NaN Numerical Computation of Cylindrical Conductor Internal Impedance 91 Table 2 Comparison of phase angles of tubular cylindrical conductor internal impedance. α∙re φ (°) Proposed MATLAB 10 -2 9.30814638898·10 -6 9.30814663686·10 -6 10 -1 9.30814668336·10 -4 9.30814668521·10 -4 10 0 9.30813198101·10 -2 9.30813198100·10 -2 10 1 9.164530090507745 9.164530090507741 10 2 44.85885196305934 44.85885196305934 10 3 44.98566888986672 NaN 10 5 44.99985675983501 NaN 10 10 44.99999999856761 NaN 10 15 45.00000000000000 NaN As can be seen from the results in Table 1 and Table 2, when computing formula (1) using MATLAB an underflow/overflow stability issue occurs for larger function parameters. These numerical instabilities are a direct consequence of the denominator consisting of subtraction of two products. When these products become identical up to the last decimal place that the program package can compute, the denominator becomes equal to zero thus resulting in a Not a Number value. The proposed numerical procedure successfully circumvents these issues as can be seen form the results of the analysis. 5. CONCLUSION In this paper an algorithm for computation of pul internal impedance of cylindrical conductor under large complex function arguments is presented. The high accuracy and stability of the algorithm was achieved by selecting a formula for pul internal impedance which does not lead to undefined values for relatively small function arguments and by scaling the modified Bessel functions present in this formula by an appropriate scaling factor. The developed algorithm represents a basis for computation of pul internal impedance of multilayered tubular cylindrical conductors which is in development. REFERENCES [1] S. Vujević, D. Lovrić, "On the numerical computation of cylindrical conductor internal impedance for complex arguments of large magnitude", In Proceedings of the Extended Abstracts of the 12th International Conference on Applied Electromagnetics (ПЕС 2015), Niš, Serbia, 2015, pp. (P1_1) 1-4. [2] P. Sarajčev, S. Vujević, "Grounding Grid Analysis: Historical Background and Classification of Methods", International Review of Electrical Engineering, vol. 4, pp. 670-683, 2009. [3] H. W. Dommel, "EMTP Theory Book, 2nd edition", Microtran Power System Analysis Corporation, 1992. [4] F. P. Dawalibi, R. D. Southey, "Analysis of electrical interference from power lines to gas pipelines - Part I: Computation methods", IEEE Transactions on Power Delivery, vol. 4, no. 3, pp. 1840-1846, 1989. [5] J. Moore, R. Pizer, "Moment Methods in Electromagnetics - Techniques and Applications", John Wiley and Sons, 2007. [6] J. A. Stratton, "Electromagnetic theory", John Wiley & Sons, 2007. [7] S. Vujević, D. Lovrić, V. Boras, "High-Accurate Numerical Computation of Internal Impedance of Cylindrical Conductors for Complex Arguments of Arbitrary Magnitude", IEEE Transactions on Electromagnetic Compatibility, vol. 56, pp. 1431-1438, 2014. [8] M. Abramowitz, I. A. Stegun, "Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables", Dover Publications, 1964.