CC BY
21
Numerical Method for Calculating Non-stationary Processes in the Non-homogeneous Electric Circuit. Direct and Reverse Problem
Berzan V., Patsiuk V., Rybakova G.
Institute of Power Engineering of the Academy of Sciences of Moldova Chisinau, Republic of Moldova
Abstract. In the paper has been proposed the conservative numerical scheme of the calculation of dynamic processes in non-homogeneous electric circuits until reaching the phase of the stationary process. The numerical calculation scheme ensures the accuracy of the solution also in the case of the process analysis in the circuits with significant loss and dissipation of energy. The proposed method is also robust to solve the inverse problem in the field of mathematical physics, thus restoring the initial parameters of the non-stationary process based on the knowledge of the distribution of the voltage and current waves in the circuit. The results of the numerical solution were compared with those obtained by the finite difference of time method (FDTD) and the Godunov scheme. Calculations of the non-stationary process were performed in the partially homogeneous circuit with energy losses and high variability of linear parameters. The circuit under consideration is similar to the stator winding of a high power generator. It has been demonstrated the possibility of restoring the initial excitation parameters, for example, due to partial discharges or a short pulse. It has been found that the proposed numerical method can be used for purposes of increasing the precision of diagnosing the current state of insulation of high power rotating electric machines as a result of solving the inverse problem of the propagation of current and voltage waves in the non-homogeneous circuit.
Keywords: non-homogeneous circuits, distributed, variability of linear parameters, numerical method, losses, solves the telegraph equations, diagnostic.
Metoda numerica de calcul a proceselor nestationare in circuitul electric neomogen.
Problema directa inversa Berzan V., Patiuc V., Ribacova G.
Institutul de Energetica al Academiei de §tiinte a Moldovei Chi§inau, Republica Moldova Rezumat. in lucrare se examineaza metoda numerica de calculul a proceselor dinamice in circuite electrice neomogene pana la atingerea fazei procesului stationar. S-a propus schema numerica conservativa de calcul. Estimarea preciziei solutiilor numerice obtinute se bazeaza pe utilizarea legii de conservare a energiei in circuit. in schema de calcul propusa, se exclude efectul de acumulare al erorilor de calcul, inclusiv in modul de reflexie §i refractie multipla in circuitul neomogen. Metoda propusa este robusta §i pentru rezolvarea problemei inverse din domeniul fizicii matematice, deci restabilirea parametrilor initiali ai procesului nestationar in baza cunoa§terii repartitiei undelor de potential §i curent in circuitul analizat pentru oarecare moment de timp. Rezultatele solutiei numerice s-au comparat cu cele obtinute prin metoda finita a dferentei de timp (FDTD) §i schema lui Godunov. S-au executat calcule a procesului nestationar in circuitul partial omogen cu pierderi de energie §i variabilitate ridicata a parametrilor lineica. Circuitul examinat este similar infa§urarii statorului unui generator de mare putere. S-a demonstrat posibilitatea restabilirii parametrilor excitatiei initiale, de exemplu, conditionata de descarcarile partiale sau de un impuls de scurta durata. S-a constatat, ca metoda numerica propusa se poate utiliza in scopuri de sporire a preciziei de diagnosticare a starii curente a izolatiei ma§inilor electrice rotative de mare putere, ca urmare a rezolvarii problemei inverse a propagarii undelor de curent §i tensiune in circuitul neomogen.
Cuvinte-cheie: circuite neomogene, parametri distribuiti §i variabili, metoda numerica, pierderi, rezolvarea ecuatiilor telegrafice, diagnosticare.
Метод численного расчета нестационарных процессов в неоднородной электрической цепи.
Прямая и обратная задача Берзан В.П., Пацюк В.И., Рыбакова Г.А.
Институт энергетики Академии наук Молдовы Кишинэу, Республика Молдова Аннотация. В статье рассматривается численный метод расчета динамических процессов в неоднородной электрической цепи до установившегося значения. Была предложена консервативная численная схема расчета, которая учитывает явление диссипации и дисперсии энергии. Оценка точности
численных решений выполнена на основе закона сохранения энергии в цепи. Предлагаемая схема расчета, исключает эффект накопления ошибок вычислений, в том числе и в режиме многократного отражения волн в неоднородной цепи. Предложенный метод является работоспособным и при решении обратной задачи математической физики, таким образом можно восстанавливать исходные параметры первоначального процесса на основе знания распределения потенциала и волны тока в цепи в некоторый момент времени после начала нестационарного процесса в цепи. Результаты численного решения были сопоставлены с результатами, полученными с использованием метода конечных разностей во времени (РБТО) и схемой Годунова. Выполнены расчеты нестационарного процесса в частично однородной цепи с потерями мощности при высокой степени дискретного изменения значений параметров линии. Исследуемая неоднородная цепь похожа с цепью обмотки статора генератора большой мощности. Была показана возможность восстановления исходных параметров процесса (возбуждения), например, обусловленного частичным разрядом в изоляции или импульсом короткой длительности, приложенного к входу цепи. Было установлено, что предлагаемый численный метод может быть использован для повышения точности диагностики текущего состояния изоляции вращающейся электрической машины большой мощности по данным решения обратной задачи о распространении волн тока и напряжения в неоднородной цепи.
Ключевые слова: неоднородная цепь, распределенные и переменные погонные параметры, численный метод, потери, телеграфные уравнения, диагностика.
I. Introduction
Development of infrastructure of energy system is accompanied not only by new interrelations and physical structure complication, but by necessity of exact and efficient power management. The last one implies determination of active power value and energy losses during the stationary and dynamical regimes of its operation.
At present in the most of cases, it is accepted to consider stationary regime of electrical mains functioning as a normal mode. As an indicator of stationary regime, the constancy by time of physical values used for quantitative characteristics of energetic processes (electrical power transmission and distribution) is accepted. It is necessary to mention, that the stationary regime in its conventional interpretation does not exist physically in the ac lines, because all quantities are functions of time in such circuits. At present, standard method of stationary processes study in AC electrical circuits is based on the quasi-steady state determination, when the variable in time quantities are replaced by their energetic equivalents for alternative voltage and current period.
Since the operation of any energy system is not stationary by time in principle, all transient processes that can appear in energy system operation should be considered as normal modes. In this case (for this kind of problem definition), the role of theoretical calculations and mathematical modeling in all structural elements of energy system become of great importance. Since voltages and currents under transient processes differ from their time functions under
steady-state regimes, this make some additional difficulties in determination of energetic parameters and characteristics. For example, if we need to determine power losses up to hundredth part of percent, then the instantaneous values of voltages and currents at all parts of the circuit must be calculated accurate within four-five significant digits. This problem becomes still more complicated for circuits with distributed parameters: here we have to calculate not only distribution of instantaneous values of currents regarding to space variable, but we have to calculate wave processes with great accuracy taking into consideration multiple reflections and refractions of the incident wave caused by non-homogeneities of the circuit. Application of analytical methods in this case leads to many difficulties. Approximate numerical methods are more acceptable here, but their application requires accuracy rating of the obtained solutions. It is necessary to mention that in many cases such accuracy rating represents the independent scientific problem that not always can be solved identically and simply. It is natural to assume that the accuracy rating of numerical solutions can be realized by their comparison with some known analytical solutions for test problems.
Analysis of the corresponding publications [1]-[9], [27]-[32] shows clearly that it is not still formed strongly valid approach to calculation of the wave processes and the power transmission energy datum in the distributed systems with variable (tunable) parameters. Generation of cluster of test problems (with solutions known with high accuracy and recognized as standard sample solutions) is also one of important
research problems in the field of electrical engineering and power engineering.
By now, one can enumerate not so many nonstationary problems for electrical circuits with distributed and lumped parameters that have been solved [1]-[9]. In the strict sense, only the problem about the rectangular voltage and current wave motion along the homogeneous semiinfinite line is solved analytically. Unfortunately, only this analytical solution can be used as a sample solution for a posteriori accuracy estimation of the approximate methods. It has been proved that satisfaction accuracy coincides to minimum 2-3 significant digits [13]. In this case the numerical and the analytical solutions (represented in graphical form) are visually congruent.
In the majority of investigations some additional simplifying conditions have been used. For example, quite often during calculation of electrical circuits only active longitudinal resistance R > 0 is taking into account, but the transverse leakage (or shunt conductance) between the direct and the inverse wires is assumed to be equal to zero: G = 0 [1],[2],[4]-[8]. Even in so powerful commercial tool as the software complex SimPowerSystems [11], utilizing analog-digital models and methods EMTP-RV [12], parameter G is not used in mathematical model. As an example, as a rare exception, instantaneous connection to the direct voltage of the cable line with nonzero leakage current through the imperfect insulation has been investigated in [3].
In spite of this, the traditional approaches and methods for calculations of the loading regimes and commutation transient processes in distributing systems (high-current long circuits, communication lines, etc.) are sufficiently intricate and can not pretended to universality. More detailed review connected with the analytical or numerical solutions of the long line evolutionary equations is in [14]-[16].
For further development of theory and practice of electrical engineering and power engineering it is necessary currently to create some thesaurus (based on the latest measuring and computational technologies) containing numerical models and sample (test) examples for electromagnetic circuits and fields. Its following approbation on physical models should contribute to overcome of the existing gap between theoretical and experimental researches in the field of electrical power transmission at long distance [16],[17]. Solutions of such model
(test) problems should be represented in a maximally simple and convenient form. As a result, any specialist (familiar with the theoretical electrical and power engineering) will use it, and will repeat these results varying the initial data at his judgment. The proposal of some problems with their solutions for such cluster of sample problems represents one of the objectives of this research.
It is known, that none of the deductive methods of calculation or forecasting "does not like" the heavy gradients (neither by time nor by space). The situation becomes more complicated when it is necessary to calculate the shock wave evolution (strong discontinuities) in the distinctly non-homogeneous medium with parameters differing in orders. For example, wave resistances in the backbone power transmission lines and in the distributing networks with cable insertions differ in 8...12 times. If we consider the Franklin's lightning rod or Faraday cage as a piecewise homogeneous long line, then the linear leakage (shunt conductance) at the separate sections changes quite in hundreds of times. Under the emergency situations (such as open-phase fault and drop) the load resistance can suddenly go down from infinitely large values (at the idling regime) until zero (at short-circuited regime). The estimation of the limited parameters of the transient processes and elaboration of the effective approaches to their calculations remains the actual problem for reliability control of operation of the energy system.
Numerical method for calculation of non-homogeneous circuits with variable parameters has been proposed in this paper. Also, numerical results of its implementation for some typical electrical engineering problems have been presented. It is necessary to mention, that this numerical method is based on ideas of the known method of characteristics and of the method of the first differential approximation [14]. In particular, proposed method gives possibility to deduce the uniform computation equations for essentially non-homogeneous parametrical structures under the connection-disconnection of loads and of other lumped systems, and with different forms of input signals (continuous or discontinuous). Comparison of accuracy on some test problems solved by different methods (including analytical method, and well-known Finite Difference Time Domain (FDTD) method, and also by Baum-Liu-Tesche (BLT) method [27] - [32]), was been carried out.
II. Problem definition
For electrical circuit with distributed parameters both the load with lumped parameters and the voltage source can be connected to the arbitrary points of the circuit. The circuit with distributed parameters can be presented as piecewise homogeneous with big gradients of linear parameters alternations at the connection sections of long line. In this research electromagnetic energy transmission is determined by conduction currents and is described by telegraph equations. To provide unique solution, a set of telegraph equations should be completed by the corresponding boundary and initial conditions.
So, this research is aimed to prove and to test novel numerical method for calculation of electrical circuits with non-homogeneous structure under different types of losses (electrical circuit under this consideration contains subcircuits with distributed parameters, and chain loops with lumped parameters). Another purpose is to elaborate criterion for estimation of accuracy for numerical solutions, and to use numerical method for solution of some sample problems.
III. Energy integral and unicity theorem
Electromagnetic energy transmission along the long line by means of conduction currents can be described by telegraph equations that represent Kirchhoffs laws for closed circuit generated by subcircuit:
di du du di
L— + C — +Ri = 0; C— + G — + Gu = 0. (1) dt dx dt dx
cs
we
Obviously, when RS = LS = 0. have short circuited regime: u = 0 , and condition RS =cx> corresponds to idling regime of the line: i = 0 (the load is disconnected). Initial conditions usually are assumed to be equal to zero (electric charge is missing in circuit before commutation).
Now let consider situation of power take-off or connection of "bucking out" systems at the intermediate points x = xn of the line. In this case, the currents and voltages (as functions of spatial variable x) can experience discontinuities of the first kind or other jumps. However, (3) does not change its form if we substitute i = i\ - i2 and u = ui - u2, where the inferior indexes refer to function values at the left and at the right of the draw-off point. Let remark that the active-reactive lumped loads can consist from the arbitrary set of parallel and series connected RLC-circuits.
Starting from the general theory [18] let obtain energy integral for hyperbolic system (i). For this purpose it is necessary to multiply the first equation of (1) by i, and the second one - by u, and then to add the obtained results:
di .du 2 ^ du di 2 n Li — +i— + Ri2 + Cu — + u — + Gu2 = 0 or dt dx dt dx
1 d 2 2 2 2 d
--( Li2 + Cu 2) + Ri2 + Gu 2 +—(iu) = 0.
2 dt dx
Integrating the last expression over 0 < x < l, 0< x< t and taking into consideration the zero initial data, we obtain the following equation:
To obtain the unique solution, the (1) should be completed by boundary and initial conditions.
Let suppose that the electrical circuit at the initial time t = 0 is connected to the external voltage or current source:
u = U0 (t) or i = I0 (t), when x = 0 , (2) and its output is closed by active-reactive load in the form of series RLC-circuit for which the following integro-differential relation takes place:
di 1 '
u = RSi + LS~t + Ii(T)d1 , when x = l. (3)
t i , i
jj( Ri2 + Gu 2)dxd x + - j( Li2 +Cu 2)dx-
0 0 2 0 t
= j [ i(0, t) u (0, t) - i(l, t) u (l, t)] d x.
(4)
S 0
The left part of this energy balance equation represents sum of its active (irreversibly converting to the heat) and reactive (reversible) components, but the right part represents the difference between energies of the source and the receiver. All components in (4) have J (joule) dimension.
It is obvious that under zero initial and boundary conditions integral equation (4) holds only for trivial solution i = u = 0 . So, the unicity theorem results from Fredholm alternative in assumption of solution existence.
IV. Numerical method for telegraph equations in the case of line with losses
After this discretization FDA-based equations can be written as follows:
(L + ax)
(C + px)
( di_ dt
\
( du_ dt '
xd 2i >
2 dt
x d2u 2 "dt7
2
/
2. A
Ri-
du dx
- — + Gu-dx
hLa d i
2 dx2 hCa d2u
Hx2
= 0;
= 0.
A set of linear equations (1) is of hyperbolic type that implies the finite velocity of electromagnetic wave propagation along the line. Velocity is determined by linear parameters of the line as follows: a = uJlc . Godunov's finite-difference scheme [33], [34] of predictor-corrector type usually is constructed on the grids with integer and half-integer points [14],[33]:
( L + xa)n
in-1/2 -I
n n-1
u —u
+(Ri)n-1/2 =0;
un-i/2 -u in —in—1 (C+x^,)n-1/ 2u-^ + -+(Gu)n-1/ 2 =0, (5)
where
in = un-1/2 un+1/2 + (izb )n—1/2 + (izb )n+1/2 ;
(Zb ) 1/2 + (zb )n+1/2
Considering this fact regarding telegraph equations, the following second order equations can be obtained:
-Îr = LC + (LG + RC)— + RGi ;
dx2 dt2 dt
d 2u
d 2u
, = LC , dx2 dt2
duu
(LG + RC)— + RGu .
dt
And after that it is possible to transform FDA-based equations to the form:
_ di du di Lx d2i Lah d2i
L— + — + Ri + ax— +------- = 0 ;
dt dx dt 2 dt2 2 dx2
^du di ^ du C x d 2u Cah d 2u
C— + — + Gu + Bx— +--2----- = 0 .
dt dx dt 2 dt2 2 dx2
un = in-1/2 in+1/2 + (u / zb )n-1/2 + (u / zb )n+1/2 . Li —--
1/(Zb )n-1/2 + 1/(Zb )n+1/2
R GL G RC
a = —I--; B = —I--; x = hn 1/2 / an 1/9.
2 2^ 2 2L n-1/2 n-1/2
As it follows from the obtained equations, coefficients a and p can be chosen in order to minimize differential additions to the initial equations:
Wave resistance of the circuit is ZB =yj L / C , weighting coefficients a, p and the time step x are chosen in order to reduce to minimum finite difference dispersion and dissipation. For the better understanding it is possible to use First Differential Approximation (FDA) method for finite-difference equations (5). Firstly, it is necessary to transform these equations from general notation [18] to the form:
n r> hLa .
( L + ax)it +ux +Ri——ixx =0;
di Lx d2i
ax
dt Lah
dt
n du C x d 2u Bx— +---
dt 2 dt2
2 dt2
2
LC —i + ( LG + RC)— + RGi
dt
•0;
Cah
LC + ( LG + RC)— + RGu
dt2 dt
(C + px)u + ix +Gu- uxx = 0
So, after the corresponding transformations we obtain:
n
n
d2i ( Lx Lh
dt2 ^ 2 Lah
2a RGi ^ 0;
di dt
ax-—( LG + RC)
du ( C x Ch
dt2
2a
du ~ôt
Px-—( LG + RC )
Cah
RGu ^ 0.
Now we can observe, that when x = h a, values of coefficients should be defined from equations:
hLa
La2
a =-(LG + RC) =-(LG + RC) =
2x
1
=21R
LG C
hCa Ca2
P= ^ ( lg +RC) = — ( LG +RC) = 2x 2
= 1 ( G + — 21 L
and matrix of wave resistances is calculated as
z = l1/2c-1/2.
If we will denote by an1/2 the maximal velocity of the electromagnetic wave propagation and replace scalar values in formulas (5) by the corresponding vectors i = (i1, i2,..., im), u = (u1, u2,..., um) and matrixes L, C, R, G, Z, a , p, then we will obtain design equations for distributed system with arbitrary number of conductors.
V. Comparison of the FDTD method with Godunov's finite-difference scheme
It is possible to apply the FDTD method to the telegraph equations (1). In order to apply this method, at first (on the domain D = {(x, t): x e [0,l], t > 0}) we generate two grids with integer rahx and half-integer co hx nodes. Grid step h over the space variable is calculated as h = l / N, and the step over the time variable is chosen according to the scheme stability condition as x = h / a, where
a = 1/VLC is velocity of electromagnetic wave propagation. Thus, we have:
Under these conditions coefficients of time derivatives annul and the remaining terms tend to zero with first order when h ^ 0 .
Decomposition of elementary cells by space coordinate x is carried out in order to provide condition x = hn-m / an-1/2 = const for any n
index. Approximations for boundary conditions of (3), and for more general kinds of boundary conditions have been described in [14].
It is necessary to mention, that finite-difference equations (5) are taking into account linear parameter changes along the longitudinal coordinate x. So, they can be easily generalized in conformity with multiphase electrical systems with branch points and with other complicate factors.
In the case of multiwire line (when L, C are symmetrical square matrixes of self and mutual inductances and capacities), wave velocities correspond to the eigenvalues of the matrix
®hx = {(xm, tn): xm = mh; tn = nx;
hx v m ' n' m ' n '
A =
0 L-1 C-1 0
m = 0, N ; n = 0,1,2,...}
® hx = {(xm-1/2, tn-1/2) : xm-1/2 = Xm ~ h / 2;
tn-1/2 = tn - x/2; m = 0,N +1; n = 1,2,3,...} '
The main idea of the FDTD method is the follow: the current function i (x, t) is calculated at the integer nodes of the grid chx, but the voltage function u (x, t) is calculated at the halfinteger nodes of the grid cc hi. In this case, the derivatives from (1) can be approximated by finite differences with second order of accuracy with respect to parameters h and x . In this way, we obtain the following finite difference scheme:
■ n+1 -n n+1/2 n+1/2 -n+1 . -n
L m m + um+1/2-um-1/2 + R im +im = 0-
m = 0, N, n = 0,1,2,...
. n+1/2 n-1/2 ■n _-n , n+1/2 n-1/2
C "m-1/2 "m-1/2 + 'm 'm-1 + G Um_1/2 +"m-1/2 _ 0-
m _ 1,N, n _ 1,2,3,...
(6)
Equation (6) have to be completed by approximation of the initial and boundary conditions. In order to obtain the second order approximation for the initial condition we assume:
m _I(*m); "£/2 =|1_^T U(*m-1/2) _
G 2C
(7)
2Ch
[I((Xm ) - I(Xm-1)]
where I ( x) and U ( x) are values of current and voltage at the initial time moment t _ 0.
Boundary condition at the input of the line (2) has the following form:
n-1/2
n-1/2
+ 2U0 (tn-1/2 ) '
(8)
and condition (3) at the output (when LS _ 1/ CS _ 0) has been transformed as:
, n-1/2 N +1/2
1
, tR tRs
1+—+—-2L hL
2 RS'N-1 -f 1+Tl-TRS SN I 2L hL
(9)
n-1/2 N-1/2
Let compare accuracy of solution of (1)-(3) by Godunov's scheme, and by FDTD method (6)-(9) with the following values of the non-dimensional parameters: L = C = 1; R = G = 0.48; Rs= 3; l = 0.7. We consider that voltage and current are equal to zero at t = 0 ( U(x) = 0 ; I ( x) = 0 ), and the sinusoidal voltage at the input point of the line is equal to U0 (t) _ sin(2rct).
Calculations have been carried out up to the time moment t = Tmax = 4.
Table 1 contains results of comparison of exact analytical solution for current at the output point of the line Ia(t) with approximate solution I^(t), obtained by the FDTD method, and with approximate solution IG(t), obtained by Godunov's scheme. There have been used the
!f -1a ||C and
G -Ia ||C are
following notations: the maximal values of the differences between the solutions, -Ia\L and ||/G-Ia\L are the
corresponding mean square deviations. The first column of Table 1 contains numbers of grid nodes over space variable.
TABLE 1
The absolute error values of the solutions obtained by exact analytical method (ia), and by approximate methods: FDTD method (//), and Godunov's scheme (ig)
N PF - IM PG - IM PF - iJL2 Pg - IM2
10 0.011108 0.026568 0.013233 0.027565
20 0.002734 0.013365 0.003232 0.013599
40 0.000681 0.006638 0.000800 0.006736
80 0.000170 0.003308 0.000200 0.003357
160 0.000043 0.001651 0.000050 0.001674
Analysis of the obtained results illustrates clearly theoretical accuracy of these two methods: decreasing (by two times) of the grid step leads to four times decreasing of the FDTD method error, and to two times decreasing of the Godunov's scheme error (the second order of accuracy for the FDTD method and the first order - for the Godunov's scheme). In such a way, for continuous solutions of (1)-(3) the FDTD method is more exact and, correspondingly, is more preferable.
Now let consider the same problem, but with condition that during some period of time at the input of the line the short-circuit occurs, i.e. the input voltage becomes zero in some period of time: U0(t) = sin(2rct) when 0 < t < 0.6 or 2.6 < t < 4, and U0(t) = 0 when 0.6 < t < 2.6.
The corresponding time dependences of voltage and current at the input of the line (x = 0) for N = 10; 20; 40 are represented in Fig. 1. Exact analytical solution is marked by thick line, and solution obtained by the FDTD method is marked by thin line. Solution obtained by the means of Godunov's scheme practically does not differ from the exact solution.
*
i iii A c A
(\ Ml \ \
\ lllll III ,J Hi, i
\fi IflffM iV 1
'ill1 H II \ \
f I'l1' v v
T r
Fig. 1. Time dependence of voltage and current at the input of the line (a) when N = 10 (b); N = 20 (c), and N = 40 (d).
Fig. 1 demonstrates that the FDTD method of discontinuous solutions leads to the appearance of large oscillations which do not decrease with decreasing of the grid step. But there is no reason to be surprised since as early as 1959 it has been proved by Godunov [34]. It has been proved in [34], that among linear finite difference schemes with second order of accuracy for the equation du/dt + 8u/dx = 0 there is no one satisfying the condition of monotony, i.e. there is no one that does not lead to appearance of oscillations during computing of discontinuous solutions.
Thus, the FDTD method, in spite of the fact that it is of second order of accuracy, can be restrictedly applied to the telegraph equations as
it does not permit the modeling of such regimes as short-circuit and idling.
Regarding the BLT method [30]-[32], it is necessary to mention that it has restricted applications since the main assumption of its application is the presence of the steady-state sinusoidal regime in the long line. Thus, this method can not be used for calculation of transient and wave processes with multiple reflections of the incident (direct) wave.
VI. The solution of the direct and inverse problem by the method of finite differences.
The solution of the direct and inverse propagation of the voltage and current waves is of interest to the field of diagnosis of the current technical state of the electrical equipments, for example of the insulation of high power rotating electric machines. To illustrate the possibilities of the finite differences method, we present the results of simulation of the propagation and analysis of the waves of voltage and current in the stator phase winding of the synchronous generator with the 30 MVA power, which includes 72 portions with linear parameters:
L = 6• 10"7 H / m , R = 10~3 Ohm / m,
wave
iLC
C = 6.5• 10"11 F / m , G1 = 2.7 • 10~7Sm / m , the propagation velocity
= 1.6 • 108 m / 5 , the length of the
portion l1 = 1.67m, characteristic impedance
Z =
L2 = 5
L- = 96.1Ohm for the portions outside
C1
and respectively: C2 = 1.214 -10~9 F / m, • 10-5 Sm/m, = 2.7 m,
the stator notches
10~6 H / m , R20 = 0.6Ohm / m , a20 = 1.283-107 m / 5 ,
G2° = 3.8
z0 = 02 _
L
—0 = 964.2 Ohm for the winding
C2
portions located in the stator notches. At the output of the winding from the notch and at its entry into the notch, the linear parameters and the wave impedance values change their values by jumps (characteristic impedance with active resistance Z0 = R0). For our investigation it is
reasonable that the values of the analyzed circuit sizes is presented in the relative unit system [13],
considering that the analyzed circuit has the lengthl _ 167m. For these conditions we have:
L _ C _ a _ Z _1 , R _ 1.74 -10-3, G1 _ 4.33 -10-3,
l _ 10
-2 .
L2 _ '
1.33, C2 _ 18.68,
G2 _ 0.609,
l2 _ 1.62 -10
R2 _ 1.04, z02 _ 0.67.
Parameters marked with index 1 refer to the front portions and the parameters denoted by index 2 refer to the phase winding portions located in the stator notches.
A. The solution of the direct problem
We apply at the entrance of the non-homogeneous circuit a pulse with the half-wave shape, the duration of which coincides with the propagation time of this wave on the front and the notch of the winding. For these conditions one can write:
u(0,t) = 0 V t>ti ;
u(0,t)=sin(nt/tl) V t<t
where tx = (lj a1) + (l2/ a2) = 0.22 /us.
For the point x = 159 m that practically
coincides with the end of the circuit, the idle mode is achieved. The propagation time of the wave in the non-homogeneous phase winding is approx. 8 ^s. In order to ensure the stability of the numerical calculation scheme with the finite difference method, the restriction for the spacetime step of the calculation grid was formulated: t = hj/aj = h2/a2 . As a result of this condition, in the numerical calculation grid, the number of nodes on the notch portion exceeds the number of nodes of the grid placed on the front portion 20 times.
Based on some computational samples having the independent number of nodes on the front part between 3 and 30, it was found that the deviations between the values of the numerical solutions obtained are very small. As a result, the numerical calculation grid with 3-4 nodes on the front of the phase winding circuit was used for the following calculations.
In fig.2 it is shows the evolution of the non-stationary process at the input of the circuit (x = 0), at its middle (x = 79.5 m) and at the end of the circuit (x = l) for the during that is equal to the wave propagation time in the circuit (= 16 /us) In curves 1 and 2 describe of the current wave evolution in the time and space. Curve 3 represents the voltage wave. For a clearer illustration of the simulated process,
curves 2 (current) and 3 (voltage) are displaced with values of 0.5 and 1.0 relative units (see Figure 2), which excludes the overlap of the current and voltage curves in fig. 2. This provides a clearer view of the current curves (1) and voltage (2) over time.
Fig.2. The transient process for current (curves 1 and 2) and voltage (curve 3) in the TBC-30 generator phase winding circuit
In fig. 3 shows the evolution of the transient process in the non-homogeneously examined circuit with losses at the input of a "lightning pulse" signal. This signal is described in time by the following relationship:
u (0, t) = A(e-at - e-pt),
for which the constants values A=27,5, a = 18,79, /3 = 20,88 were selected from the condition that the energy of this impulse applied at the input of the phase winding of the TBC-30 generator coincides with the value of the sinusoidal impulse (Figure 2), which was equal to Wsm = 0.15u.r.
The results of the simulation of the transient process in the winding with losses at the application of the "lightning impulse" pulse are shown in fig.3.
Fig. 3. The transient process (current curves 1 and
2) and voltage (curve 3) in the phase winding circuit of the TBC-30 generator when applying of the "lightning impulse"
The character of the amplitude attenuation of the current waves (1) and of the voltage waves which appeared in the non-homogeneous circuit with losses due to the application of the "lightning impulse" is shown in fig. 4.
Fig. 4. Distribution of current waves (1) and voltage (2) in the circuit at time t = 16 ^s
B. The solution for the reverse problem
For the moment t = t0 > 0 we know the spatial distribution of the current and the voltage in the examined circuit. The boundary conditions at the ends of the circuit are also known V t e (0, t0). The problem lies in determining the solution for the domain 0<x<1 V t< t0 This is about getting the reverse solution to the problem. For this, in the calculus relations presented by the scheme in the finite differences it is necessary to assign to the variables x and t negative values, so the variable x is replaced by the variable -x and the variable t by the variable -t.
To approve this algorithm, we will executed in the following test calculations. We will examine two variants. The first variant refers to the circuit with distributed parameters with very low losses, which is close to an ideal circuit. Parameters of this circuit: L = C = 1,
R = 1.74 •10"3, G = 4.33 • 103. The second variant has circuit with energy losses. Losses are determined by the parameters that have the characteristic values for winding portion located in the stator's notch:
L = C = 1, R = 1.56,G = 4.33. Length of the circuit: l = 167m.
In fig. 5 is the presented of the numerical solution for the low loss circuit and in fig. 6 for the high loss circuit.
The impulses at the entry in the low loss circuit, including the return signal (Fig. 5), practically does not differ and have the values of
the amplitudes and the areas of the impulses very close curve.
1.0 0.8 0.6 0.4 0.2 0
' ' '
... 1;3 /12
/
/
0.5
1.0
1.5
2.0
x
Fig. 5. The currents (curves 1; 3) for the time t = 0.2 (1) and t = 1.8 (3-restored function), t = 2.0 (2) in the circuit with L=C=1, R=1.7410-3,
G=4.3310-
i 0.8
0.6
0.4
0.2
0 -0.1
0.5
1.0
1.5
2.0
x
Fig.6. The currents (curves 1; 3) for the time t =
0.2 (1) and t = 1.8 (3-restored function), t = 2.0 (2) in the circuit with L=C=1, R=1,56, G=0,408
In the case of the circuit with losses the amplitude has essentially decreased as a result of the energy dissipation in the circuit (fig. 6). The impulse noted (1) corresponds to the time point t = 0.2 .
In the non-homogeneous circuit, such as the stator phase of high power electric machines the process of impulse propagation is more complex (fig.7). This complexity is preserved in and the case of obtaining the solution for the inverse problem.
In order to obtain the solution of the inverse problem it is necessary to know the distribution of the current and voltage waves (initial conditions) at the time from which the restoration process of the inverse time is started.
In the case of verifying the correctness of the numerical solution for the inverse problem, it is natural to use as initial condition the numerical solution of the direct impulse propagation problem with known parameters in the circuit, e.g. for the time t. In fig. 6 shows the restored function of the current (curve 3), using the information obtained for time t = 2.0.
In fig. 7 is the presented of the numerical solution of direct problem of the propagation the pulse with the form of "lightning impulse" over the time interval 0 < t < 0.5 and the numerical solution of the inverse problem in the time interval 0.5 > t > 0.
0.25 0.5 0.25 t
Fig.7. The currents (curves 1,2,3) and the voltages (curves 4,5,6) in the points x = 0 (1; 4); 0.25 (2.5) and 0.47 (3; 6) for case of solution of the direct problem (0.5> t > 0) and for the time (0.5> t >0) we have the solution of inverse the problem
For estimating the accuracy of the reversed restoration process it is used as a criterion maintain the energy balance in the examined circuit during propagation (direct problem) and during the process of restoring the initial impulse (inverse problem). The results of the calculations are presentation in fig. 7 for the current and voltage waves and in fig. 8 the energy evolution in the magnetic field WL (2), the electric field WC (3), the sum of the energy accumulated in the magnetic and the electric field WLC=WL+WC (1), the dissipation of the energy in the active resistance of the conductor WR (5), the losses in the electrical insulation WG (6), the sum of the energy transformed in the heat WRG=WR+WG (4) for the direct problem and the inverse problem. The curve (7) shows the energy balance of the process in the circuit Ws= WL + WC + WR + WG, which is in good accordance with the law of conservation.
0.25 0.5 0.25 t
Fig.8. The evolution of the components (1-6) and the sum energy (7) for the direct solution and the inverse solution at the propagation of the "lightning pulse" signal in the stator phase of the generator TBC-30
The constant value of the circuit energy for both the numerical solution of the direct problem and the inverse problem is a convincing argument that the finite differences method allows to obtain the correct solution to simulate of the dynamic processes in the inversed time.
In the non-homogeneous circuits with high losses in which in the connection area of the portions with different values of the linear parameters there are jump changes of the values of these parameters, the precision of the restored functions of the current and voltage waves decreases, but precision remains satisfactory for solving different engineering problem, For example, to diagnose of the insulation of rotating electric machines.
The comparative analysis of the numerical solutions for the direct problem and the reverse problem confirms the satisfactory coincidence of the propagation signals in the homogeneous and non-homogeneous circuits obtained with the finite differences method. This is confirmed by relatively small deviations of the instantaneous values of the original signal and the restored signal as well as the symmetry of the energy function to the return point of the process (point of the entry pulse in the circuit). With the increase of losses in the circuit and the occurrence of the phenomenon of multiple refractions and of reflections of the waves in the jump points of the values of linear parameters of the circuit, the precision of the restoration of the original signal decreases.
However, this deviation is not very significant and will not strongly influence the results obtained in order to solve practical problems, for example the diagnosis of electrical insulation according to the intensity of partial discharges [35-37].
VII. Conclusion
1. Since any steady-state process is always preceded by transient process, their computations must be realized in the same consecution. Elaborated numerical scheme, named Albatross, is conservative with zero finite-difference dissipation and minimal dispersion. These properties result to the fact that computational error does not accumulate, and it gives possibility to provide transparent calculations of nonstationary solutions without loss of accuracy at large time intervals corresponding to 300...500 electromagnetic waves run along the line length right up to steady-state regime. At the
same time parameters of the line, of generator and of the load can change instantaneously modeling the load droppings or load pickups, the emergency situations as short-circuits, breaks in circuit, jump changes of the voltage, lightning strokes, etc.
2. Finite-difference-based scheme allows not only to repeat (with any order of accuracy) solutions obtained by classical methods, but it gives possibility to extend essentially class of solvable problems from theory of electrical circuits with variable parameters (in comparison with known methods).
3. It has been established and proved by numerical experiments that well-known FDTD method and BLT method elaborated for Maxwell equations are not applicable or give a poor accuracy in the case of discontinuous solutions in transient processes.
4. The results of investigation of the dynamic processes and interaction of the current and voltage pulses and wave in the lines with the losses, confirmated of the robustness of the finite differences method in the case of solving direct an inverse task in non-homogeneous circuits. The method qualitatively and quantitatively faithfully reflects the dynamics of the processes not only in the circuits with constant parameters, but also for the circuit the values of the parameters of which are functions of the spatial coordinate.
5. There was presented the method of solving the direct and inverse problems of impulses and waves propagation in non-homogeneous circuits with losses and the jumping of the linear parameters of the circuits. The method described can be used to diagnose the technical state of insulation of electrical and power equipment.
References
[1] K.A. Krug, Transient Processes in Linear Electrical Circuits, Moscow, 1948, 344 p. (in Russian).
[2] S. Haiasi, Waves in Electric Transmission Lines, Moscow-Leningrad, 1960, 343 p. (in Russian).
[3] Y.G. Kaganov, Electrical Circuits with Distributed Parameters and Chain Schemes, Moscow, 1990, 248 p. (in Russian).
[4] A.K .Mann and V.K. Spirodonov, "Wave method of determination of distance to the place of defect of the cable line," in Proc. 1959 of the VNIIE, Moscow, vol. 8, pp. 28-43 (in Russian).
[5] V.A. Vershkov, K.T. Nahapetean, O.V. Ol'shevschii, et al., "Complex tests of half-wave electric line at a 500 kV grid of the European part
of the USSR," Electricestvo (Electricity), 1968, no. 8, pp.10-16 (in Russian).
[6] M.V. Kostenko, M.I. Gumerova, A.N. Danilin, et al., Wave Processes and Overvoltage in Underground Grids, Leningrad, Energoatomizdat, 1991, 232 p. (in Russian).
[7] B.V. Efimov, Storm Overvoltage at Air Grids, Ed. of the Kolsky Branch of the Russian Academy of Sc., 2000, 134 p. (in Russian).
[8] G. Dragan, N. Golovanov, C. Mazzeti, et al., High Voltage Techno-logies, Vol. II, Bucure§ti, Ed. AGIR, 2001, 732 p. (in Romanian).
[9] J.B. Nitsch and S.V. Tkachenko, "Propagation of current waves along quasi-periodic thin-wire structures: taking radiation losses into account", Radio Science Bulletin, 2007, no. 322, pp. 1940.
[10] Neplan: by BCP Switzerland. Power Systems Engineering, www.neplan.com
[11] SimPowerSystems: www.mathworks.com/products/simulink
[12] EMTP-RV. www.emtp.com
[13] V. Berzan and V Rimschi, Nonstationary Processes in Inhomogeneous Electric Circuits, Chisinau, Ed. of the ASM, 1998, 412 p. (in Romanian).
[14] V. Rimschi, V. Berzan, V. Patiuc, et al., Wave phenomena at inhomogeneous structures, Vol. 5, Theory and methods of calculation of electric circuits, electromagnetic fields and protective shells, Chisinau, Ed. of the ASM, 2008, 664 p. (in Russian).
[15] M. Tirsu, V. Berzan, V. Rimschi, et al., "Research on influence of high-voltage cable un-homogeneities on process of short waves distribution," Electric Power Systems Research, 2008, vol. 78, no. 12, pp. 2046-2052.
[16] V. Rimschi, V. Berzan, and V. Patsyuk, "Theory and calculation of the line circuits with distributed and lumped parameters". Energy Technologies, 2009, vol. 42, no 2, pp. 12-22.
[17] K.S. Demircian and P.A. Butirin, "Problems of security and development of electrical and power branches in Russian Federation," Trans. of the Russian Academy of Sc., Power Engineering, 2008, no. 1, pp. 5-17 (in Russian).
[18] A.A. Samarskii, Basics of Numerical Methods, Moscow, Ed. Nauca, 1987, 465 p. (in Russian).
[19] V. Patsyuk, "Non-sine voltages and sine currents at the input of disconnected line with losses," Problems of the Regional Energetics, 2008, no. 3, pp. 47-60 (in Russian).
[20] G.N. Alexandrov, Transmission of the Electric Energy, Sankt Petersburg, Ed. of Technical University, 2007, 412 p. (in Russian).
[21] V.A. Soldatov and V.M. Postolati, Calculation and Optimization of Parameters and Regimes of Controlled Multiwires Lines, Chisinau, Ed. Shtiintsa, 1990, 239 p. (in Russian).
[22] E.V. Dmitriev, A.M. Gaçimov, A.K. Nair, et al., "Regarding calculation of parameters of the model of transmission line with consideration of corona of wires," Electricity (Elektrihcestvo),
2007, no. 12, pp. 2- 14 (in Russian).
[23] Yu.S. Beliakov."Regarding problem of identification of parameters of the air transmission line," Electricity (Electricestvo),
2008, no 6, pp. 18-23 (in Russian).
[24] RD153-34.3-35.125-99: Instructions regarding prevention of the 6-1150 kV electric lines from storms and internal overvoltages, Sankt Petersburg, Mintopenergo of RF, 1999, 355 p. (in Russian).
[25] E.M. Bazelian and Yu.P. Raizer, Physics of Storms and Preven-tion from Storms, Moscow, Ed. Fizmatlit, 2001, 320 p. (in Russian).
[26] V. Patsyuk, "Moving capacity discharge on the long transmission line," Problems of the Regional Energetics, 2008, no. 3, pp. 85-100 (in Russian).
[27] Kane Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Trans. on Antennas and Propagation, 2001, vol. 14, pp. 302-307.
[28] A Taflove and M.E. Brodwin, "Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell's equations," IEEE Trans. on Microwave Theory and Techniques, 2004, vol. 23. pp. 623-630.
[29] A. Taflove and M.E. Brodwin, "Computation of the electromagnetic fields and induced temperatures within a model of the microwave-irradiated human eye," IEEE Trans. on Microwave Theory and Techniques, 2004, vol. 23. p. 888-896.
[30] C.E. Baum, "The theory of electromagnetic interference control," Modern Radio Science, 1990, pp. 87-101.
[31] C.E.Baum, T.K. Liu, and F.M. Tesche, "On the analysis of general multiconductor transmissionline networks," Electromagnetic Topology for the Analysis and Design of Complex Systems. Fast Electrical and Optical Measurements, Ed. Martinus Nijhoff,1986, p.467-547.
[32] J.P. Parmantier, X. Ferrieres, S. Bertuol, and C.E. Baum, "Various ways to think of the resolution of the BLT equation with an LU technique," Interaction Notes, note 535, 1998.
[33] V. Patsyuk, Mathematical Methods for Electrical Circuits and Fi-elds, Chisinau, Ed. of the State University of Moldova, 2009, 442 p.
[34] S.K. Godunov, "Difference method of numerical calculation of discontinuous solution of hydrodynamic equations". Matematicheskii Sbornik, 1959, vol. 47 (89), no. 3, pp. 271-306 (in Russian).
[35] V.Berzan, M.Tar§u. Method of technical diagnosis of the winding insulation state of the electric machine stator. Patent MD 1383 G2 2000.09.30.
[36] V.Berzan, A. Barladeanu, M.Tar§u, F. Engster. Method for detecting the location of defects in insulation of stator slits of high power electric machines. Patent MD 1460 G2 of 2000.11.30.
[37] M.Tir§u, V.P., Berzan., V.X Rimschi., P. M. Postolache. Research on influence of highvoltage cable un-homogeneities on process of short waves distribution. Canada, ELECTRIC POWER SYSTEMS RESEARCH (EPSR), published by Elsevier (ISSN: 03787796), http://dx.doi.org/10.1016/j.epsr.2008.06. 011
About authors.
Berzan V.P. Doctor of Science, Deputy Director of the Power Engineering Institute of the ASM. Fields of scientific interest: energy, steady and transient processes in electrical circuits, mathematical modeling, diagnostics of energy equipment. E-
mail:[email protected]
Patsyuk V. Doctor of Physics and Mathematics, Associate Professor of the State University of Moldova, Leading Researcher of the Institute of Energy of the Academy of Sciences of Moldova. Areas of scientific interests: mathematical physics, numerical analysis, theoretical mechanics and theoretical electrical engineering. E-mail: [email protected]
Rybakova Galina. Doctor of Physics and Mathematics, Associate Professor of the State University of Moldova, Leading Researcher of the Institute of Energy of the Academy of Sciences of Moldova. Areas of scientific interests: mathematical physics, numerical analysis, mechanics of a deformed body. E-mail: [email protected]
